Télécharger algop1.eso

Retour à la liste

algop1
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
 
 
 
 
 
 

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