Télécharger sjonl2.eso

Retour à la liste

sjonl2
C SJONL2    SOURCE    PASCAL    21/09/24    21:15:01     11108          
      SUBROUTINE SJONL2(sigi, depst, vari, xmat, trac, jtrac ,
     >                  tras, jtras, trat, jtrat, sigf, varf,
     >                  depsp, ireturn)
c --- SJONL2 ---
c computing plastic flow in joint element
c under 2D plane (strain) conditions
c
c former version 13-06-00 ylp.
c update         15-03-01 ylp.
c                a shear unloading parameter has been added --> beta
*
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 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(2), ddepsp(2)
         dimension xmat(*)
c xmat 1 --> Ks     elastic shear modulus
c      2 --> Kt     elastic traction modulus
c      3 --> rho?   mass per unit volume
c      4 --> alpha? thermal coefficient
c      5 --> edge   location of the edge of the friction cone
c      6 --> cplg   make damage dependencies
c      7 --> beta   control shear unloading
         dimension trac(2,jtrac), tras(2,jtras), trat(2,jtrat)
c 1 --> strain
c 2 --> stress
c WARNING original version
c dimension TRAS(2,NTRAT), TRAT(2,NTRAS) ?????
c_______________________________________________________________________
c
*       save iloop
*       data iloop/0/
*       iloop = iloop + 1
c       do i=1,2
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
c look at xmat(6) and deduce elastic/plastic...
c
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
         Icompr  = 4
         Istrain = 1
         Istress = 2
c unloading constitutive parameters
c the unloading shear stiffness is computed by using a linear
c interpolation between the elastic stiffness and the
c damaged stiffness... Beta is equal to 0 to obtain a
c simple elastic unloading and is equal to 1 to get the
c damaged stiffness
       Beta = xmat(7)
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(Ipress) = depst(Ipress)/INCR
       ddepst(Ishear) = depst(Ishear)/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
c     [1
c       if ( JTP .ne. 0 ) then
         Aipt = vari(Ipress)
c          \--> initial plastic tensile strain
c       print *,'initial plastic tensile strain',Aipt
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 check if Aipt has a positive value or zero
 
c       [2
         if ( ireturn .ne. 0 ) return
c       2]
c                            \--> ERROR
c       print *,'tensile stress RTS', RTS
c       endif
c     1]
c     [1
c       if ( JSP .ne. 0 ) then
           Aips = vari(Ishear)
c          \--> initial plastic shear strain
c       print *,'initial plastic shear strain', Aips
c compute the initial cohesion residual
c         print *,tras(Istress,1)
         call yofxcu(Aips, tras, jtras, RSS, dRSS, ireturn)
c       [2
           if (ireturn .ne. 0) return
c                             \--> ERROR
c       2]
c       print *,'initial cohesion residual',RSS
c       endif
c     1]
c       if ( JCP .eq. 0 ) then
           Aipc = vari(Icompr)
c          \--> initial plastic compression strain
c compute the initial compressive strength
         call yofxcu(Aipc, trac, jtrac, RCS, dRCS, ireturn)
c       print *,'Aipc RCS ',Aipc, RCS
c       [2
           if (ireturn .ne. 0) return
c                             \--> ERROR
c       2]
c       endif
c     1]
c______________________________________________________________________
c
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
       Aitt = vari(Itotal)
c        \--> initial total tensile strain
c compute final total tensile strain...
       Aftt = Aitt + ddepst(Ipress)
*      write(6,*) 'Aitt, Aftt =',Aitt, Aftt
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
** mise 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'
*            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     [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 damage 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 )
c                 \--> damaged stiffness
             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 damage 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...
c          print *,ubound(trac)
         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
             Ksd = RSS / (Aips + TS0 / Ks0)
c         [3
           if ( RSS .eq. TS0 ) then
c                            \--> still in the elastic regime
             Ks = Ksd
             else
c             \--> non linear regime...
c                  the stiffness is computed by using an
c                  interpolation between the elastic stiffness
c                  and the damaged stiffness
               Ks  = ( ( 1.0D0 - Beta ) * Ks0 ) + ( Beta * Ksd )
             endif
c         3]
           endif
c update the shear damage variable
         if (JST .eq. 0) Dshear = RSS / TS0
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(Ipress) + ( ddepst(Ipress) * Kt * Dshear)
            endif
c        2]
          Tfs = sigi(Ishear) + ( ddepst(Ishear) * Ks * Dtensi)
       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(Ipress) + ( ddepst(Ipress) * Kc )
          endif
c        2]
          Tfs = sigi(Ishear) + ( ddepst(Ishear) * Ks * Dcompr)
         endif
c     1]
 
c       print *,'stress state elastic prediction'
c       print *,'tension    compression    shear'
c       print *,Tft, Tfc, Tfs
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
c      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(Tfs) .gt. (RSS*Chom) ) then
c                                       \--> YIELDING IN SHEAR
c                print *,'YIELDING IN SHEAR'
c compute the total plastic shear strain...
                 Afps = Aips + ((abs(Tfs)/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)
c                  print *,'updated shear stress ',Tfs
                 endif
c             5]
             else
c              \--> already got some damage due to tension
               Afpstmp = Aips+ (abs(Tfs)/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(Tfs) .gt. (Chom*RSSTMP) ) then
c                                       \--> YIELDING IN SHEAR
                   Afps = Afpstmp
c                        print *,'YIELDING IN SHEAR 2'
                     RSS = RSSTMP
                   Tfs = sign(Chom*RSS , Tfs)
c                    print *,'updated shear stress ',Tfs
                   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(Tfs) .gt. (RSS*Chom) ) then
c                                       \--> YIELDING IN SHEAR
c                print *,'YIELDING IN SHEAR'
c compute the total plastic shear strain...
                 Afps = Aips + ((abs(Tfs)/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)
c                  print *,'updated shear stress ',Tfs
                 endif
c             5]
             else
c              \--> already got some damage due to compression
               Afpstmp = Aips+ (abs(Tfs)/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(Tfs) .gt. (Chom*RSSTMP) ) then
c                                       \--> YIELDING IN SHEAR
                   Afps = Afpstmp
c                        print *,'YIELDING IN SHEAR 2'
                     RSS = RSSTMP
                   Tfs = sign(Chom*RSS , Tfs)
c                    print *,'updated shear stress ',Tfs
                   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 !
c           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(Ipress) = Tft
           sigi(Ipress) = Tft
         else
           sigf(Ipress) = Tfc
           sigi(Ipress) = Tfc
         endif
c     1]
         sigf(Ishear) = Tfs
         sigi(Ishear) = Tfs
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(Ipress) = varf(Ipress) - vari(Ipress)
         else
           ddepsp(Ipress) = varf(Icompr) - vari(Icompr)
         endif
c     1]
         ddepsp(Ishear) = varf(Ishear) - vari(Ishear)
         vari(Ipress) = Afpt
         vari(Ishear) = Afps
         vari(Icompr) = Afpc
         vari(Itotal) = Aftt
         depsp(Ipress) = depsp(Ipress) + ddepsp(Ipress)
       depsp(Ishear) = depsp(Ishear) + ddepsp(Ishear)
c END OF MAIN LOOP
c    0]
      enddo
c        print *,'output variables :'
c        print *,' --> stress state - normal shear'
c        print *, sigf(Ipress), sigf(Ishear)
c        print *,'internal variables - normal - shear - total'
c        print *, varf(Ipress), varf(Ishear), varf(Itotal)
c        print *,'plastic strain - normal - shear'
c        print *, depsp(Ipress), depsp(Ishear)
c  And, that's all folks !!!!
c       print *,'sjonl2 ended normaly...'
      end
 
 
 
 
 
 

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