Télécharger algop2.eso

Retour à la liste

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 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(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 equation
c     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                   -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.(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
        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       **************************
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)
 
        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) 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
 
          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       ************************************
 
        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