1 : C ALGOP1CP SOURCE AM 12/02/28 21:15:00 7291 2 : Subroutine algoplas1cp(rkt,rkc,esigmae,H66,P2,pi2,alpha, 3 : . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot, 4 : . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6, 5 : . parahot3,idimpara3,deps,esigmaf,rkappat,rkappac, 6 : . lnoconv,lcp,U) 7 : 8 : c SUBROUTINE added by THG 9 : 9 : c THIS SUBROUTINE IS CALLED BY DPRAN3D 10 : c the system is solved by the method of Broyden (Newton-Raphson secante) 11 : 12 : c rkt x function_t(x(1),x(2)) + (1-rkt) x(1) = 0 13 : c rkc x function_c(x(1),x(2)) + (1-rkc) x(2) = 0 14 : 15 : c with 16 : c rkt and rkc = 0 or 1 17 : c x(1)=Dlambdat 18 : c x(2)=Dlambdac 19 : 20 : IMPLICIT REAL*8 (A-B,D-H,O-Z) 21 : implicit integer (I-K,M,N) 22 : implicit logical (L) 23 : implicit character*10 (C) 24 : 25 : **** dimension B(i2,i2) 26 : dimension B(2,2) 27 : c jacobian matrix for the system of equations. 28 : **** dimension Binv(i2,i2) 29 : dimension Binv(2,2) 30 : c invert of the jacobian matrix for the system of equations. 31 : **** dimension df(i2) 32 : dimension df(2) 33 : **** dimension dfc(i6) 34 : dimension dfc(6) 35 : c derivative of the compression plasticity function 36 : **** dimension dft(i6) 37 : dimension dft(6) 38 : c derivative of the tension plasticity function 39 : **** dimension dx(i2) 40 : dimension dx(2) 41 : c increment of the solution 42 : **** dimension f(i2) 43 : dimension f(2) 44 : **** dimension fn(i2) 45 : dimension fn(2) 46 : **** dimension H66(i6,i6) 47 : dimension H66(6,6) 48 : c matrix of elasticity (sometimes denoted D) 49 : **** dimension x(i2) 50 : dimension x(2) 51 : c solution to be found (and initial guess at the entry of the subr.) 52 : **** dimension xnew(i2) 53 : dimension xnew(2) 54 : **** dimension esigmae(i6) 55 : dimension esigmae(6) 56 : c trial stress vector, = elastic predictor 57 : **** dimension trav1(i1),trav2(i2),trav2b(i2),trav2c(i2), 58 : **** . trav6(i6),trav22(i2,i2) 59 : dimension trav1(1),trav2(2),trav2b(2),trav2c(2), 60 : . trav6(6),trav22(2,2) 61 : **** dimension P2(i6,i6),U(i6,i6),pi2(i6),v2loc(i2,i2) 62 : dimension P2(6,6),U(6,6),pi2(6),v2loc(2,2) 63 : **** dimension deps(i6),esigmain(6),dsigma(6),esigmaf(6) 64 : dimension deps(6),esigmain(6),dsigma(6),esigmaf(6) 65 : dimension parahot3(idimpara3) 66 : 67 : DATA itermax /90/, tolx /1.e-4/ 68 : 69 : ***** 70 : * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS) 71 : IF(I1.GT.1) PRINT *, ' ALGOPLAS1CP - ERREUR I1 = ', I1, ' > 1 ' 72 : IF(I2.GT.2) PRINT *, ' ALGOPLAS1CP - ERREUR I2 = ', I2, ' > 2 ' 73 : IF(I6.GT.6) PRINT *, ' ALGOPLAS1CP - ERREUR I6 = ', I6, ' > 6 ' 74 : ***** 75 : 76 : i7 = 7 77 : i8 = 8 78 : i9 = 9 79 : i10=10 80 : i11=11 81 : i12=12 82 : i13=13 83 : i14=14 84 : i15=15 85 : i16=16 86 : i17=17 87 : i18=18 88 : i19=19 89 : i20=20 90 : i21=21 91 : i22=22 92 : i23=23 93 : i24=24 94 : i25=25 95 : i26=26 96 : i27=27 97 : i28=28 98 : i29=29 99 : i30=30 100 : i31=31 101 : i32=32 102 : i33=33 103 : i34=34 104 : i35=35 105 : i36=36 106 : i37=37 107 : i38=38 108 : i39=39 109 : 110 : r1 = 1. 111 : r0 = 0. 112 : r2 = 2. 113 : r10 = 10. 114 : 115 : lerror = .false. 116 : lproblem = .false. 117 : G = Eo/(r2*(r1+poison)) 118 : lfirstiter = .true. 119 : 120 : 87 continue 121 : 122 : c initialization of fn 123 : FN(i1) = r0 124 : FN(i2) = r0 125 : 126 : c initialization of the activated yield surfaces 127 : lftact = .false. 128 : lfcact = .false. 129 : 130 : c initialization counter loop plane stress 131 : iderin = 0 132 : iter = 0 133 : 134 : c initialization 135 : do iloc=1,6 136 : esigmain(iloc) = esigmae(iloc) 137 : esigmaf(iloc) = esigmae(iloc) 138 : enddo 139 : rkappat = rkappait 140 : rkappac = rkappaic 141 : test = 0.d0 142 : dxprevt = 0.0d0 143 : dxprevc = 0.0d0 144 : 145 : c calculate a first f[x(1),x(2),esigmae,rkappait,rkappaic] 145 : call find_f3d(f,esigmae,rkappait,rkappaic,x(i1),x(i2),alpha, 146 : . pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot, 147 : . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 148 : if (lerror) then 149 : return 150 : endif 151 : 152 : if ((f(i1).lt.10.).and.(f(i2).lt.10.)) then 153 : return 154 : endif 155 : if (f(1).lt.10.) rkt = 0.d0 156 : 157 : f(i1) = rkt*f(i1) 158 : f(i2) = rkc*f(i2) 159 : 160 : if (f(i1).gt.r0) lftact = .true. 161 : if (f(i2).gt.r0) lfcact = .true. 162 : 163 : fn(i1) = f(i1) 164 : fn(i2) = f(i2) 165 : 166 : c ================================================================= 167 : c ! ! 168 : c ! loop "Newton-Raphson ! 169 : c ! ==================== ! 170 : c ================================================================= 171 : 172 : 11 iter=iter+i1 173 : 10 continue 174 : 175 : if (iter.lt.2) then 176 : 177 : c =================================================================== 178 : c ! ! 179 : c ! initialization of JACOBIAN B ! 180 : c ! ============================ ! 181 : c =================================================================== 182 : 183 : c initialization of Jacobian B in x = (0;0) by the following equation 184 : c B(11) = - (dft/dsigma H66 dft/dsigma + dtaut) 185 : c B(12) = - (dft/dsigma H66 dfc/dsigma) 186 : c B(21) = - (dfc/dsigma H66 dft/dsigma) 187 : c B(22) = - (dfc/dsigma H66 dfc/dsigma + dtauc) 188 : 188 : call dsigRank(esigmaf,dft,i3,i6,parahot3,idimpara3, 189 : . rkappait,lcp) 190 : c ******** 191 : do iloc=i1,i6 192 : DFT(iloc) = rkt * DFT(iloc) 193 : enddo 194 : 194 : call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2,pi2,i1,i6,lcp) 195 : c *********** 196 : do iloc=i1,i6 197 : DFC(iloc) = rkc * DFC(iloc) 198 : enddo 199 : 199 : call mulATB(dft,H66,trav6,i6,i1,i6,i6) 200 : c ****** 200 : call mulAB(trav6,dft,trav1,i1,i6,i6,i1) 201 : c ***** 202 : df11 = trav1(i1) 203 : 203 : call mulAB(trav6,dfc,trav1,i1,i6,i6,i1) 204 : c ***** 205 : df12 = trav1(i1) 206 : df21 = df12 207 : 207 : call mulATB(dfc,H66,trav6,i6,i1,i6,i6) 208 : c ****** 208 : call mulAB(trav6,dfc,trav1,i1,i6,i6,i1) 209 : c ***** 210 : df22 = trav1(i1) 211 : 212 : c derivatives of the hardening functions 213 : dtaut = dfhardct(rkappat) 214 : dtauc = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc) 215 : 216 : B(i1,i1) = -(df11+rkt*dtaut) +(r1-rkt) 217 : B(i1,i2) = -df12 218 : B(i2,i1) = -df21 219 : B(i2,i2) = -(df22+rkc*dtauc) +(r1-rkc) 220 : 221 : c =================================================================== 222 : c ! ! 223 : c ! invert the initial JACOBIAN ! 224 : c ! =========================== ! 225 : c =================================================================== 226 : 227 : dtmB = B(i1,i1)*B(i2,i2)-B(i1,i2)*B(i2,i1) 228 : if (dtmB.eq.r0) then 229 : lerror = .true. 230 : return 231 : else 232 : Binv(i1,i1) = B(i2,i2)/dtmB 233 : Binv(i1,i2) =-B(i1,i2)/dtmB 234 : Binv(i2,i1) =-B(i2,i1)/dtmB 235 : Binv(i2,i2) = B(i1,i1)/dtmB 236 : endif 237 : 238 : else 239 : 240 : c ********************************************************** 241 : c ! UPDATE THE INVERT JACOBIAN BY PAR SHERMAN-MORRISON ! 242 : c ********************************************************** 243 : 244 : c ! YIELD FUNCTIONS ! 245 : c 246 : c calculate f[0.d0,0.d0,esigmaf,rkappat,rkappac] 246 : call find_f3d(fn,esigmaf,rkappat,rkappac,0.d0,0.d0, 247 : . alpha,pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot, 248 : . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 249 : if (lerror) then 250 : return 251 : endif 252 : 253 : c Reactivation of yield surface 254 : if (fn(i1).gt.r10) rkt=1. 255 : if (fn(i2).gt.r10) rkc=1. 256 : 257 : fn(i1) = rkt*fn(i1) 258 : fn(i2) = rkc*fn(i2) 259 : 260 : rloc = xnew(i1) 261 : xnew(i1) = rkt*rloc 262 : c xnew(1) = dlambdat 263 : rloc = xnew(i2) 264 : xnew(i2) = rkc*rloc 265 : c xnew(2) = dlambdac 266 : 267 : DF(i1) = FN(i1) - F(i1) 268 : DF(i2) = FN(i2) - F(i2) 269 : 270 : c calculation of the new Jacobian Binv_i+1 271 : c Now, we apply the formula of Sherman-Morrison 272 : c Binv.df 272 : call mulab(Binv,DF,TRAV2,i2,i2,i2,i1) 273 : c dx - Binv.df 274 : TRAV2b(i1) = DX(i1)-TRAV2(i1) 275 : TRAV2b(i2) = DX(i2)-TRAV2(i2) 276 : c dx.Binv 276 : call mulATB(dx,Binv,trav2c,i2,i1,i2,i2) 277 : c ****** 278 : c (dx-Binv.df).(dx.Binv) 278 : call mulAB(trav2b,trav2c,trav22,i2,i1,i1,i2) 279 : c ***** 280 : c dx.Binv.df 280 : call mulATB(dx,trav2,trav1,i2,i1,i2,i1) 281 : c ****** 282 : rloc = r1/trav1(i1) 283 : do jloc=i1,i2 284 : do iloc=i1,i2 285 : TRAV22(iloc,jloc) = rloc*TRAV22(iloc,jloc) 286 : v2loc(iloc,jloc) = BINV(iloc,jloc) 287 : BINV(iloc,jloc) = v2loc(iloc,jloc) + TRAV22(iloc,jloc) 288 : end do 289 : end do 290 : 291 : X(i1) = XNEW(i1) 292 : X(i2) = XNEW(i2) 293 : 294 : F(i1) = FN(i1) 295 : F(i2) = FN(i2) 296 : 297 : endif 298 : 299 : c ****************************** 300 : c ! CALCULATION OF DX_i ! 301 : c ****************************** 302 : c -1 303 : c DX_i = - J_i F(X_i) 304 : c 305 : TRAV22(i1,i1) = -BINV(i1,i1) 306 : TRAV22(i1,i2) = -BINV(i1,i2) 307 : TRAV22(i2,i1) = -BINV(i2,i1) 308 : TRAV22(i2,i2) = -BINV(i2,i2) 308 : call mulab(trav22,F,DX,i2,i2,i2,i1) 309 : if (dabs(dx(i1)) .lt. 1.E-16) dx(i1) = r0 310 : if (dabs(dx(i2)) .lt. 1.E-16) dx(i2) = r0 311 : if (((dabs(dx(i1)) .gt. 1.).and.(rkt.eq.1.)).or. 312 : . ((dabs(dx(i2)) .gt. 1.).and.(rkc.eq.1.))) then 313 : lnoconv = .true. 314 : lerror = .true. 315 : return 316 : endif 317 : if ( (dabs(dx(i1)) .gt. 1.).and.(rkt.eq.0.) ) 318 : . dx(1) = 1. * dx(1)/dabs(dx(i1)) 319 : if ( (dabs(dx(i2)) .gt. 1.).and.(rkc.eq.0.) ) 320 : . dx(2) = 1. * dx(2)/dabs(dx(i2)) 321 : 322 : if (iter.gt.2) then 323 : if (dabs(dx(1)).gt.dabs((1.50d0 * dxprevt))) then 324 : rloc = dx(1) 325 : dx(1) = (1.50d0 * dxprevt) * rloc/(dabs(rloc)) 326 : endif 327 : if (dabs(dx(2)).gt.dabs((1.50d0 * dxprevc))) then 328 : rloc = dx(2) 329 : dx(2) = (1.50d0 * dxprevc) * rloc/(dabs(rloc)) 330 : endif 331 : endif 332 : 333 : dxprevt = dabs(dx(1)) 334 : dxprevc = dabs(dx(2)) 335 : 336 : c ***************************** 337 : c ! CALCULATION OF X_i+1 ! 338 : c ***************************** 339 : c 340 : c X_i+1 = X_i + DX_i 341 : c 342 : XNEW(i1) = X(i1) + DX(i1) 343 : XNEW(i2) = X(i2) + DX(i2) 344 : 345 : c The strain hardening parameters can only grow, they cannot decrease. 346 : 347 : if ((xnew(i1).le.r0).and.(xnew(i2).le.r0)) then 348 : lnoconv = .true. 349 : lerror = .true. 350 : return 351 : endif 352 : 353 : if (xnew(i1).eq.r0) then 354 : rkt = r0 355 : endif 356 : if (xnew(i2).eq.r0) then 357 : rkc = r0 358 : endif 359 : c new calculation of dx if we have reactivated a yield surface 360 : if ((xnew(i1).lt.r0).and.(rkc.eq.1.)) then 361 : rkt = r0 362 : if (lfirstiter) then 363 : x(1) = 0.d0 364 : x(2) = 0.d0 365 : lfirstiter = .false. 366 : go to 87 367 : else if (iter.lt.70) then 368 : go to 11 369 : else 370 : lerror = .true. 371 : lnoconv = .true. 372 : return 373 : endif 374 : endif 375 : if ((xnew(i2).lt.r0).and.(rkt.eq.1.)) then 376 : rkc = r0 377 : if (lfirstiter) then 378 : x(1) = 0.d0 379 : x(2) = 0.d0 380 : lfirstiter = .false. 381 : go to 87 382 : else if (iter.lt.70) then 383 : go to 11 384 : else 385 : lerror = .true. 386 : lnoconv = .true. 387 : return 388 : endif 389 : endif 390 : 391 : if (xnew(i1).gt.r0) rkt=1. 392 : if (xnew(i2).gt.r0) rkc=1. 393 : 394 : rloc = xnew(i1) 395 : xnew(i1) = rkt*rloc 396 : c xnew(1) = dlambdat 397 : rloc = xnew(i2) 398 : xnew(i2) = rkc*rloc 399 : c xnew(2) = dlambdac 400 : 401 : c dxcor(rected) for calculation of dsigma 402 : c dx = xnew-x because if xnew<0 => xnew=0 and the dx must be corrected 403 : dxcor1 = xnew(1) - x(1) 404 : dxcor2 = xnew(2) - x(2) 405 : 406 : c ************************** 407 : c ! CALCULATION OF DSIGMA ! 408 : c ************************** 409 : c 409 : call sigf3D_implicit(esigmain,dxcor1,dxcor2,esigmaf, 410 : . H66,P2,pi2,ag,i6, 411 : . fc,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 412 : if (lerror) then 413 : return 414 : endif 415 : do i=1,6 416 : dsigma(i) = esigmaf(i) - esigmain(i) 417 : enddo 418 : 419 : c *************************************** 420 : c ! PLANE STRESS ITERATIONS ! 421 : c *************************************** 422 : c 423 : bb = abs(dsigma(3)) 424 : esigmaf(3) = esigmain(3)-dsigma(3) 425 : esigmaf(1) = esigmain(1)-dsigma(3)*poison/(1.d0-poison) 426 : esigmaf(2) = esigmain(2)-dsigma(3)*poison/(1.d0-poison) 427 : if (bb.gt.(10.d0)) then 428 : if (iderin.eq.0) then 429 : iter = iter-1 430 : endif 431 : iderin = iderin+1 432 : if (iderin.gt.50) then 433 : write(6,*) ' problem internal iterations plane stress' 434 : return 435 : endif 436 : go to 10 437 : else 438 : continue 439 : endif 440 : 441 : c ************************************* 442 : c ! UPDATE OF THE STRESSES ! 443 : c ************************************* 444 : c 445 : do i=1,6 446 : esigmaf(i) = esigmaf(i) + dsigma(i) 447 : end do 448 : esigmaf(3) = 0.d0 449 : do i=1,6 450 : esigmain(i) = esigmaf(i) 451 : end do 452 : 453 : c ************************************************* 454 : c ! UPDATE OF THE HARDENING VARIABLES ! 455 : c ************************************************* 456 : X(i1) = XNEW(i1) 457 : X(i2) = XNEW(i2) 458 : 459 : c update of the hardening parameters 460 : rkappait= parahot3(idimpara3-37) 461 : rkappaic= parahot3(idimpara3-36) 462 : rkappat = rkappait+x(1) 463 : rkappac = rkappaic+x(2) 464 : 465 : c ************************* 466 : c ! YIELD FUNCTIONS ! 467 : c ************************* 468 : c 468 : call find_f3d(fn,esigmaf,rkappat,rkappac,0.d0,0.d0, 469 : . alpha,pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot, 470 : . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 471 : if (lerror) then 472 : return 473 : endif 474 : if (fn(i1).gt.r10) rkt=1. 475 : if (fn(i2).gt.r10) rkc=1. 476 : 477 : fn(i1) = rkt*fn(i1) 478 : fn(i2) = rkc*fn(i2) 479 : 480 : c ********************************** 481 : c ! TEST TO END THE LOOP ! 482 : c ********************************** 483 : if ((dabs(dx(i1)).lt.(1.E-10)).and.(dabs(dx(i2)).lt.(1.E-10)). 484 : . and.(abs(fn(i1)).le.10.).and.(abs(fn(i2)).le.10.)) then 485 : return 486 : endif 487 : test = r0 488 : if (iter.ne.i1) then 489 : c The plasticity functions must be satisfied with a precision of 10 Pa, i.e. 0.000010 N/m² 490 : if ((dabs(fn(i1)).le.10.).and.(dabs(fn(i2)).le.10.)) then 491 : c convergence has been obtained 492 : return 493 : endif 494 : 495 : do k=1,2 496 : if (dx(k).ne.r0) then 497 : temp = abs(dx(k)/(1.E-10)) 498 : if (temp.gt.test) test = temp 499 : endif 500 : end do 501 : if (test.lt.tolx) then 502 : return 503 : endif 504 : endif 505 : c END OF THE LOOP EITHER BECAUSE THE EQUATION IS SOLVED 506 : c OR BECAUSE THE X VARIABLE DOES NOT VARY ANY MORE 507 : c OTHERWISE, GO ON WITH THE ITERATIONS 508 : c ************************************ 509 : 510 : if (iter.lt.70) then 511 : go to 11 512 : else 513 : lerror = .true. 514 : lnoconv = .true. 515 : endif 516 : 517 : RETURN 518 : END 519 : 520 : 521 :
© Cast3M 2003 - All rights reserved.
Disclaimer