Numérotation des lignes :

algop2cp
C ALGOP2CP  SOURCE    CB215821  16/04/21    21:15:07     8920      Subroutine Algoplas2cp(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 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) = 0c     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(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 parahot3(idimpara3),rk(2),xtest(2),dsiglt(6),dsiglc(6)      dimension B11(1),B12(1),B21(1),B22(1)       DATA itermax /90/, tolx /1.e-4/  ******     MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS)      IF(I1.GT.1) PRINT *, ' ALGOPLAS2CP - ERREUR  I1 = ', I1, ' > 1 '      IF(I2.GT.2) PRINT *, ' ALGOPLAS2CP - ERREUR  I2 = ', I2, ' > 2 '      IF(I6.GT.6) PRINT *, ' ALGOPLAS2CP - 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.      lonechance = .true.      G = Eo/(r2*(r1+poison))      rk(1) = rkt      rk(2) = rkc c     initialization of fn      FN(i1) = r0      FN(i2) = r0 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      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     !                  loop "Newton-Raphson                         !c     !                  ====================                         !c     =================================================================    11 iter=iter+i1   10   continue  c       *************************************c       !    CALCULATION OF THE 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.          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       !  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-18) dx(i1) = r0        if (ABS(dx(i2)) .lt. 1.E-18) dx(i2) = r0        if (((ABS(dx(i1)) .gt. 1.).and.(rk(1).eq.1.)).or.     .      ((ABS(dx(i2)) .gt. 1.).and.(rk(2).eq.1.))) then          lnoconv = .true.          lerror = .true.          return        endif        if ( (ABS(dx(i1)) .gt. 1.).and.(rk(1).eq.0.) )     .     dx(1) = 1. * dx(1)/ABS(dx(i1))        if ( (ABS(dx(i2)) .gt. 1.).and.(rk(2).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         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))        esigmaf(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 in 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 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       !      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)         X(i1) = XNEW(i1)        X(i2) = XNEW(i2)         F(i1) = FN(i1)        F(i2) = FN(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(rk(1)*fn(i1)).le.10.).and.     .        (ABS(rk(2)*fn(i2)).le.10.)) thenc         convergence has been obtained            return          endifc          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       ************************************       if (iter.lt.70) then        go to 11      else        if ((ABS(rk(1)*fn(i1)).le.10000.).and.     .      (ABS(rk(2)*fn(i2)).le.10000.)     .   .and.(lonechance)) then          iter = 2          lonechance = .false.          go to 11        else          lerror = .true.          lnoconv = .true.        endif      endif       RETURN      END      

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