1 : C ALGOP2CP SOURCE AM 12/02/28 21:15:07 7291 2 : Subroutine Algoplas2cp(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 Newton 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 : 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(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),esigmain(6),dsigma(6),esigmaf(6) 63 : dimension parahot3(idimpara3),rk(2),xtest(2),dsiglt(6),dsiglc(6) 64 : dimension B11(1),B12(1),B21(1),B22(1) 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 *, ' ALGOPLAS2CP - ERREUR I1 = ', I1, ' > 1 ' 72 : IF(I2.GT.2) PRINT *, ' ALGOPLAS2CP - ERREUR I2 = ', I2, ' > 2 ' 73 : IF(I6.GT.6) PRINT *, ' ALGOPLAS2CP - 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 : r0 = 0. 113 : r2 = 2. 114 : r10 = 10. 115 : 116 : lerror = .false. 117 : lproblem = .false. 118 : lonechance = .true. 119 : G = Eo/(r2*(r1+poison)) 120 : rk(1) = rkt 121 : rk(2) = rkc 122 : 123 : c initialization of fn 124 : FN(i1) = r0 125 : FN(i2) = r0 126 : 127 : c initialization counter loop plane stress 128 : iderin = 0 129 : iter = 0 130 : 131 : c initialization 132 : do iloc=1,6 133 : esigmain(iloc) = esigmae(iloc) 134 : esigmaf(iloc) = esigmae(iloc) 135 : enddo 136 : rkappat = rkappait 137 : rkappac = rkappaic 138 : test = 0.d0 139 : xnew(1) = 0.0d0 140 : xnew(2) = 0.0d0 141 : 142 : c calculate a first f[x(1),x(2),esigmae,rkappait,rkappaic] 142 : call find_f3d(f,esigmae,rkappait,rkappaic,x(i1),x(i2),alpha, 143 : . pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot, 144 : . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 145 : if (lerror) then 146 : return 147 : endif 148 : 149 : if ((f(i1).lt.10.).and.(f(i2).lt.10.)) then 150 : return 151 : endif 152 : if (f(1).lt.10.) rk(1) = 0.d0 153 : 154 : fn(i1) = f(i1) 155 : fn(i2) = f(i2) 156 : 157 : c ================================================================= 158 : c ! ! 159 : c ! loop "Newton-Raphson ! 160 : c ! ==================== ! 161 : c ================================================================= 162 : 163 : 11 iter=iter+i1 164 : 10 continue 165 : 166 : 167 : c ************************************* 168 : c ! CALCULATION OF THE JACOBIAN ! 169 : c ************************************* 170 : 170 : call dsigRank(esigmaf,dft,i3,i6,parahot3,idimpara3, 171 : . rkappait,lcp) 171 : call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2,pi2,i1,i6,lcp) 172 : 172 : call dsiglambda(xnew(1),xnew(2),P2,pi2,ag,fc,fb,lcp,rkappait, 173 : . parahot3,idimpara3,U,H66,esigmaf,dsiglt,dsiglc) 174 : 174 : call mulATB(dft,dsiglt,B11,6,1,6,1) 174 : call mulATB(dft,dsiglc,B12,6,1,6,1) 174 : call mulATB(dfc,dsiglt,B21,6,1,6,1) 174 : call mulATB(dfc,dsiglc,B22,6,1,6,1) 175 : 176 : dtaut = dfhardct(rkappat) 177 : dtauc = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc) 178 : 179 : B(1,1) = B11(1) - dtaut 180 : B(1,2) = rk(1) * rk(2) * B12(1) 181 : B(2,1) = rk(2) * rk(1) * B21(1) 182 : B(2,2) = B22(1) - (1.-alpha) * dtauc 183 : 184 : dtmB = B(i1,i1)*B(i2,i2)-B(i1,i2)*B(i2,i1) 185 : if (dtmB.eq.r0) then 186 : lerror = .true. 187 : return 188 : else 189 : Binv(i1,i1) = B(i2,i2)/dtmB 190 : Binv(i1,i2) =-B(i1,i2)/dtmB 191 : Binv(i2,i1) =-B(i2,i1)/dtmB 192 : Binv(i2,i2) = B(i1,i1)/dtmB 193 : endif 194 : 195 : c ************************* 196 : c ! CALCULATION OF DX_i ! 197 : c ************************* 198 : c -1 199 : c DX_i = - J_i F(X_i) 200 : c 201 : TRAV22(i1,i1) = -BINV(i1,i1) 202 : TRAV22(i1,i2) = -BINV(i1,i2) 203 : TRAV22(i2,i1) = -BINV(i2,i1) 204 : TRAV22(i2,i2) = -BINV(i2,i2) 204 : call mulab(trav22,F,DX,i2,i2,i2,i1) 205 : if (dabs(dx(i1)) .lt. 1.E-18) dx(i1) = r0 206 : if (dabs(dx(i2)) .lt. 1.E-18) dx(i2) = r0 207 : if (((dabs(dx(i1)) .gt. 1.).and.(rk(1).eq.1.)).or. 208 : . ((dabs(dx(i2)) .gt. 1.).and.(rk(2).eq.1.))) then 209 : lnoconv = .true. 210 : lerror = .true. 211 : return 212 : endif 213 : if ( (dabs(dx(i1)) .gt. 1.).and.(rk(1).eq.0.) ) 214 : . dx(1) = 1. * dx(1)/dabs(dx(i1)) 215 : if ( (dabs(dx(i2)) .gt. 1.).and.(rk(2).eq.0.) ) 216 : . dx(2) = 1. * dx(2)/dabs(dx(i2)) 217 : 218 : c ************************* 219 : c ! CALCULATION OF X_i+1 ! 220 : c ************************* 221 : c 222 : c X_i+1 = X_i + DX_i 223 : c 224 : do iloc=1,2 225 : Xtest(iloc) = X(iloc) + DX(iloc) 226 : if ((Xtest(iloc).le.(0.)).and.(rk(iloc).eq.1)) then 227 : do jloc=1,2 228 : if ((rk(jloc).eq.1).and.(jloc.ne.iloc)) 229 : . DX(jloc) = - f(jloc)*BINV(jloc,jloc) 230 : enddo 231 : endif 232 : enddo 233 : 234 : XNEW(i1) = X(i1) + DX(i1) 235 : XNEW(i2) = X(i2) + DX(i2) 236 : 237 : XNEW(i1) = max(0.D0,XNEW(i1)) 238 : XNEW(i2) = max(0.D0,XNEW(i2)) 239 : 240 : do iloc=1,2 241 : if (XNEW(iloc).eq.(0.)) then 242 : rk(iloc) = 0. 243 : else 244 : rk(iloc) = 1. 245 : endif 246 : enddo 247 : 248 : c The strain hardening parameters can only grow, they cannot decrease. 249 : if ((xnew(i1).le.r0).and.(xnew(i2).le.r0)) then 250 : lnoconv = .true. 251 : lerror = .true. 252 : return 253 : endif 254 : 255 : dxcor1 = xnew(1) - x(1) 256 : dxcor2 = xnew(2) - x(2) 257 : 258 : c ************************** 259 : c ! CALCULATION OF DSIGMA ! 260 : c ************************** 261 : c 261 : call sigf3D_implicit(esigmain,dxcor1,dxcor2,esigmaf,H66, 262 : . P2,pi2,ag,i6,fc,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 263 : if (lerror) then 264 : return 265 : endif 266 : do i=1,6 267 : dsigma(i) = esigmaf(i) - esigmain(i) 268 : enddo 269 : 270 : c *********************************** 271 : c ! PLANE STRESS ITERATIONS ! 272 : c *********************************** 273 : c 274 : bb = abs(dsigma(3)+esigmaf(3)) 275 : esigmaf(3) = -dsigma(3) 276 : esigmaf(1) = esigmain(1)-dsigma(3)*poison/(1.d0-poison) 277 : esigmaf(2) = esigmain(2)-dsigma(3)*poison/(1.d0-poison) 278 : if (bb.gt.(10.d0)) then 279 : if (iderin.eq.0) then 280 : iter = iter-1 281 : endif 282 : iderin = iderin+1 283 : if (iderin.gt.50) then 284 : write(6,*) ' problem in internal iterations plane stress' 285 : return 286 : endif 287 : go to 10 288 : else 289 : continue 290 : endif 291 : 292 : c ********************************* 293 : c ! UPDATE OF THE STRESSES ! 294 : c ********************************* 295 : c 296 : do i=1,6 297 : esigmaf(i) = esigmaf(i) + dsigma(i) 298 : end do 299 : esigmaf(3) = 0.d0 300 : do i=1,6 301 : esigmain(i) = esigmaf(i) 302 : end do 303 : 304 : c ********************************************* 305 : c ! UPDATE OF THE HARDENING PARAMETERS ! 306 : c ********************************************* 307 : X(i1) = XNEW(i1) 308 : X(i2) = XNEW(i2) 309 : 310 : c update of the hardening parameters 311 : rkappait= parahot3(idimpara3-37) 312 : rkappaic= parahot3(idimpara3-36) 313 : rkappat = rkappait+x(1) 314 : rkappac = rkappaic+x(2) 315 : 316 : c **************************** 317 : c ! YIELD FUNCTIONS ! 318 : c **************************** 319 : c 320 : c calculate f[0.d0,0.d0,esigmaf,rkappat,rkappac] 320 : call find_f3d(fn,esigmaf,rkappat,rkappac,0.d0,0.d0, 321 : . alpha,pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot, 322 : . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) 323 : if (lerror) then 324 : return 325 : endif 326 : if (fn(i1).gt.r10) rk(1)=1. 327 : if (fn(i2).gt.r10) rk(2)=1. 328 : 329 : rkt = rk(1) 330 : rkc = rk(2) 331 : 332 : X(i1) = XNEW(i1) 333 : X(i2) = XNEW(i2) 334 : 335 : F(i1) = FN(i1) 336 : F(i2) = FN(i2) 337 : 338 : c ******************************* 339 : c ! TEST TO END THE LOOP ! 340 : c ******************************* 341 : 342 : test = r0 343 : if (iter.ne.i1) then 344 : 345 : c The plasticity functions must be satisfied with a precision of 10 Pa, i.e. 0.000010 N/m² 346 : if ((dabs(rk(1)*fn(i1)).le.10.).and. 347 : . (dabs(rk(2)*fn(i2)).le.10.)) then 348 : c convergence has been obtained 349 : return 350 : endif 351 : c 352 : do k=1,2 353 : if (dx(k).ne.r0) then 354 : temp = abs(rk(k)*dx(k)/(1.E-10)) 355 : if (temp.gt.test) test = temp 356 : endif 357 : end do 358 : if (test.lt.tolx) then 359 : if ((dabs(rk(1)*fn(i1)).le.100.).and. 360 : . (dabs(rk(2)*fn(i2)).le.100.)) then 361 : return 362 : else 363 : lnoconv = .true. 364 : return 365 : endif 366 : endif 367 : endif 368 : c END OF THE LOOP EITHER BECAUSE THE EQUATION IS SOLVED 369 : c OR BECAUSE THE X VARIABLE DOES NOT VARY ANY MORE 370 : c OTHERWISE, GO ON WITH THE ITERATIONS 371 : c ************************************ 372 : 373 : if (iter.lt.70) then 374 : go to 11 375 : else 376 : if ((dabs(rk(1)*fn(i1)).le.10000.).and. 377 : . (dabs(rk(2)*fn(i2)).le.10000.) 378 : . .and.(lonechance)) then 379 : iter = 2 380 : lonechance = .false. 381 : go to 11 382 : else 383 : lerror = .true. 384 : lnoconv = .true. 385 : endif 386 : endif 387 : 388 : RETURN 389 : END 390 : 391 : 392 :
© Cast3M 2003 - All rights reserved.
Disclaimer