1 : C ALGOP1 SOURCE AM 12/02/28 21:15:06 7291 2 : Subroutine algoplas1(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,lnoconv,lcp,U) 6 : 7 : c SUBROUTINE added by THG 8 : 8 : c THIS SUBROUTINE IS CALLED BY DPRAN3D 9 : c the system is solved by the method of Broyden (Newton-Raphson secante) 10 : 11 : c rkt x function_t(x(1),x(2)) + (1-rkt) x(1) = 0 12 : c rkc x function_c(x(1),x(2)) + (1-rkc) x(2) = 0 13 : 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 : c trial stress vector, = elastic predictor 56 : **** dimension trav1(i1),trav2(i2),trav2b(i2),trav2c(i2), 57 : **** . trav6(i6),trav22(i2,i2) 58 : dimension trav1(1),trav2(2),trav2b(2),trav2c(2), 59 : . trav6(6),trav22(2,2) 60 : **** dimension P2(i6,i6),U(i6,i6),pi2(i6),v2loc(i2,i2) 61 : dimension P2(6,6),U(6,6),pi2(6),v2loc(2,2) 62 : dimension deps(i6) 63 : dimension parahot3(idimpara3) 64 : 65 : DATA itermax /90/, tolx /1.e-4/ 66 : 67 : ***** 68 : * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS) 69 : IF(I1.GT.1) PRINT *, ' ALGOPLAS1 - ERREUR I1 = ', I1, ' > 1 ' 70 : IF(I2.GT.2) PRINT *, ' ALGOPLAS1 - ERREUR I2 = ', I2, ' > 2 ' 71 : IF(I6.GT.6) PRINT *, ' ALGOPLAS1 - ERREUR I6 = ', I6, ' > 6 ' 72 : ***** 73 : 74 : 75 : i7 = 7 76 : i8 = 8 77 : i9 = 9 78 : i10=10 79 : i11=11 80 : i12=12 81 : i13=13 82 : i14=14 83 : i15=15 84 : i16=16 85 : i17=17 86 : i18=18 87 : i19=19 88 : i20=20 89 : i21=21 90 : i22=22 91 : i23=23 92 : i24=24 93 : i25=25 94 : i26=26 95 : i27=27 96 : i28=28 97 : i29=29 98 : i30=30 99 : i31=31 100 : i32=32 101 : i33=33 102 : i34=34 103 : i35=35 104 : i36=36 105 : i37=37 106 : i38=38 107 : i39=39 108 : 109 : r1 = 1. 110 : r2 = 2. 111 : r0 = 0. 112 : r10 = 10. 113 : i0 = 0 114 : i1 = 1 115 : 116 : lerror = .false. 117 : lproblem = .false. 118 : G = Eo/(r2*(r1+poison)) 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 calculate a first f[x(1),x(2),esigmae,rkappait,rkappaic] 130 : call find_f3d(f,esigmae,rkappait,rkappaic,x(i1),x(i2),alpha, 131 : . pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot, 132 : . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 133 : if (lerror) then 134 : return 135 : endif 136 : 137 : f(i1) = rkt*f(i1) 138 : f(i2) = rkc*f(i2) 139 : 140 : if (f(i1).gt.r0) lftact = .true. 141 : if (f(i2).gt.r0) lfcact = .true. 142 : 143 : if ((f(i1).lt.10.).and.(f(i2).lt.10.)) then 144 : return 145 : endif 146 : 147 : fn(i1) = f(i1) 148 : fn(i2) = f(i2) 149 : 150 : c ================================================================= 151 : c ! ! 152 : c ! initialization of Jacobian B ! 153 : c ! ============================ ! 154 : c ================================================================= 155 : 156 : c initialization of Jacobian B in x = (0;0) by the following equation 157 : c B(11) = - (dft/dsigma H66 dft/dsigma + dtaut) 158 : c B(12) = - (dft/dsigma H66 dfc/dsigma) 159 : c B(21) = - (dfc/dsigma H66 dft/dsigma) 160 : c B(22) = - (dfc/dsigma H66 dfc/dsigma + dtauc) 161 : 161 : call dsigRank(esigmae,dft,i3,i6,parahot3,idimpara3,rkappait,lcp) 162 : c ******** 163 : do iloc=i1,i6 164 : DFT(iloc) = rkt * DFT(iloc) 165 : enddo 166 : 166 : call dsigDrucker(fc,fb,esigmae,dfc,alpha,P2,pi2,i1,i6,lcp) 167 : c *********** 168 : do iloc=i1,i6 169 : DFC(iloc) = rkc * DFC(iloc) 170 : enddo 171 : 171 : call mulATB(dft,H66,trav6,i6,i1,i6,i6) 172 : c ****** 172 : call mulAB(trav6,dft,trav1,i1,i6,i6,i1) 173 : c ***** 174 : df11 = trav1(i1) 175 : 175 : call mulAB(trav6,dfc,trav1,i1,i6,i6,i1) 176 : c ***** 177 : df12 = trav1(i1) 178 : 179 : df21 = df12 180 : 180 : call mulATB(dfc,H66,trav6,i6,i1,i6,i6) 181 : c ****** 181 : call mulAB(trav6,dfc,trav1,i1,i6,i6,i1) 182 : c ***** 183 : df22 = trav1(i1) 184 : 185 : c derivatives of the hardening functions 186 : dtaut = dfhardct(rkappait) 187 : dtauc = dfhardcc(rkappaic,fc,fc0,rkappa1,ac,bc) 188 : 189 : B(i1,i1) = -(df11+rkt*dtaut) +(r1-rkt) 190 : B(i1,i2) = -df12 191 : B(i2,i1) = -df21 192 : B(i2,i2) = -(df22+rkc*dtauc) +(r1-rkc) 193 : 194 : c ================================================================= 195 : c ! ! 196 : c ! invert the initial JACOBIAN ! 197 : c ! =========================== ! 198 : c ================================================================= 199 : 200 : dtmB = B(i1,i1)*B(i2,i2)-B(i1,i2)*B(i2,i1) 201 : if (dtmB.eq.r0) then 202 : lerror = .true. 203 : lnoconv = .true. 204 : c print *,'ERROR BROYDEN dtmB=0' 205 : return 206 : else 207 : Binv(i1,i1) = B(i2,i2)/dtmB 208 : Binv(i1,i2) =-B(i1,i2)/dtmB 209 : Binv(i2,i1) =-B(i2,i1)/dtmB 210 : Binv(i2,i2) = B(i1,i1)/dtmB 211 : endif 212 : 213 : c ================================================================= 214 : c ! ! 215 : c ! loop "Newton-Raphson ! 216 : c ! ==================== ! 217 : c ================================================================= 218 : 219 : do iter=i1,itermax 220 : 221 : c ************************* 222 : c ! CALCULATION OF DX_i ! 223 : c ************************* 224 : c -1 225 : c DX_i = - J_i F(X_i) 226 : c 227 : TRAV22(i1,i1) = -BINV(i1,i1) 228 : TRAV22(i1,i2) = -BINV(i1,i2) 229 : TRAV22(i2,i1) = -BINV(i2,i1) 230 : TRAV22(i2,i2) = -BINV(i2,i2) 230 : call mulab(trav22,F,DX,i2,i2,i2,i1) 231 : if (dabs(dx(i1)) .lt. 1.E-16) dx(i1) = r0 232 : if (dabs(dx(i2)) .lt. 1.E-16) dx(i2) = r0 233 : if (((dabs(dx(i1)) .gt. 1.).and.(rkt.eq.1.)).or. 234 : . ((dabs(dx(i2)) .gt. 1.).and.(rkc.eq.1.))) then 235 : lnoconv = .true. 236 : lerror = .true. 237 : return 238 : endif 239 : if ( (dabs(dx(i1)) .gt. 1.).and.(rkt.eq.0.) ) 240 : . dx(1) = 1. * dx(1)/dabs(dx(i1)) 241 : if ( (dabs(dx(i2)) .gt. 1.).and.(rkc.eq.0.) ) 242 : . dx(2) = 1. * dx(2)/dabs(dx(i2)) 243 : 244 : c *************************** 245 : c ! CALCULATION OF X_i+1 ! 246 : c *************************** 247 : c 248 : c X_i+1 = X_i + DX_i 249 : c 250 : XNEW(i1) = X(i1) + DX(i1) 251 : XNEW(i2) = X(i2) + DX(i2) 252 : 253 : 254 : * PRINT *,'ALGOPLAS1 X(I1),DX(I1),XNEW(I1)',X(I1),DX(I1),XNEW(I1) 255 : * PRINT *,'ALGOPLAS1 X(I2),DX(I2),XNEW(I2)',X(I2),DX(I2),XNEW(I2) 256 : 257 : 258 : 259 : 260 : 261 : c The strain hardening parameters can only grow, they cannot decrease. 262 : 114 if ((xnew(i1).lt.r0).and.(xnew(i2).lt.r0)) then 263 : lnoconv = .true. 264 : return 265 : endif 266 : c From Simo et Hugues, Computational Inelasticity, p.212 procedure 2 267 : if (xnew(i1).lt.r0) then 268 : rkt = r0 269 : go to 87 270 : endif 271 : if (xnew(i2).lt.r0) then 272 : rkc = r0 273 : go to 87 274 : endif 275 : if ((xnew(i1).ge.r0).and.(lftact)) rkt=1. 276 : if ((xnew(i2).ge.r0).and.(lfcact)) rkc=1. 277 : 278 : rloc = xnew(i1) 279 : xnew(i1) = rkt*rloc 280 : c xnew(2) = dlambdat 281 : rloc = xnew(i2) 282 : xnew(i2) = rkc*rloc 283 : c xnew(2) = dlambdac 284 : 285 : c ********************************************************** 286 : c ! UPDATE OF THE INVERT JACOBIAN BY SHERMAN-MORRISON ! 287 : c ********************************************************** 288 : 289 : if ((dx(i1).ne.r0).or.(dx(i2).ne.r0)) then 290 : rloc = xnew(i1) 291 : xnew(i1) = rkt*rloc 292 : c xnew(1) = dlambdat 293 : rloc = xnew(i2) 294 : xnew(i2) = rkc*rloc 295 : c xnew(2) = dlambdac 296 : 297 : c calculate f[xnew(1),xnew(2),esigmae,rkappait,rkappaic] 297 : call find_f3d(fn,esigmae,rkappait,rkappaic,xnew(i1),xnew(i2), 298 : . alpha,pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot, 299 : . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 300 : 301 : if (lerror) then 302 : return 303 : endif 304 : 305 : c test activated yield surfaces 306 : if ((fn(i1).gt.r10).and.(lftact)) rkt = r1 307 : if ((fn(i2).gt.r10).and.(lfcact)) rkc = r1 308 : 309 : fn(i1) = rkt*fn(i1) 310 : fn(i2) = rkc*fn(i2) 311 : 312 : DF(i1) = FN(i1) - F(i1) 313 : DF(i2) = FN(i2) - F(i2) 314 : 315 : c calculation of the new Jacobian Binv_i+1 316 : c Now, we apply the formula of Sherman-Morrison 317 : c Binv.df 317 : call mulab(Binv,DF,TRAV2,i2,i2,i2,i1) 318 : c dx - Binv.df 319 : TRAV2b(i1) = DX(i1)-TRAV2(i1) 320 : TRAV2b(i2) = DX(i2)-TRAV2(i2) 321 : c dx.Binv 321 : call mulATB(dx,Binv,trav2c,i2,i1,i2,i2) 322 : c ****** 323 : c (dx-Binv.df).(dx.Binv) 323 : call mulAB(trav2b,trav2c,trav22,i2,i1,i1,i2) 324 : c ***** 325 : c dx.Binv.df 325 : call mulATB(dx,trav2,trav1,i2,i1,i2,i1) 326 : c ****** 327 : rloc = r1/trav1(i1) 328 : do jloc=i1,i2 329 : do iloc=i1,i2 330 : TRAV22(iloc,jloc) = rloc*TRAV22(iloc,jloc) 331 : v2loc(iloc,jloc) = BINV(iloc,jloc) 332 : BINV(iloc,jloc) = v2loc(iloc,jloc) + TRAV22(iloc,jloc) 333 : end do 334 : end do 335 : endif 336 : 337 : X(i1) = XNEW(i1) 338 : X(i2) = XNEW(i2) 339 : 340 : c ****************************** 341 : c ! TEST TO END THE LOOP ! 342 : c ****************************** 343 : 344 : test = r0 345 : if (iter.ne.i1) then 346 : c The plasticity functions must be satisfied with a precision of 10 Pa, i.e. 0.000010 N/m² 347 : if ((dabs(fn(i1)).le.10.).and.(dabs(fn(i2)).le.10.)) then 348 : c convergence has been obtained 349 : return 350 : endif 351 : 352 : do k=1,2 353 : if (x(k).ne.r0) then 354 : temp = abs(dx(k)/(1.E-10)) 355 : if (temp.gt.test) test = temp 356 : endif 357 : end do 358 : 359 : if (test.lt.tolx) then 360 : lnoconv = .true. 361 : return 362 : endif 363 : endif 364 : c END OF THE LOOP EITHER BECAUSE THE EQUATION IS SOLVED 365 : c OR BECAUSE THE X VARIABLE DOES NOT VARY ANY MORE 366 : c OTHERWISE, GO ON WITH THE ITERATIONS 367 : c ************************************ 368 : 369 : F(i1) = FN(i1) 370 : F(i2) = FN(i2) 371 : 372 : end do 373 : write(2,1) test,f(i1),f(i2),ntot 374 : 1 format( ' itermax exceeded in Broyden, test:',e7.2, 375 : . ', ft:',E8.2,' N/m², fc:',e8.2,' N/m²'/ 376 : . ' ntot = ',i6) 377 : lerror = .true. 378 : lnoconv = .true. 379 : 380 : RETURN 381 : END 382 : 383 : 384 :
© Cast3M 2003 - All rights reserved.
Disclaimer