Sources Esope | Cast3M

dsigrank
C DSIGRANK  SOURCE    FANDEUR   14/03/25    21:15:10     7993
      SUBROUTINE dsigRank(sigma,dgt,i3,i6,parahot3,idimpara3,
     .                    rkappait,lcp)
 
c     This subroutine calculates the derivative of the Rankine
c     plastic function with respect to sigma
 
c     gt = ft (because associated plasticity)
c        = S1 = largest principal stress
 
      IMPLICIT REAL*8 (A-B,D-H,O-Z)
      implicit integer (I-K,M,N)
      implicit logical (L)
      implicit character*10 (C)
      character*20 cloc1,cloc2
 
      dimension sigma(i6),dgt(i6)
      dimension rcossigmapr(3,3),parahot3(idimpara3)
      dimension P1(3,3),pi1(3),vloc(3),vloc1(1)
 
*****
*     MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS)
      IF(I3.GT.3) PRINT *, ' DSIGRANK - ERREUR  I3 = ', I3, ' > 3 '
*****
      i1 = 1
      i2 = 2
      i4 = 4
      i5 = 5
      r1 = 1.d0
      r2 = 2.d0
      r4 = 4.d0
 
      if (lcp) then
c       in plane stress, calculation of the derivative of the principal stress
        sigmax  = sigma(1)
        sigmay  = sigma(2)
        sigmaxy = sigma(4)
        diff = sigmax-sigmay
        racine = diff*diff + r4*sigmaxy*sigmaxy
        cloc1 = 'racine'
        cloc2 = 'dsigrank'
        call check_sqrt(racine,cloc1,cloc2,lerror)
c       **** **********
        if (Lerror) then
          return
        endif
        racine = sqrt(racine)
        rtest = max(sigmax,sigmay)
        if (racine.gt.abs(rtest)*1.e-6) then
            dgt(1) = ( diff/racine+r1)/r2
            dgt(2) = (-diff/racine+r1)/r2
            dgt(4) = r2*sigmaxy/racine
        else
c       APEX of the Rankine surface: sigma1 = sigma2 ; sigma12 = 0
          dgt(1) = sqrt(r2)/r2
          dgt(2) = dgt(1)
          dgt(4) = 0.D0
        endif
        dgt(3) = 0.d0
        dgt(5) = 0.d0
        dgt(6) = 0.d0
      else
c       in full 3D, calculation of the derivative of the principal stress
        call sigmapr3ETC(sigma,S1,S2,S3,rcossigmapr,lcp)
c            =============
c       algorithm to find the eigen values and eigen vectors of symmetric matrix sigma
c       APEX test
c       Carefull with precision of the real*8
        sig1 = sigma(1)/10.d6
        sig2 = sigma(2)/10.d6
        sig3 = sigma(3)/10.d6
        sig4 = sigma(4)/10.d6
        sig5 = sigma(5)/10.d6
        sig6 = sigma(6)/10.d6
        rapex = abs(sigma(1))+abs(sigma(2))+abs(sigma(3)) -
     .          3.*(abs(sigma(1)*sigma(2)*sigma(3))**(1.d0/3.d0))
        rapex = rapex + 9.*(abs(sigma(4))+abs(sigma(5))+
     .                 abs(sigma(6)))
        if (rapex.le.3.d-4) then
c         APEX of the Rankine surface
          dgt(1)=sqrt(3.d0)/3.d0
          dgt(2)=dgt(1)
          dgt(3)=dgt(1)
          dgt(4)=0.d0
          dgt(5)=0.d0
          dgt(6)=0.d0
        else
c         Calculation of dgt/dsig = dgt/ds1 * ds1/dsig
c         dgt/dsig is the first column of the transposed matrix of axis change
c             between global axis and local axis of the principal stresses
          rkx = rcossigmapr(i1,i1)
          rmx = rcossigmapr(i2,i1)
          rnx = rcossigmapr(i3,i1)
          rky = rcossigmapr(i1,i2)
          rmy = rcossigmapr(i2,i2)
          rny = rcossigmapr(i3,i2)
          rkz = rcossigmapr(i1,i3)
          rmz = rcossigmapr(i2,i3)
          rnz = rcossigmapr(i3,i3)
 
c         Sxx, Syy, Szz, Sxy, Sxz, Syz
          dgt(i1) = rkx*rkx
          dgt(i2) = rmx*rmx
          dgt(i3) = rnx*rnx
          dgt(i4) = 2.0d0*rkx*rmx
          dgt(i5) = 2.0d0*rkx*rnx
          dgt(i6) = 2.0d0*rmx*rnx
        endif
      endif
 
      RETURN
      END