Numérotation des lignes :

algop2
C ALGOP2    SOURCE    CB215821  16/04/21    21:15:07     8920      Subroutine algoplas2(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 DPRAN3Dc     the system is solved by the method Newton c           rkt x function_t(x(1),x(2)) + (1-rkt) x(1) = 0c           rkc x function_c(x(1),x(2)) + (1-rkc) x(2) = 0 c     withc           rkt and rkc = 0 or 1c           x(1)=Dlambdatc           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(6),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),rk(2),esigmaf(6),dsiglc(6)      dimension B11(1),B12(1),B21(1),B22(1),xtest(2),dsiglt(6)       DATA itermax /90/, tolx /1.e-4/  ******     MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS)      IF(I1.GT.1) PRINT *, ' ALGOPLAS2 - ERREUR  I1 = ', I1, ' > 1 '      IF(I2.GT.2) PRINT *, ' ALGOPLAS2 - ERREUR  I2 = ', I2, ' > 2 '      IF(I6.GT.6) PRINT *, ' ALGOPLAS2 - 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))      rk(1) = rkt      rk(2) = rkc   87  continue c     initialization of fn      FN(i1) = r0      FN(i2) = r0 c     initialization      do iloc=1,6        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     !             initialization of JACOBIAN B                  !c     !             ============================                  !c     ============================================================= c     initialization of the Jacobian B by the following equationc     B(11) = - (dft/dsigma dsigma/dlambdat - dtaut)c     B(12) = - (rk(1)*rk(2)*dft/dsigma dsigma/dlambdac)c     B(21) = - (rk(2)*rk(1)*dfc/dsigma dsigmaLdlambdat)c     B(22) = - (dfc/dsigma dsigmaLdlambdac - (1-alpha)*dtauc)       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.        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 DE DX_i !c       **************************c                   -1c       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.(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       *************************cc       X_i+1 = X_i + DX_ic        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 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       !   UPDATE OF THE STRESSES   !c       ******************************         call sigf3D_implicit(esigmae,xnew(1),xnew(2),esigmaf,H66,P2,     .      pi2,ag,i6,fc,fb,parahot3,idimpara3,lerror,lcp,esigmae,U)        if (lerror) then          return        endif c       **************************c       !     YIELD FUNCTIONS    !c       **************************cc       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)         F(i1) = FN(i1)        F(i2) = FN(i2) c       *************************************c       !   UPDATE OF THE INVERT 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.          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       !    TEST END OF THE LOOP    !c       ******************************         test = r0        if (iter.ne.i1) thenc         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.)) thenc         convergence has been obtained            return          endif           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        endifc       END OF THE LOOP EITHER BECAUSE THE EQUATION IS SOLVEDc                           OR BECAUSE THE X VARIABLE DOES NOT VARY ANY MOREc       OTHERWISE, GO ON WITH THE ITERATIONSc       ************************************         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      

© Cast3M 2003 - Tous droits réservés.
Mentions légales