Numérotation des lignes :

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 DPRAN3Dc     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) = 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)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 equationc     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                   -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-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       ***************************cc         X_i+1 = X_i + DX_ic        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        endifc       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*rlocc       xnew(2) = dlambdat        rloc = xnew(i2)        xnew(i2) = rkc*rlocc       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*rlocc         xnew(1) = dlambdat          rloc = xnew(i2)          xnew(i2) = rkc*rlocc         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+1c         Now, we apply the formula of Sherman-Morrisonc                    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) thenc         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.)) thenc           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        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