C ALGOP1 SOURCE CB215821 16/04/21 21:15:06 8920 Subroutine algoplas1(rkt,rkc,esigmae,H66,P2,pi2,alpha, . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot, . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6, . parahot3,idimpara3,deps,lnoconv,lcp,U) c SUBROUTINE added by THG c THIS SUBROUTINE IS CALLED BY DPRAN3D c the system is solved by the method of Broyden (Newton-Raphson secante) c rkt x function_t(x(1),x(2)) + (1-rkt) x(1) = 0 c rkc x function_c(x(1),x(2)) + (1-rkc) x(2) = 0 c with c rkt and rkc = 0 or 1 c x(1)=Dlambdat c x(2)=Dlambdac IMPLICIT REAL*8 (A-B,D-H,O-Z) implicit integer (I-K,M,N) implicit logical (L) implicit character*10 (C) **** dimension B(i2,i2) dimension B(2,2) c jacobian matrix for the system of equations. **** dimension Binv(i2,i2) dimension Binv(2,2) c invert of the jacobian matrix for the system of equations. **** dimension df(i2) dimension df(2) **** dimension dfc(i6) dimension dfc(6) c derivative of the compression plasticity function **** dimension dft(i6) dimension dft(6) c derivative of the tension plasticity function **** dimension dx(i2) dimension dx(2) c increment of the solution **** dimension f(i2) dimension f(2) **** dimension fn(i2) dimension fn(2) **** dimension H66(i6,i6) dimension H66(6,6) c matrix of elasticity (sometimes denoted D) **** dimension x(i2) dimension x(2) c solution to be found (and initial guess at the entry of the subr.) **** dimension xnew(i2) dimension xnew(2) dimension esigmae(i6) c trial stress vector, = elastic predictor **** dimension trav1(i1),trav2(i2),trav2b(i2),trav2c(i2), **** . trav6(i6),trav22(i2,i2) dimension trav1(1),trav2(2),trav2b(2),trav2c(2), . trav6(6),trav22(2,2) **** dimension P2(i6,i6),U(i6,i6),pi2(i6),v2loc(i2,i2) dimension P2(6,6),U(6,6),pi2(6),v2loc(2,2) dimension deps(i6) dimension parahot3(idimpara3) DATA itermax /90/, tolx /1.e-4/ ***** * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS) IF(I1.GT.1) PRINT *, ' ALGOPLAS1 - ERREUR I1 = ', I1, ' > 1 ' IF(I2.GT.2) PRINT *, ' ALGOPLAS1 - ERREUR I2 = ', I2, ' > 2 ' IF(I6.GT.6) PRINT *, ' ALGOPLAS1 - ERREUR I6 = ', I6, ' > 6 ' ***** i7 = 7 i8 = 8 i9 = 9 i10=10 i11=11 i12=12 i13=13 i14=14 i15=15 i16=16 i17=17 i18=18 i19=19 i20=20 i21=21 i22=22 i23=23 i24=24 i25=25 i26=26 i27=27 i28=28 i29=29 i30=30 i31=31 i32=32 i33=33 i34=34 i35=35 i36=36 i37=37 i38=38 i39=39 r1 = 1. r2 = 2. r0 = 0. r10 = 10. i0 = 0 i1 = 1 lerror = .false. lproblem = .false. G = Eo/(r2*(r1+poison)) 87 continue c initialization of fn FN(i1) = r0 FN(i2) = r0 c initialization of the activated yield surfaces lftact = .false. lfcact = .false. c calculate a first f[x(1),x(2),esigmae,rkappait,rkappaic] call find_f3d(f,esigmae,rkappait,rkappaic,x(i1),x(i2),alpha, . pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot, . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) if (lerror) then return endif f(i1) = rkt*f(i1) f(i2) = rkc*f(i2) if (f(i1).gt.r0) lftact = .true. if (f(i2).gt.r0) lfcact = .true. if ((f(i1).lt.10.).and.(f(i2).lt.10.)) then return endif fn(i1) = f(i1) fn(i2) = f(i2) c ================================================================= c ! ! c ! initialization of Jacobian B ! c ! ============================ ! c ================================================================= c initialization of Jacobian B in x = (0;0) by the following equation c B(11) = - (dft/dsigma H66 dft/dsigma + dtaut) c B(12) = - (dft/dsigma H66 dfc/dsigma) c B(21) = - (dfc/dsigma H66 dft/dsigma) c B(22) = - (dfc/dsigma H66 dfc/dsigma + dtauc) call dsigRank(esigmae,dft,i3,i6,parahot3,idimpara3,rkappait,lcp) c ******** do iloc=i1,i6 DFT(iloc) = rkt * DFT(iloc) enddo call dsigDrucker(fc,fb,esigmae,dfc,alpha,P2,pi2,i1,i6,lcp) c *********** do iloc=i1,i6 DFC(iloc) = rkc * DFC(iloc) enddo call mulATB(dft,H66,trav6,i6,i1,i6,i6) c ****** call mulAB(trav6,dft,trav1,i1,i6,i6,i1) c ***** df11 = trav1(i1) call mulAB(trav6,dfc,trav1,i1,i6,i6,i1) c ***** df12 = trav1(i1) df21 = df12 call mulATB(dfc,H66,trav6,i6,i1,i6,i6) c ****** call mulAB(trav6,dfc,trav1,i1,i6,i6,i1) c ***** df22 = trav1(i1) c derivatives of the hardening functions dtaut = dfhardct(rkappait) dtauc = dfhardcc(rkappaic,fc,fc0,rkappa1,ac,bc) B(i1,i1) = -(df11+rkt*dtaut) +(r1-rkt) B(i1,i2) = -df12 B(i2,i1) = -df21 B(i2,i2) = -(df22+rkc*dtauc) +(r1-rkc) c ================================================================= c ! ! c ! invert the initial JACOBIAN ! c ! =========================== ! c ================================================================= dtmB = B(i1,i1)*B(i2,i2)-B(i1,i2)*B(i2,i1) if (dtmB.eq.r0) then lerror = .true. lnoconv = .true. c print *,'ERROR BROYDEN dtmB=0' return else Binv(i1,i1) = B(i2,i2)/dtmB Binv(i1,i2) =-B(i1,i2)/dtmB Binv(i2,i1) =-B(i2,i1)/dtmB Binv(i2,i2) = B(i1,i1)/dtmB endif c ================================================================= c ! ! c ! loop "Newton-Raphson ! c ! ==================== ! c ================================================================= do iter=i1,itermax c ************************* c ! CALCULATION OF DX_i ! c ************************* c -1 c DX_i = - J_i F(X_i) c TRAV22(i1,i1) = -BINV(i1,i1) TRAV22(i1,i2) = -BINV(i1,i2) TRAV22(i2,i1) = -BINV(i2,i1) TRAV22(i2,i2) = -BINV(i2,i2) call mulab(trav22,F,DX,i2,i2,i2,i1) if (ABS(dx(i1)) .lt. 1.E-16) dx(i1) = r0 if (ABS(dx(i2)) .lt. 1.E-16) dx(i2) = r0 if (((ABS(dx(i1)) .gt. 1.).and.(rkt.eq.1.)).or. . ((ABS(dx(i2)) .gt. 1.).and.(rkc.eq.1.))) then lnoconv = .true. lerror = .true. return endif if ( (ABS(dx(i1)) .gt. 1.).and.(rkt.eq.0.) ) . dx(1) = 1. * dx(1)/ABS(dx(i1)) if ( (ABS(dx(i2)) .gt. 1.).and.(rkc.eq.0.) ) . dx(2) = 1. * dx(2)/ABS(dx(i2)) c *************************** c ! CALCULATION OF X_i+1 ! c *************************** c c X_i+1 = X_i + DX_i c XNEW(i1) = X(i1) + DX(i1) XNEW(i2) = X(i2) + DX(i2) * PRINT *,'ALGOPLAS1 X(I1),DX(I1),XNEW(I1)',X(I1),DX(I1),XNEW(I1) * PRINT *,'ALGOPLAS1 X(I2),DX(I2),XNEW(I2)',X(I2),DX(I2),XNEW(I2) c The strain hardening parameters can only grow, they cannot decrease. 114 if ((xnew(i1).lt.r0).and.(xnew(i2).lt.r0)) then lnoconv = .true. return endif c From Simo et Hugues, Computational Inelasticity, p.212 procedure 2 if (xnew(i1).lt.r0) then rkt = r0 go to 87 endif if (xnew(i2).lt.r0) then rkc = r0 go to 87 endif if ((xnew(i1).ge.r0).and.(lftact)) rkt=1. if ((xnew(i2).ge.r0).and.(lfcact)) rkc=1. rloc = xnew(i1) xnew(i1) = rkt*rloc c xnew(2) = dlambdat rloc = xnew(i2) xnew(i2) = rkc*rloc c xnew(2) = dlambdac c ********************************************************** c ! UPDATE OF THE INVERT JACOBIAN BY SHERMAN-MORRISON ! c ********************************************************** if ((dx(i1).ne.r0).or.(dx(i2).ne.r0)) then rloc = xnew(i1) xnew(i1) = rkt*rloc c xnew(1) = dlambdat rloc = xnew(i2) xnew(i2) = rkc*rloc c xnew(2) = dlambdac c calculate f[xnew(1),xnew(2),esigmae,rkappait,rkappaic] call find_f3d(fn,esigmae,rkappait,rkappaic,xnew(i1),xnew(i2), . alpha,pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot, . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) if (lerror) then return endif c test activated yield surfaces if ((fn(i1).gt.r10).and.(lftact)) rkt = r1 if ((fn(i2).gt.r10).and.(lfcact)) rkc = r1 fn(i1) = rkt*fn(i1) fn(i2) = rkc*fn(i2) DF(i1) = FN(i1) - F(i1) DF(i2) = FN(i2) - F(i2) c calculation of the new Jacobian Binv_i+1 c Now, we apply the formula of Sherman-Morrison c Binv.df call mulab(Binv,DF,TRAV2,i2,i2,i2,i1) c dx - Binv.df TRAV2b(i1) = DX(i1)-TRAV2(i1) TRAV2b(i2) = DX(i2)-TRAV2(i2) c dx.Binv call mulATB(dx,Binv,trav2c,i2,i1,i2,i2) c ****** c (dx-Binv.df).(dx.Binv) call mulAB(trav2b,trav2c,trav22,i2,i1,i1,i2) c ***** c dx.Binv.df call mulATB(dx,trav2,trav1,i2,i1,i2,i1) c ****** rloc = r1/trav1(i1) do jloc=i1,i2 do iloc=i1,i2 TRAV22(iloc,jloc) = rloc*TRAV22(iloc,jloc) v2loc(iloc,jloc) = BINV(iloc,jloc) BINV(iloc,jloc) = v2loc(iloc,jloc) + TRAV22(iloc,jloc) end do end do endif X(i1) = XNEW(i1) X(i2) = XNEW(i2) c ****************************** c ! TEST TO END THE LOOP ! c ****************************** test = r0 if (iter.ne.i1) then c The plasticity functions must be satisfied with a precision of 10 Pa, i.e. 0.000010 N/m² if ((ABS(fn(i1)).le.10.).and.(ABS(fn(i2)).le.10.)) then c convergence has been obtained return endif do k=1,2 if (x(k).ne.r0) then temp = abs(dx(k)/(1.E-10)) if (temp.gt.test) test = temp endif end do if (test.lt.tolx) then lnoconv = .true. return endif endif c END OF THE LOOP EITHER BECAUSE THE EQUATION IS SOLVED c OR BECAUSE THE X VARIABLE DOES NOT VARY ANY MORE c OTHERWISE, GO ON WITH THE ITERATIONS c ************************************ F(i1) = FN(i1) F(i2) = FN(i2) end do write(2,1) test,f(i1),f(i2),ntot 1 format( ' itermax exceeded in Broyden, test:',e7.2, . ', ft:',E8.2,' N/m², fc:',e8.2,' N/m²'/ . ' ntot = ',i6) lerror = .true. lnoconv = .true. RETURN END