C ALGOP2CP SOURCE CB215821 16/04/21 21:15:07 8920 Subroutine Algoplas2cp(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,esigmaf,rkappat,rkappac, . lnoconv,lcp,U) c SUBROUTINE added by THG c THIS SUBROUTINE IS CALLED BY DPRAN3D c the system is solved by the method Newton 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) dimension esigmae(6) 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),esigmain(6),dsigma(6),esigmaf(6) dimension parahot3(idimpara3),rk(2),xtest(2),dsiglt(6),dsiglc(6) dimension B11(1),B12(1),B21(1),B22(1) DATA itermax /90/, tolx /1.e-4/ ***** * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS) IF(I1.GT.1) PRINT *, ' ALGOPLAS2CP - ERREUR I1 = ', I1, ' > 1 ' IF(I2.GT.2) PRINT *, ' ALGOPLAS2CP - ERREUR I2 = ', I2, ' > 2 ' IF(I6.GT.6) PRINT *, ' ALGOPLAS2CP - 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. r0 = 0. r2 = 2. r10 = 10. lerror = .false. lproblem = .false. lonechance = .true. G = Eo/(r2*(r1+poison)) rk(1) = rkt rk(2) = rkc c initialization of fn FN(i1) = r0 FN(i2) = r0 c initialization counter loop plane stress iderin = 0 iter = 0 c initialization do iloc=1,6 esigmain(iloc) = esigmae(iloc) esigmaf(iloc) = esigmae(iloc) enddo rkappat = rkappait rkappac = rkappaic test = 0.d0 xnew(1) = 0.0d0 xnew(2) = 0.0d0 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 if ((f(i1).lt.10.).and.(f(i2).lt.10.)) then return endif if (f(1).lt.10.) rk(1) = 0.d0 fn(i1) = f(i1) fn(i2) = f(i2) c ================================================================= c ! ! c ! loop "Newton-Raphson ! c ! ==================== ! c ================================================================= 11 iter=iter+i1 10 continue c ************************************* c ! CALCULATION OF THE JACOBIAN ! c ************************************* call dsigRank(esigmaf,dft,i3,i6,parahot3,idimpara3, . rkappait,lcp) call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2,pi2,i1,i6,lcp) call dsiglambda(xnew(1),xnew(2),P2,pi2,ag,fc,fb,lcp,rkappait, . parahot3,idimpara3,U,H66,esigmaf,dsiglt,dsiglc) call mulATB(dft,dsiglt,B11,6,1,6,1) call mulATB(dft,dsiglc,B12,6,1,6,1) call mulATB(dfc,dsiglt,B21,6,1,6,1) call mulATB(dfc,dsiglc,B22,6,1,6,1) dtaut = dfhardct(rkappat) dtauc = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc) B(1,1) = B11(1) - dtaut B(1,2) = rk(1) * rk(2) * B12(1) B(2,1) = rk(2) * rk(1) * B21(1) B(2,2) = B22(1) - (1.-alpha) * dtauc dtmB = B(i1,i1)*B(i2,i2)-B(i1,i2)*B(i2,i1) if (dtmB.eq.r0) then lerror = .true. 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 ! 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-18) dx(i1) = r0 if (ABS(dx(i2)) .lt. 1.E-18) dx(i2) = r0 if (((ABS(dx(i1)) .gt. 1.).and.(rk(1).eq.1.)).or. . ((ABS(dx(i2)) .gt. 1.).and.(rk(2).eq.1.))) then lnoconv = .true. lerror = .true. return endif if ( (ABS(dx(i1)) .gt. 1.).and.(rk(1).eq.0.) ) . dx(1) = 1. * dx(1)/ABS(dx(i1)) if ( (ABS(dx(i2)) .gt. 1.).and.(rk(2).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 do iloc=1,2 Xtest(iloc) = X(iloc) + DX(iloc) if ((Xtest(iloc).le.(0.)).and.(rk(iloc).eq.1)) then do jloc=1,2 if ((rk(jloc).eq.1).and.(jloc.ne.iloc)) . DX(jloc) = - f(jloc)*BINV(jloc,jloc) enddo endif enddo XNEW(i1) = X(i1) + DX(i1) XNEW(i2) = X(i2) + DX(i2) XNEW(i1) = max(0.D0,XNEW(i1)) XNEW(i2) = max(0.D0,XNEW(i2)) do iloc=1,2 if (XNEW(iloc).eq.(0.)) then rk(iloc) = 0. else rk(iloc) = 1. endif enddo c The strain hardening parameters can only grow, they cannot decrease. if ((xnew(i1).le.r0).and.(xnew(i2).le.r0)) then lnoconv = .true. lerror = .true. return endif dxcor1 = xnew(1) - x(1) dxcor2 = xnew(2) - x(2) c ************************** c ! CALCULATION OF DSIGMA ! c ************************** c call sigf3D_implicit(esigmain,dxcor1,dxcor2,esigmaf,H66, . P2,pi2,ag,i6,fc,fb,parahot3,idimpara3,lerror,lcp,esigmae,U) if (lerror) then return endif do i=1,6 dsigma(i) = esigmaf(i) - esigmain(i) enddo c *********************************** c ! PLANE STRESS ITERATIONS ! c *********************************** c bb = abs(dsigma(3)+esigmaf(3)) esigmaf(3) = -dsigma(3) esigmaf(1) = esigmain(1)-dsigma(3)*poison/(1.d0-poison) esigmaf(2) = esigmain(2)-dsigma(3)*poison/(1.d0-poison) if (bb.gt.(10.d0)) then if (iderin.eq.0) then iter = iter-1 endif iderin = iderin+1 if (iderin.gt.50) then write(6,*) ' problem in internal iterations plane stress' return endif go to 10 else continue endif c ********************************* c ! UPDATE OF THE STRESSES ! c ********************************* c do i=1,6 esigmaf(i) = esigmaf(i) + dsigma(i) end do esigmaf(3) = 0.d0 do i=1,6 esigmain(i) = esigmaf(i) end do c ********************************************* c ! UPDATE OF THE HARDENING PARAMETERS ! c ********************************************* X(i1) = XNEW(i1) X(i2) = XNEW(i2) c update of the hardening parameters rkappait= parahot3(idimpara3-37) rkappaic= parahot3(idimpara3-36) rkappat = rkappait+x(1) rkappac = rkappaic+x(2) c **************************** c ! YIELD FUNCTIONS ! c **************************** c c calculate f[0.d0,0.d0,esigmaf,rkappat,rkappac] call find_f3d(fn,esigmaf,rkappat,rkappac,0.d0,0.d0, . 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 if (fn(i1).gt.r10) rk(1)=1. if (fn(i2).gt.r10) rk(2)=1. rkt = rk(1) rkc = rk(2) X(i1) = XNEW(i1) X(i2) = XNEW(i2) F(i1) = FN(i1) F(i2) = FN(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(rk(1)*fn(i1)).le.10.).and. . (ABS(rk(2)*fn(i2)).le.10.)) then c convergence has been obtained return endif c do k=1,2 if (dx(k).ne.r0) then temp = abs(rk(k)*dx(k)/(1.E-10)) if (temp.gt.test) test = temp endif end do if (test.lt.tolx) then if ((ABS(rk(1)*fn(i1)).le.100.).and. . (ABS(rk(2)*fn(i2)).le.100.)) then return else lnoconv = .true. return endif 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 ************************************ if (iter.lt.70) then go to 11 else if ((ABS(rk(1)*fn(i1)).le.10000.).and. . (ABS(rk(2)*fn(i2)).le.10000.) . .and.(lonechance)) then iter = 2 lonechance = .false. go to 11 else lerror = .true. lnoconv = .true. endif endif RETURN END