Télécharger sjonl3.eso

Retour à la liste

sjonl3
C SJONL3    SOURCE    CHAT      06/06/01    21:20:08     5450
      SUBROUTINE SJONL3(sigi, depst, vari, xmat, tras, jtras ,
     >                  trat, jtrat, trac, jtrac, sigf, varf,
     >                  depsp, ireturn)
c --- SJONL3 ---
c computing plastic flow in joint element
c under 3D plane (strain) conditions
c
c beat release 24-7-00 ylp.
c sigi    IN  initial stress
c depst   IN  total strain increment
c vari    IN  initial internal variable
c xmat    IN  material properties
c tras    IN  shear plastic softening behaviour
c jtras   IN  <tras> number of points
c trat    IN  tensile plastic softening behaviour
c jtrat   IN  <trat> number of points
c trac    IN  compression plastic behaviour
c jtrac   IN  <trac> number of points
c sigf    OUT final stress
c varf    OUT final internal variable
c depsp   OUT plastic strain increment
c ireturn OUT routine control integer
 
         IMPLICIT REAL*8(A-H,K-Z)
 
         IMPLICIT INTEGER(I-J)
         dimension sigi(*), depst(*), vari(*),
     >           sigf(*), depsp(*), varf(*)
         dimension ddepst(3), ddepsp(3)
         dimension xmat(*)
       dimension tras(2,jtras), trat(2,jtrat), trac(2,jtrac)
c 1 --> strain
c 2 --> stress
c_______________________________________________________________________
c
*       save iloop
*       data iloop/0/
c        print *,'debug -- sjonl3d -- check everything'
c
c       do i=1,3
c        print *,'depst ',depst(i)
c        print *,'vari ', vari(i)
c        print *,'sigi ', sigi(i)
c        print *,'sigf ', sigf(i)
c        print *,'varf ',varf(i)
c        print *,'depst ',depst(i)
c       enddo
c       do i=1,6
c        print *,'xmat ', xmat(i)
c       enddo
c        print *,'tras ', tras ,' jtras ', jtras
c      print *,'trat ', trat ,' jtrat ', jtrat
c        print *,'trac ', trac ,' jtrac ', jtrac
c
c        print *,'ireturn ', ireturn
*       iloop = iloop + 1
c look at xmat(6) and deduce elastic/ plastic...
c ...coupled/uncoupled variables.
c all variables Jxxx are integer used as logical with the
c following rule (0 is true, 1 is false)
c JCP : use plastic law in compression
c JSP : use plastic law in shear
c JTP : use plastic law in tension
c JST : couple shear and tension
c JSC : couple shear and compression
       call sjoncpl(xmat(6),JCP,JSP,JTP,JST,JSC)
*       print *,'DEBUG --> entree dans sjonl2'
*       print *,'DepsN DepsS ' ,depst(2), depst(1)
* set some usefull variables
         ZERO = 0.0d0
         ONE  = 1.0d0
c set the initial location of the edge of the cone
       ETi = xmat(5)
c index definition
       Ipress  = 2
         Ishear  = 1
         Itotal  = 3
 
         Ishea1  = 1
         Ishea2  = 2
         Inorm   = 3
         Icompr  = 4
 
         Istrain = 1
         Istress = 2
c number of increments
         INCR = 5
c______________________________________________________________________
c
c elastic initialisation of the joint stiffness --> K
c K[z]
c with z = t tensile ; s shear ; c compression
       Kt0 = xmat(Ipress)
         Ks0 = xmat(Ishear)
         Kc0 = Kt0
c       print *,'Deps - tension ',depst(Ipress)
c       print *,'Deps - shear   ',depst(Ishear)
c split the displacement increment in few smaller increments
       ddepst(Inorm)  = depst(Inorm)/INCR
         ddepst(Ishea1)  = depst(Ishea1)/INCR
         ddepst(Ishea2)  = depst(Ishea2)/INCR
c STARTING MAIN LOOP
c    [0
      do iflag =1 , INCR
c______________________________________________________________________
c
c initial residual strength
c RTS --> tensile
c RSS --> shear
c
c internal variables initialisation --> Alpha
c A[xyz]
c with x = i initial ; f final ; k intermediate computational value
c      y = p plastic ; t total
c      z = t tensile ; s shear
       Aipt = vari(Ipress)
c          \--> initial plastic tensile strain
         Aips = vari(Ishear)
c          \--> initial plastic shear strain
         Aipc = vari(Icompr)
c          \--> initial plastic compression strain
c______________________________________________________________________
c
c initial residual strength
c RTS --> tensile
c RSS --> shear
c RCS --> compression
c compute the  tensile stress RTS corresponding to the initial
c plastic tensile strain with <trat> evolution
       call yofxcu(Aipt ,trat, jtrat, RTS, dRTS, ireturn)
c           \--> Aipt initial plastic tensile strain
c                RTS  OUT
c                dRTS = dTkt/dAkps useless value
c     [1
         if (ireturn .ne. 0) return
c                             \--> ERROR
c     1]
c compute the initial cohesion residual
       call yofxcu(Aips, tras, jtras, RSS, dRSS, ireturn)
c     [1
         if (ireturn .ne. 0) return
c                             \--> ERROR
c     1]
c comput the initial compression residual
       call yofxcu(Aipc, trac, jtrac, RCS, dRCS, ireturn)
c     [1
         if (ireturn .ne. 0) return
c                             \--> ERROR
c     1]
c______________________________________________________________________
c
 
       Aitt = vari(Itotal)
c        \--> initial total tensile strain
c compute final total tensile strain...
       Aftt = Aitt + ddepst(Inorm)
c ... and the initial tensile strength...
       TT0 = trat(Istress,1)
c ... and the elastic compressive strength
       TC0 = trac(Istress,1)
c       print *,'TC0 ',TC0
c check one or two things if entering sjonl2 for first time
c     [1
* mis en commentaire par TC
*       if ( iloop .eq. 1) then
c       [2
*           if ( TT0 .eq. ZERO ) then
*             print *,'cut-off tensile strength is zero'
*             print *,'denied option'
*             return
*         endif
c       2]
c       [2
*         if ( TT0 .gt. ETi ) then
*             print *,'WARNING cut-off ouside the friction cone...'
*             print *,'... tensile strength at the edge of the cone'
*           print *,'set to cut-off tensile strength'
*             return
*           endif
c       2]
*       endif
*  fin de mise en commentaire par TC
c     1]
c initialization of the damage variables
       Dtensi = 1.0d0
         Dshear = 1.0d0
         Dcompr = 1.0d0
c WARNING the damage variables are understand here as :
c D = 1 means no damage
c D = 0 means full damage...
c which is exactly the opposite of the usual definition !!!
c compute the predicted elastic stress state
c using the secant stiffness
c
c first, look at the tension modulus...
c maximal (initial) tensile stress initialisation
c N.B.: it is thus assumed that the behaviour in tension is elastic
c plastic, and that the peak value features the end of the elastic
c regime.
c check whether tension or compression region is reached
c     [1
       if ( Aftt .gt. ZERO ) then
c                             \--> tensile region reached...
c       [2
         if ( JTP .ne. 0 ) then
c               \--> use elastic law
             Kt = Kt0
           else
c          \--> use damaged elastic predictor
c         [3
             if ( RTS .eq. ZERO ) then
c                              \--> zero residual tensile strength
               Kt = ZERO
c                  \--> tensile modulus set to zero
             else
c           \--> still got residual tensile strength
              Kt = RTS / ( Aipt + TT0 / Kt0 )
             endif
c         3]
c update tensile damage variable to enable coupling effect
c         [3
           if ( JST .eq. 0) Dtensi = RTS / TT0
c         3]
         endif
c       2]
c           print *,'Dtensi ' ,Dtensi
         JCT = 0
c JCT is an integer value used as a logical to control the location
c of the stress state with respect to tension or compression:
c tension     --> JCT = 0
c compression --> JCT = 1
       else
c       \--> compression region reached...
c       [2
         if ( JCP .ne. 0 ) then
c                \--> use elastic law
             Kc = Kc0
         else
c         \--> use damaged elastic predictor
c         [3
           if ( RCS .eq. ZERO ) then
               Kc = ZERO
             else
             Kc = RCS / (Aipc + TC0 / Kc0 )
             endif
c         3]
c update compression damage variable
c          Dcompr = Kc / Kc0 --> too fast
c ... try something else...
         endif
c       2]
         JCT = 1
         endif
c     ]1
c ... then, look at the shear modulus
       TS0 = tras(Istress,1)
c     [1
       if (JSP .ne. 0 ) then
           Ks = Ks0
         else
c       [2
         if ( RSS .eq. ZERO ) then
c                           \--> zero residual shear strength
c                                or opened bonds
             Ks = ZERO
c               \--> shear modulus set to zero
c Such case shall only appear if the bonds are opened.
c otherwise, residual cohesion is generally observed
           else
c          \--> still got residual shear strength
             Ks = RSS / (Aips + TS0 / Ks0)
           endif
c update the shear damage variable
         if (JST .eq. 0) Dshear = RSS / TS0
c(Ks / Ks0)
c       print *,'Dshear ' , Dshear
c       2]
       endif
c     1]
 
c     [1
       if (JCT .eq. 0) then
c check if previous stress state was in compression
c        [2
          if ( (Aitt * Aftt) .lt. ZERO ) then
               Tft = Kt * Aftt * Dtensi
            else
             Tft = sigi(Inorm) + ( ddepst(Inorm) * Kt * Dshear)
            endif
c        2]
          Tfs1 = sigi(Ishea1) + ( ddepst(Ishea1) * Ks * Dtensi)
          Tfs2 = sigi(Ishea2) + ( ddepst(Ishea2) * Ks * Dtensi)
            Tfs0 = sqrt((Tfs1**2)+(Tfs2**2))
       else
c check if previous stress state was in tension
c        [2
          if ( (Aitt * Aftt) .lt. ZERO ) then
               Tfc = Kc * Aftt * Dcompr
            else
c              Tfc = sigi(Ipress) + ( ddepst(Ipress) * Kc * Dshear)
            Tfc = sigi(Inorm) + ( ddepst(Inorm) * Kc )
          endif
c        2]
          Tfs1 = sigi(Ishea1) + ( ddepst(Ishea1) * Ks * Dcompr)
          Tfs2 = sigi(Ishea2) + ( ddepst(Ishea2) * Ks * Dcompr)
            Tfs0 = sqrt((Tfs1**2)+(Tfs2**2))
         endif
c     1]
 
c       print *,'stress state elastic prediction'
c       print *,'tension    compression    shear'
c       print *,Tft, Tfc, Tfs1, Tfs2
c______________________________________________________________________
c
c before checking if any kind of yielding occured, update the
c final internal variables in the case of simple elastic regime.
       Afpt = Aipt
         Afps = Aips
         Afpc = Aipc
c______________________________________________________________________
c
c now, check yielding
c N.B.: note that each mode might be activated independently
c without any restriction!!!!
c first, look at the tension/compression mode...
c     [1
       if ( JCT .eq. 0 ) then
c       [2
         if ( JTP .eq. 0 ) then
c         [3
           if ( abs(Tft) .gt. RTS*Dshear ) then
c                                   \--> YIELDING IN TENSION
c      print *,'yielding in tension...'
c            Afpt = Aftt - TT0 / Kt0
             Afpt = Aipt + (Tft-RTS*Dshear)/Kt
c       print *,'Afpt ',Afpt
c              \--> updated plastic tensile strain
*      print *,'updated plastic tensile strain ', Afpt
c compute the updated residual tensile strength...
               call yofxcu(Afpt, trat, jtrat, RTS ,dRTS, ireturn)
c                   \--> Afpt updated plastic tensile strain
c                        RTS  updated residual tensile strength
c                       dTkt dummy
             Tft = RTS * Dshear
c               \--> compute the final tensile stress
c      print *,'residual tensile strength ',Tft
           endif
c         3]
         endif
c       2]
       else
c       [2
         if ( JCP .eq. 0 ) then
c          if ( abs(Tfc) .gt. abs(RCS)*Dshear) then
c         [3
           if ( abs(Tfc) .gt. abs(RCS) ) then
c                                        \--> YIELDING IN COMPRESSION
c            print *,'yielding in compression'
c              Afpc = abs(Aftt) - TC0 / Kc0
c            Afpc = Aipc + ((Tfc*-1.0d0)-RCS*Dshear)/Kc
             Afpc = Aipc + ((Tfc*(-1.0d0))-RCS)/Kc
c              \--> updated plastic compression strain
c             print *,'updated plastic compression strain ', Afpc
c compute the updated residual compression strength
             call yofxcu(Afpc, trac, jtrac, RCS , dRCS, ireturn)
c then, go back to the compression domain (negative stress)
               Tfc = RCS * (-1.0d0)
c              print *,'residual compressive strength ',RCS
           endif
c         3]
           endif
c       2]
         endif
c     1]
c ... finally, look at the shear mode...
c compute the tensile value at the edge using an homothetic transform.
c the cut-off value plays the role of a drag-stress
c       ETf = RTS + ETi - TT0
c ... definitely not relevant!!!!
c ... try something else... like a simple homothetic rule related to
c the location of the residual cut-off with regard to the initial
c location of the edge of the friction cone
       ETf = ETi * RTS / TT0
c the behaviour in shear is assumed to follow a
c homothetic rule with respect to the behaviour in pure shear
c the coefficient of this homothetic transform is given by...
c     [1
       if ( JSP .eq. 0 ) then
c       [2
         if ( ETf .ne. ZERO ) then
c         [3
           if ( JCT .eq. 0 ) then
c                 \--> shear yielding with tension
             Chom = (ETf - Tft)/(ETf)*Dtensi
c           print *,'Coef. hom. ', Chom
c WARNING the original version used
c Chom = (RTS - Tft)/(RTS +Pnor)
c xmat(5) or PNOR is a simple numerical value used is order to
c avoid division by zero if the residual cohesion vanishes.
c           print *,'???', Dtensi, ONE
c           [4
             if ( Dtensi .eq. ONE) then
c                        \--> no tensile damage
c             [5
                 if ( abs(Tfs0) .gt. (RSS*Chom) ) then
c                                       \--> YIELDING IN SHEAR
c                print *,'YIELDING IN SHEAR'
c compute the total plastic shear strain...
                 Afps = Aips + ((abs(Tfs0)/Chom)-RSS)/Ks
c                  print *,'plasti shear strain ', Afps
c ...and the final cohesion residual
                 call yofxcu(Afps,tras,jtras,RSS,dRSS,ireturn)
c               [6
                 if (ireturn .ne. 0) return
c                                    \--> ERROR
c               6]
                 Tfs = sign(Chom*RSS , Tfs)
                           if (Tfs0 .eq. ZERO) then
                       Tfs1 = ZERO
                       Tfs2 = ZERO
                     else
                                   Tfs1= Tfs*Tfs1/Tfs0
                           Tfs2= Tfs*Tfs2/Tfs0
                     endif
c                    print *,'updated shear stress ',Tfs
c                    print *,'Tfs1 Tfs2 ',Tfs1,Tfs2
                 endif
c             5]
             else
c              \--> already got some damage due to tension
               Afpstmp = Aips+ (abs(Tfs0)/Chom-RSS)/(Ks*Dtensi)
c                print *, Afpstmp
c             [5
                 if ( Afpstmp .gt. ZERO ) then
                   call yofxcu(Afpstmp,tras,jtras,RSSTMP,dRSS,ireturn)
c               [6
                 if (ireturn .ne. 0) return
c                                    \--> ERROR
c               6]
c                print *,'RSSTMP', RSSTMP
c               [6
                 if ( abs(Tfs0) .gt. (Chom*RSSTMP) ) then
c                                       \--> YIELDING IN SHEAR
                   Afps = Afpstmp
                         print *,'YIELDING IN SHEAR 2'
                     RSS = RSSTMP
                   Tfs = sign(Chom*RSS , Tfs)
                     if (Tfs0 .eq. ZERO) then
                       Tfs1 = ZERO
                       Tfs2 = ZERO
                     else
                                   Tfs1= Tfs*Tfs1/Tfs0
                           Tfs2= Tfs*Tfs2/Tfs0
                     endif
                     print *,'updated shear stress ',Tfs
                     print *,'Tfs1 Tfs2 ',Tfs1,Tfs2
                   endif
c               6]
                 endif
c             5]
             endif
c           4]
           else
c          \--> check shear yielding in compression region
               Chom = (ETf - Tfc)/(ETf)*Dcompr
c              print *,'Coef. hom. ', Chom
c            print *,'???', Dcompr, ONE
c           [4
             if ( Dcompr .eq. ONE) then
c                        \--> no crushing damage
c             [5
                 if ( abs(Tfs0) .gt. (RSS*Chom) ) then
c                                       \--> YIELDING IN SHEAR
c                print *,'YIELDING IN SHEAR'
c compute the total plastic shear strain...
                 Afps = Aips + ((abs(Tfs0)/Chom)-RSS)/Ks
c                  print *,'plasti shear strain ', Afps
c ...and the final cohesion residual
                 call yofxcu(Afps,tras,jtras,RSS,dRSS,ireturn)
c               [6
                 if (ireturn .ne. 0) return
c                                    \--> ERROR
c               6]
                 Tfs = sign(Chom*RSS , Tfs)
                           if (Tfs0 .eq. ZERO) then
                       Tfs1 = ZERO
                       Tfs2 = ZERO
                     else
                                   Tfs1= Tfs*Tfs1/Tfs0
                           Tfs2= Tfs*Tfs2/Tfs0
                     endif
c                  print *,'updated shear stress ',Tfs
c                    print *,'Tfs1 Tfs2 ',Tfs1,Tfs2
                 endif
c             5]
             else
c              \--> already got some damage due to compression
               Afpstmp = Aips+ (abs(Tfs0)/Chom-RSS)/(Ks*Dcompr)
c                print *, Afpstmp
c             [5
                 if ( Afpstmp .gt. ZERO ) then
                   call yofxcu(Afpstmp,tras,jtras,RSSTMP,dRSS,ireturn)
c               [6
                 if (ireturn .ne. 0) return
c                                    \--> ERROR
c               6]
c                print *,'RSSTMP', RSSTMP
c               [6
                 if ( abs(Tfs0) .gt. (Chom*RSSTMP) ) then
c                                       \--> YIELDING IN SHEAR
                   Afps = Afpstmp
c                        print *,'YIELDING IN SHEAR 2'
                     RSS = RSSTMP
                   Tfs = sign(Chom*RSS , Tfs)
                                 if (Tfs0 .eq. ZERO) then
                       Tfs1 = ZERO
                       Tfs2 = ZERO
                     else
                                   Tfs1= Tfs*Tfs1/Tfs0
                           Tfs2= Tfs*Tfs2/Tfs0
                     endif
c                    print *,'updated shear stress ',Tfs
c                    print *,'Tfs1 Tfs2 ',Tfs1,Tfs2
                   endif
c               6]
                 endif
c             5]
             endif
c           4]
             endif
c         3]
         else
c        \--> the edge has reach the origin... no homothetic rule
c             can be computed !
           print *,'EDGE ZERO -- ERROR'
         endif
c       2]
       endif
c     1]
c______________________________________________________________________
c
c final stage --> updating the output variables...
c ... first the final stress...
c     [1
       if ( JCT .eq. 0 ) then
         sigf(Inorm) = Tft
           sigi(Inorm) = Tft
         else
           sigf(Inorm) = Tfc
           sigi(Inorm) = Tfc
         endif
c     1]
         sigf(Ishea1) = Tfs1
         sigi(Ishea1) = Tfs1
         sigf(Ishea2) = Tfs2
         sigi(Ishea2) = Tfs2
c ... then the final internal variables...
       varf(Ipress) = Afpt
         varf(Ishear) = Afps
         varf(Icompr) = Afpc
         varf(Itotal) = Aftt
c ... and finally, the plastic strain increments.
c     [1
       if ( JCT .eq. 0 ) then
         ddepsp(Inorm) = varf(Ipress) - vari(Ipress)
         else
           ddepsp(Inorm) = varf(Icompr) - vari(Icompr)
         endif
c     1]
         ddepsp(Ishea1) = varf(Ishear) - vari(Ishear)
       ddepsp(Ishea2) = varf(Ishear) - vari(Ishear)
         vari(Ipress) = Afpt
         vari(Ishear) = Afps
         vari(Icompr) = Afpc
         vari(Itotal) = Aftt
         depsp(Inorm)  = depsp(Inorm)  + ddepsp(Inorm)
       depsp(Ishea1) = depsp(Ishea1) + ddepsp(Ishea1)
       depsp(Ishea2) = depsp(Ishea2) + ddepsp(Ishea2)
c END OF MAIN LOOP
c    0]
      enddo
c        print *,'output variables :'
c        print *,' --> stress state - normal shear'
c        print *, sigf(Inorm), sigf(Ishea1), sigf(Ishea2)
c        print *,'internal variables - normal - comp - shear - total'
c        print *, varf(Ipress), varf(Icompr),varf(Ishear), varf(Itotal)
c        print *,'plastic strain - normal - shear'
c        print *, depsp(Inorm), depsp(Ishea1),depsp(Ishea2)
c  And, that's all folks !!!!
c       print *,'sjonl2 ended normaly...'
      end
 
 
 
 

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