C ALGOP1CP SOURCE FANDEUR 14/03/25 21:15:04 7993 Subroutine algoplas1cp(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 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) 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 deps(6),esigmain(6),dsigma(6),esigmaf(6) dimension parahot3(idimpara3) DATA itermax /90/, tolx /1.e-4/ ***** * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS) IF(I1.GT.1) PRINT *, ' ALGOPLAS1CP - ERREUR I1 = ', I1, ' > 1 ' IF(I2.GT.2) PRINT *, ' ALGOPLAS1CP - ERREUR I2 = ', I2, ' > 2 ' IF(I6.GT.6) PRINT *, ' ALGOPLAS1CP - 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. G = Eo/(r2*(r1+poison)) lfirstiter = .true. 87 continue c initialization of fn FN(i1) = r0 FN(i2) = r0 c initialization of the activated yield surfaces lftact = .false. lfcact = .false. 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 dxprevt = 0.0d0 dxprevc = 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.) rkt = 0.d0 f(i1) = rkt*f(i1) f(i2) = rkc*f(i2) if (f(i1).gt.r0) lftact = .true. if (f(i2).gt.r0) lfcact = .true. fn(i1) = f(i1) fn(i2) = f(i2) c ================================================================= c ! ! c ! loop "Newton-Raphson ! c ! ==================== ! c ================================================================= 11 iter=iter+i1 10 continue if (iter.lt.2) then 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(esigmaf,dft,i3,i6,parahot3,idimpara3, . rkappait,lcp) c ******** do iloc=i1,i6 DFT(iloc) = rkt * DFT(iloc) enddo call dsigDrucker(fc,fb,esigmaf,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(rkappat) dtauc = dfhardcc(rkappac,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. 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 else c ********************************************************** c ! UPDATE THE INVERT JACOBIAN BY PAR SHERMAN-MORRISON ! c ********************************************************** c ! YIELD FUNCTIONS ! 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 c Reactivation of yield surface if (fn(i1).gt.r10) rkt=1. if (fn(i2).gt.r10) rkc=1. fn(i1) = rkt*fn(i1) fn(i2) = rkc*fn(i2) rloc = xnew(i1) xnew(i1) = rkt*rloc c xnew(1) = dlambdat rloc = xnew(i2) xnew(i2) = rkc*rloc c xnew(2) = dlambdac 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 X(i1) = XNEW(i1) X(i2) = XNEW(i2) F(i1) = FN(i1) F(i2) = FN(i2) 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-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)) if (iter.gt.2) then r_z = 1.50d0 * dxprevt if (abs(dx(1)).gt.abs(r_z)) then rloc = dx(1) dx(1) = r_z * (rloc/abs(rloc)) endif r_z = 1.50d0 * dxprevc if (abs(dx(2)).gt.abs(r_z)) then rloc = dx(2) dx(2) = r_z * (rloc/abs(rloc)) endif endif dxprevt = abs(dx(1)) dxprevc = abs(dx(2)) 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) 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 if (xnew(i1).eq.r0) then rkt = r0 endif if (xnew(i2).eq.r0) then rkc = r0 endif c new calculation of dx if we have reactivated a yield surface if ((xnew(i1).lt.r0).and.(rkc.eq.1.)) then rkt = r0 if (lfirstiter) then x(1) = 0.d0 x(2) = 0.d0 lfirstiter = .false. go to 87 else if (iter.lt.70) then go to 11 else lerror = .true. lnoconv = .true. return endif endif if ((xnew(i2).lt.r0).and.(rkt.eq.1.)) then rkc = r0 if (lfirstiter) then x(1) = 0.d0 x(2) = 0.d0 lfirstiter = .false. go to 87 else if (iter.lt.70) then go to 11 else lerror = .true. lnoconv = .true. return endif endif if (xnew(i1).gt.r0) rkt=1. if (xnew(i2).gt.r0) rkc=1. rloc = xnew(i1) xnew(i1) = rkt*rloc c xnew(1) = dlambdat rloc = xnew(i2) xnew(i2) = rkc*rloc c xnew(2) = dlambdac c dxcor(rected) for calculation of dsigma c dx = xnew-x because if xnew<0 => xnew=0 and the dx must be corrected 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) = esigmain(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 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 VARIABLES ! 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 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) rkt=1. if (fn(i2).gt.r10) rkc=1. fn(i1) = rkt*fn(i1) fn(i2) = rkc*fn(i2) c ********************************** c ! TEST TO END THE LOOP ! c ********************************** if ((abs(dx(i1)).lt.1.E-10).and.(abs(dx(i2)).lt.1.E-10).and. & (abs(fn(i1)).le.10.).and.(abs(fn(i2)).le.10.)) then return endif 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 (dx(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 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 ************************************ if (iter.lt.70) then go to 11 else lerror = .true. lnoconv = .true. endif RETURN END