Télécharger algop2cp.eso

Retour à la liste

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 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(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                   -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.(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       *************************
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
 
        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       ****************************
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)
 
        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.)) then
c         convergence has been obtained
            return
          endif
c
          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       ************************************
 
      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