1 : C ALGOP2 SOURCE AM 12/02/28 21:15:08 7291 2 : Subroutine algoplas2(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 Newton 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 : c with 15 : c rkt and rkc = 0 or 1 16 : c x(1)=Dlambdat 17 : c x(2)=Dlambdac 18 : 19 : IMPLICIT REAL*8 (A-B,D-H,O-Z) 20 : implicit integer (I-K,M,N) 21 : implicit logical (L) 22 : implicit character*10 (C) 23 : 24 : **** dimension B(i2,i2) 25 : dimension B(2,2) 26 : c jacobian matrix for the system of equations. 27 : **** dimension Binv(i2,i2) 28 : dimension Binv(2,2) 29 : c invert of the jacobian matrix for the system of equations. 30 : **** dimension df(i2) 31 : dimension df(2) 32 : **** dimension dfc(i6) 33 : dimension dfc(6) 34 : c derivative of the compression plasticity function 35 : **** dimension dft(i6) 36 : dimension dft(6) 37 : c derivative of the tension plasticity function 38 : **** dimension dx(i2) 39 : dimension dx(2) 40 : c increment of the solution 41 : **** dimension f(i2) 42 : dimension f(2) 43 : **** dimension fn(i2) 44 : dimension fn(2) 45 : **** dimension H66(i6,i6) 46 : dimension H66(6,6) 47 : c matrix of elasticity (sometimes denoted D) 48 : **** dimension x(i2) 49 : dimension x(2) 50 : c solution to be found (and initial guess at the entry of the subr.) 51 : **** dimension xnew(i2) 52 : dimension xnew(2) 53 : **** dimension esigmae(i6) 54 : dimension esigmae(6) 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(6),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),rk(2),esigmaf(6),dsiglc(6) 64 : dimension B11(1),B12(1),B21(1),B22(1),xtest(2),dsiglt(6) 65 : 66 : DATA itermax /90/, tolx /1.e-4/ 67 : 68 : 69 : ***** 70 : * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS) 71 : IF(I1.GT.1) PRINT *, ' ALGOPLAS2 - ERREUR I1 = ', I1, ' > 1 ' 72 : IF(I2.GT.2) PRINT *, ' ALGOPLAS2 - ERREUR I2 = ', I2, ' > 2 ' 73 : IF(I6.GT.6) PRINT *, ' ALGOPLAS2 - ERREUR I6 = ', I6, ' > 6 ' 74 : ***** 75 : 76 : 77 : i7 = 7 78 : i8 = 8 79 : i9 = 9 80 : i10=10 81 : i11=11 82 : i12=12 83 : i13=13 84 : i14=14 85 : i15=15 86 : i16=16 87 : i17=17 88 : i18=18 89 : i19=19 90 : i20=20 91 : i21=21 92 : i22=22 93 : i23=23 94 : i24=24 95 : i25=25 96 : i26=26 97 : i27=27 98 : i28=28 99 : i29=29 100 : i30=30 101 : i31=31 102 : i32=32 103 : i33=33 104 : i34=34 105 : i35=35 106 : i36=36 107 : i37=37 108 : i38=38 109 : i39=39 110 : 111 : r1 = 1. 112 : r2 = 2. 113 : r0 = 0. 114 : r10 = 10. 115 : i0 = 0 116 : i1 = 1 117 : 118 : lerror = .false. 119 : lproblem = .false. 120 : G = Eo/(r2*(r1+poison)) 121 : rk(1) = rkt 122 : rk(2) = rkc 123 : 124 : 87 continue 125 : 126 : c initialization of fn 127 : FN(i1) = r0 128 : FN(i2) = r0 129 : 130 : c initialization 131 : do iloc=1,6 132 : esigmaf(iloc) = esigmae(iloc) 133 : enddo 134 : rkappat = rkappait 135 : rkappac = rkappaic 136 : test = 0.d0 137 : xnew(1) = 0.0d0 138 : xnew(2) = 0.0d0 139 : 140 : c calculate a first f[x(1),x(2),esigmae,rkappait,rkappaic] 140 : call find_f3d(f,esigmae,rkappait,rkappaic,x(i1),x(i2),alpha, 141 : . pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot, 142 : . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 143 : if (lerror) then 144 : return 145 : endif 146 : 147 : if ((f(i1).lt.10.).and.(f(i2).lt.10.)) then 148 : return 149 : endif 150 : 151 : if (f(1).lt.10.) rk(1) = 0.d0 152 : 153 : fn(i1) = f(i1) 154 : fn(i2) = f(i2) 155 : 156 : c ============================================================= 157 : c ! ! 158 : c ! initialization of JACOBIAN B ! 159 : c ! ============================ ! 160 : c ============================================================= 161 : 162 : c initialization of the Jacobian B by the following equation 163 : c B(11) = - (dft/dsigma dsigma/dlambdat - dtaut) 164 : c B(12) = - (rk(1)*rk(2)*dft/dsigma dsigma/dlambdac) 165 : c B(21) = - (rk(2)*rk(1)*dfc/dsigma dsigmaLdlambdat) 166 : c B(22) = - (dfc/dsigma dsigmaLdlambdac - (1-alpha)*dtauc) 167 : 167 : call dsigRank(esigmaf,dft,i3,i6,parahot3,idimpara3, 168 : . rkappait,lcp) 168 : call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2,pi2,i1,i6,lcp) 169 : 169 : call dsiglambda(xnew(1),xnew(2),P2,pi2,ag,fc,fb,lcp,rkappait, 170 : . parahot3,idimpara3,U,H66,esigmaf,dsiglt,dsiglc) 171 : 171 : call mulATB(dft,dsiglt,B11,6,1,6,1) 171 : call mulATB(dft,dsiglc,B12,6,1,6,1) 171 : call mulATB(dfc,dsiglt,B21,6,1,6,1) 171 : call mulATB(dfc,dsiglc,B22,6,1,6,1) 172 : 173 : dtaut = dfhardct(rkappat) 174 : dtauc = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc) 175 : 176 : B(1,1) = B11(1) - dtaut 177 : B(1,2) = rk(1) * rk(2) * B12(1) 178 : B(2,1) = rk(2) * rk(1) * B21(1) 179 : B(2,2) = B22(1) - (1.-alpha) * dtauc 180 : 181 : dtmB = B(i1,i1)*B(i2,i2)-B(i1,i2)*B(i2,i1) 182 : if (dtmB.eq.r0) then 183 : lerror = .true. 184 : print *,'ERROR BROYDEN dtmB=0' 185 : return 186 : else 187 : Binv(i1,i1) = B(i2,i2)/dtmB 188 : Binv(i1,i2) =-B(i1,i2)/dtmB 189 : Binv(i2,i1) =-B(i2,i1)/dtmB 190 : Binv(i2,i2) = B(i1,i1)/dtmB 191 : endif 192 : 193 : c ================================================================= 194 : c ! ! 195 : c ! loop "Newton-Raphson ! 196 : c ! ==================== ! 197 : c ================================================================= 198 : 199 : do iter=i1,itermax 200 : 201 : c ************************** 202 : c ! CALCULATION OF DE DX_i ! 203 : c ************************** 204 : c -1 205 : c DX_i = - J_i F(X_i) 206 : c 207 : TRAV22(i1,i1) = -BINV(i1,i1) 208 : TRAV22(i1,i2) = -BINV(i1,i2) 209 : TRAV22(i2,i1) = -BINV(i2,i1) 210 : TRAV22(i2,i2) = -BINV(i2,i2) 210 : call mulab(trav22,F,DX,i2,i2,i2,i1) 211 : if (dabs(dx(i1)) .lt. 1.E-18) dx(i1) = r0 212 : if (dabs(dx(i2)) .lt. 1.E-18) dx(i2) = r0 213 : if (((dabs(dx(i1)) .gt. 1.).and.(rkt.eq.1.)).or. 214 : . ((dabs(dx(i2)) .gt. 1.).and.(rkc.eq.1.))) then 215 : lnoconv = .true. 216 : lerror = .true. 217 : return 218 : endif 219 : if ((dabs(dx(i1)).gt.1.).and.(rkt.eq.0.)) 220 : . dx(1) = 1. * dx(1)/dabs(dx(i1)) 221 : if ((dabs(dx(i2)).gt.1.).and.(rkc.eq.0.)) 222 : . dx(2) = 1. * dx(2)/dabs(dx(i2)) 223 : 224 : c ************************* 225 : c ! CALCULATION OF X_i+1 ! 226 : c ************************* 227 : c 228 : c X_i+1 = X_i + DX_i 229 : c 230 : do iloc=1,2 231 : Xtest(iloc) = X(iloc) + DX(iloc) 232 : if ((Xtest(iloc).le.(0.)).and.(rk(iloc).eq.1)) then 233 : do jloc=1,2 234 : if ((rk(jloc).eq.1).and.(jloc.ne.iloc)) 235 : . DX(jloc) = - f(jloc)*BINV(jloc,jloc) 236 : enddo 237 : endif 238 : enddo 239 : 240 : XNEW(i1) = X(i1) + DX(i1) 241 : XNEW(i2) = X(i2) + DX(i2) 242 : XNEW(i1) = max(0.D0,XNEW(i1)) 243 : XNEW(i2) = max(0.D0,XNEW(i2)) 244 : 245 : do iloc=1,2 246 : if (XNEW(iloc).eq.(0.)) then 247 : rk(iloc) = 0. 248 : else 249 : rk(iloc) = 1. 250 : endif 251 : enddo 252 : 253 : c The strain hardening parameters can only grow, they cannot decrease. 254 : if ((xnew(i1).le.r0).and.(xnew(i2).le.r0)) then 255 : lnoconv = .true. 256 : lerror = .true. 257 : return 258 : endif 259 : 260 : c ********************************************** 261 : c ! UPDATE OF THE HARDENING PARAMETERS ! 262 : c ********************************************** 263 : 264 : X(i1) = XNEW(i1) 265 : X(i2) = XNEW(i2) 266 : c update of the hardening parameters 267 : rkappait= parahot3(idimpara3-37) 268 : rkappaic= parahot3(idimpara3-36) 269 : rkappat = rkappait+x(1) 270 : rkappac = rkappaic+x(2) 271 : 272 : c ****************************** 273 : c ! UPDATE OF THE STRESSES ! 274 : c ****************************** 275 : 275 : call sigf3D_implicit(esigmae,xnew(1),xnew(2),esigmaf,H66,P2, 276 : . pi2,ag,i6,fc,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 277 : if (lerror) then 278 : return 279 : endif 280 : 281 : c ************************** 282 : c ! YIELD FUNCTIONS ! 283 : c ************************** 284 : c 285 : c calculate f[0.d0,0.d0,esigmaf,rkappat,rkappac] 285 : call find_f3d(fn,esigmaf,rkappat,rkappac,0.d0,0.d0, 286 : . alpha,pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot, 287 : . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 288 : if (lerror) then 289 : return 290 : endif 291 : 292 : if (fn(i1).gt.r10) rk(1)=1. 293 : if (fn(i2).gt.r10) rk(2)=1. 294 : 295 : rkt = rk(1) 296 : rkc = rk(2) 297 : 298 : F(i1) = FN(i1) 299 : F(i2) = FN(i2) 300 : 301 : c ************************************* 302 : c ! UPDATE OF THE INVERT JACOBIAN ! 303 : c ************************************* 304 : 304 : call dsigRank(esigmaf,dft,i3,i6,parahot3,idimpara3, 305 : . rkappait,lcp) 305 : call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2,pi2,i1,i6,lcp) 306 : 306 : call dsiglambda(xnew(1),xnew(2),P2,pi2,ag,fc,fb,lcp,rkappait, 307 : . parahot3,idimpara3,U,H66,esigmaf,dsiglt,dsiglc) 308 : 308 : call mulATB(dft,dsiglt,B11,6,1,6,1) 308 : call mulATB(dft,dsiglc,B12,6,1,6,1) 308 : call mulATB(dfc,dsiglt,B21,6,1,6,1) 308 : call mulATB(dfc,dsiglc,B22,6,1,6,1) 309 : 310 : dtaut = dfhardct(rkappat) 311 : dtauc = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc) 312 : 313 : B(1,1) = B11(1) - dtaut 314 : B(1,2) = rk(1) * rk(2) * B12(1) 315 : B(2,1) = rk(2) * rk(1) * B21(1) 316 : B(2,2) = B22(1) - (1.-alpha) * dtauc 317 : 318 : dtmB = B(i1,i1)*B(i2,i2)-B(i1,i2)*B(i2,i1) 319 : if (dtmB.eq.r0) then 320 : lerror = .true. 321 : print *,'ERROR BROYDEN dtmB=0' 322 : return 323 : else 324 : Binv(i1,i1) = B(i2,i2)/dtmB 325 : Binv(i1,i2) =-B(i1,i2)/dtmB 326 : Binv(i2,i1) =-B(i2,i1)/dtmB 327 : Binv(i2,i2) = B(i1,i1)/dtmB 328 : endif 329 : 330 : c ****************************** 331 : c ! TEST END OF THE LOOP ! 332 : c ****************************** 333 : 334 : test = r0 335 : if (iter.ne.i1) then 336 : c The plasticity functions must be satisfied with a precision of 10 Pa, i.e. 0.000010 N/m² 337 : if ((dabs(rk(1)*fn(i1)).le.10.).and. 338 : . (dabs(rk(2)*fn(i2)).le.10.)) then 339 : c convergence has been obtained 340 : return 341 : endif 342 : 343 : do k=1,2 344 : if (dx(k).ne.r0) then 345 : temp = abs(rk(k)*dx(k)/(1.E-10)) 346 : if (temp.gt.test) test = temp 347 : endif 348 : end do 349 : if (test.lt.tolx) then 350 : if ((dabs(rk(1)*fn(i1)).le.100.).and. 351 : . (dabs(rk(2)*fn(i2)).le.100.)) then 352 : return 353 : else 354 : lnoconv = .true. 355 : return 356 : endif 357 : endif 358 : endif 359 : c END OF THE LOOP EITHER BECAUSE THE EQUATION IS SOLVED 360 : c OR BECAUSE THE X VARIABLE DOES NOT VARY ANY MORE 361 : c OTHERWISE, GO ON WITH THE ITERATIONS 362 : c ************************************ 363 : 364 : F(i1) = FN(i1) 365 : F(i2) = FN(i2) 366 : 367 : end do 368 : write(2,1) test,f(i1),f(i2),ntot 369 : 1 format( ' itermax exceeded in Broyden, test:',e7.2, 370 : . ', ft:',E8.2,' N/m², fc:',e8.2,' N/m²'/ 371 : . ' ntot = ',i6) 372 : lerror = .true. 373 : lnoconv = .true. 374 : 375 : RETURN 376 : END 377 : 378 : 379 :
© Cast3M 2003 - All rights reserved.
Disclaimer