sigeps56.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    sigeps56   ../engine/source/materials/mat/mat056/sigeps56.F
!||--- called by ------------------------------------------------------
!||    mulaw      ../engine/source/materials/mat_share/mulaw.F90
!||    mulaw8     ../engine/source/materials/mat_share/mulaw8.F90
!||--- calls      -----------------------------------------------------
!||    finter2    ../engine/source/tools/curve/vinter.F
!||    vinter     ../engine/source/tools/curve/vinter.F
!||    vinter2    ../engine/source/tools/curve/vinter.F
!||====================================================================
      SUBROUTINE sigeps56(
     1     NEL    ,NUPARAM,NUVAR   ,MFUNC   ,KFUNC   ,NPF    ,
     2     TF     ,TIME   ,TIMESTEP,UPARAM  ,RHO0    ,RHO    ,
     3     VOLUME ,EINT   ,
     4     EPSPXX ,EPSPYY ,EPSPZZ  ,EPSPXY  ,EPSPYZ  ,EPSPZX ,
     5     DEPSXX ,DEPSYY ,DEPSZZ  ,DEPSXY  ,DEPSYZ  ,DEPSZX ,
     6     EPSXX  ,EPSYY  ,EPSZZ   ,EPSXY   ,EPSYZ   ,EPSZX  ,
     7     SIGOXX ,SIGOYY ,SIGOZZ  ,SIGOXY  ,SIGOYZ  ,SIGOZX ,
     8     SIGNXX ,SIGNYY ,SIGNZZ  ,SIGNXY  ,SIGNYZ  ,SIGNZX ,
     9     SIGVXX ,SIGVYY ,SIGVZZ  ,SIGVXY  ,SIGVYZ  ,SIGVZX ,
     A     SOUNDSP,VISCMAX,UVAR    ,OFF     ,NGL     ,IPT    ,
     B     IPM    ,MAT    ,EPSP    ,IPLA    ,YLD     ,PLA    ,
     C    DPLA    ,AMU    )
CGW087     B     PM     ,MAT    )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C---------+---------+---+---+--------------------------------------------
C VAR     | SIZE    |TYP| RW| DEFINITION
C---------+---------+---+---+--------------------------------------------
C NEL     |  1      | I | R | SIZE OF THE ELEMENT GROUP NEL 
C NUPARAM |  1      | I | R | SIZE OF THE USER PARAMETER ARRAY
C NUVAR   |  1      | I | R | NUMBER OF USER ELEMENT VARIABLES
C---------+---------+---+---+--------------------------------------------
C MFUNC   |  1      | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
C KFUNC   | NFUNC   | I | R | FUNCTION INDEX not used
C NPF     |  *      | I | R | FUNCTION ARRAY   
C TF      |  *      | F | R | FUNCTION ARRAY 
C---------+---------+---+---+--------------------------------------------
C TIME    |  1      | F | R | CURRENT TIME
C TIMESTEP|  1      | F | R | CURRENT TIME STEP
C UPARAM  | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
C RHO0    | NEL     | F | R | INITIAL DENSITY
C RHO     | NEL     | F | R | DENSITY
C VOLUME  | NEL     | F | R | VOLUME
C EINT    | NEL     | F | R | TOTAL INTERNAL ENERGY
C EPSPXX  | NEL     | F | R | STRAIN RATE XX
C EPSPYY  | NEL     | F | R | STRAIN RATE YY
C ...     |         |   |   |
C DEPSXX  | NEL     | F | R | STRAIN INCREMENT XX
C DEPSYY  | NEL     | F | R | STRAIN INCREMENT YY
C ...     |         |   |   |
C EPSXX   | NEL     | F | R | STRAIN XX
C EPSYY   | NEL     | F | R | STRAIN YY
C ...     |         |   |   |
C SIGOXX  | NEL     | F | R | OLD ELASTO PLASTIC STRESS XX 
C SIGOYY  | NEL     | F | R | OLD ELASTO PLASTIC STRESS YY
C ...     |         |   |   |    
C---------+---------+---+---+--------------------------------------------
C SIGNXX  | NEL     | F | W | NEW ELASTO PLASTIC STRESS XX
C SIGNYY  | NEL     | F | W | NEW ELASTO PLASTIC STRESS YY
C ...     |         |   |   |
C SIGVXX  | NEL     | F | W | VISCOUS STRESS XX
C SIGVYY  | NEL     | F | W | VISCOUS STRESS YY
C ...     |         |   |   |
C SOUNDSP | NEL     | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
C VISCMAX | NEL     | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
C---------+---------+---+---+--------------------------------------------
C UVAR    |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
C OFF     | NEL     | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
C---------+---------+---+---+--------------------------------------------
#include      "param_c.inc"
#include      "parit_c.inc"
#include      "scr17_c.inc"
#include      "com01_c.inc"
#include      "com08_c.inc"
#include      "units_c.inc"
C
      INTEGER NEL, NUPARAM, NUVAR,IPT,IPLA,
     .   NGL(NEL),MAT(NEL),IPM(NPROPMI,*)
      my_real
     .   TIME,TIMESTEP,UPARAM(*),
     .   RHO(NEL),RHO0(NEL),VOLUME(NEL),EINT(NEL),
     .   EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
     .   EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
     .   DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
     .   DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
     .   EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
     .   EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
     .   SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
     .   sigoxy(nel),sigoyz(nel),sigozx(nel),
     .   epsp(nel),amu(nel)
C-----------------------------------------------
C   O U T P U T   A r g u m e n t s
C-----------------------------------------------
      my_real
     .    signxx(nel),signyy(nel),signzz(nel),
     .    signxy(nel),signyz(nel),signzx(nel),
     .    sigvxx(nel),sigvyy(nel),sigvzz(nel),
     .    sigvxy(nel),sigvyz(nel),sigvzx(nel),
     .    soundsp(nel),viscmax(nel),dpla(nel)
C-----------------------------------------------
C   I N P U T   O U T P U T   A r g u m e n t s 
C-----------------------------------------------
      my_real 
     .        uvar(nel,nuvar), off(nel),  yld(nel),pla(nel)
C-----------------------------------------------
C   VARIABLES FOR FUNCTION INTERPOLATION 
C-----------------------------------------------
      INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
      my_real
     .   FINTER2, TF(*)
      EXTERNAL FINTER2
C        Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
C        Y       : y = f(x)
C        X       : x
C        DYDX    : f'(x) = dy/dx
C        IFUNC(J): FUNCTION INDEX
C              J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
C        NPF,TF  : FUNCTION PARAMETER
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,IADBUFV,J1,J2,JJ(MVSIZ),NFUNC,
     .        nrate,ipos1(mvsiz),ipos2(mvsiz),iad1(mvsiz),
     .        ilen1(mvsiz),iad2(mvsiz),ilen2(mvsiz),
     .        ifunc(100),nfuncv,pfun,
     .        iposp(mvsiz),iadp(mvsiz),ilenp(mvsiz),nfuncm,nratem,
     .        ipflag,iparam,npar,
     .        nindx,indx(mvsiz),mx
      my_real
     .        e,nu,p,dav,vm,r,fac,epst,ep1,ep2,ep3,ep4,ep5,ep6,
     .        e1,e2,e3,e4,e5,e6,c,cc,d,y,yp,e42,e52,e62,epst2,
     .        c1,g,g2,g3,
     .        epsmax,rate(mvsiz,2),yfac(mvsiz,2),
     .        y1(mvsiz),y2(mvsiz),h(mvsiz),dydx1(mvsiz),
     .        dydx2(mvsiz),fail(mvsiz),epsr1,
     .        epsr2,p0(mvsiz),pfac(mvsiz),
     .        dfdp(mvsiz)
c     .        ONE,ZERO,EM20
c      DATA ZERO/0.0/,ONE/1.0/,EM20/1.E-20/
C-----------------------------------------------
C     USER VARIABLES INITIALIZATION
C-----------------------------------------------
      IF (ivector/=1) THEN
       mx = mat(1)
       nfunc  = ipm(10,mx)
       DO j=1,nfunc
         ifunc(j)=ipm(10+j,mx)
       ENDDO
      ELSE
        mx = mat(1)
        nfuncm = 0
        nfuncv = ipm(10,mx)
        nfuncm = max(nfuncm,nfuncv)
 
        DO j=1,nfuncm
          IF(nfuncv>=j) THEN
            ifunc(j)=ipm(10+j,mx)
          ENDIF
        ENDDO
      ENDIF
C
      nratem = 0
C
      mx = mat(1)
      iadbufv = ipm(7,mx)-1
      e   = uparam(iadbufv+2)
      g   = uparam(iadbufv+5)
      g2  = two*g
      g3  = three*g
      nu  = uparam(iadbufv+6)
      c1  = e/three/(one - two*nu)
C
      nrate = nint(uparam(iadbufv+1))
      nratem = max(nratem,nrate)
      epsmax=uparam(iadbufv+6+2*nrate+1)
      IF(epsmax==zero)THEN
        IF(tf(npf(ifunc(1)+1)-1)==zero)THEN
         epsmax=tf(npf(ifunc(1)+1)-2)
        ELSE
         epsmax= ep30
        ENDIF
      ENDIF
      epsr1 =uparam(iadbufv+6+2*nrate+2)
      IF(epsr1==zero)epsr1=ep30
      epsr2 =uparam(iadbufv+6+2*nrate+3)
      IF(epsr2==zero)epsr2=twoep30
      ipflag = 0
      mx = mat(1)
      npar   = ipm(9,mx)
      iparam = 15 + 2*nrate
      DO i=1,nel
        pfac(i) = one
        pfun = 0
        IF (npar>=iparam) THEN
          j = nrate+1
          ifunc(j)=ipm(10+j,mx)
          pfun=nint(uparam(iadbufv+iparam))
          IF (pfun>0) ipflag = ipflag + 1
         ENDIF
      ENDDO
C
      IF (isigi==0) THEN
      IF(time==zero)THEN
        DO i=1,nel
         uvar(i,1)=zero
         uvar(i,2)=zero
         DO j=1,nrate
           uvar(i,j+2)=zero
         ENDDO
         IF (pfun>0) uvar(i,nrate+5)=zero
        ENDDO
      ENDIF
      ENDIF
C-----------------------------------------------
C
      DO i=1,nel
C
        pla(i) = uvar(i,1)
        p0(i) = -(sigoxx(i)+sigoyy(i)+sigozz(i))*third
        dav = (depsxx(i)+depsyy(i)+depszz(i))*third
        signxx(i)=sigoxx(i)+p0(i)+g2*(depsxx(i)-dav)
        signyy(i)=sigoyy(i)+p0(i)+g2*(depsyy(i)-dav)
        signzz(i)=sigozz(i)+p0(i)+g2*(depszz(i)-dav)
        signxy(i)=sigoxy(i)+g *depsxy(i)
        signyz(i)=sigoyz(i)+g *depsyz(i)
        signzx(i)=sigozx(i)+g *depszx(i)
C
        soundsp(i) = sqrt((c1+four*g/three)/rho0(i))
        viscmax(i) = zero
        dpla(i) = zero       
      ENDDO
C-------------------
C     STRAIN RATE
C-------------------
      DO i=1,nel
C  087        DAV = (EPSPXX(I)+EPSPYY(I)+EPSPZZ(I))/3.
C  087        E1 = EPSPXX(I) - DAV
C  087        E2 = EPSPYY(I) - DAV
C  087        E3 = EPSPZZ(I) - DAV
C  087        E4 = 0.5*EPSPXY(I)
C  087        E5 = 0.5*EPSPYZ(I)
C  087        E6 = 0.5*EPSPZX(I)
C  087        EPSP(I) =.5*(E1**2+E2**2+E3**2) +E4**2+E5**2+E6**2
C  087        EPSP(I) = SQRT(3.*EPSP(I))/1.5
C-------------------
C     STRAIN principal 1, 4 newton iterations 
C-------------------
        dav = (epsxx(i)+epsyy(i)+epszz(i)) * third
        e1 = epsxx(i) - dav
        e2 = epsyy(i) - dav
        e3 = epszz(i) - dav
        e4 = half*epsxy(i)
        e5 = half*epsyz(i)
        e6 = half*epszx(i)
C        -y = (e1-x)(e2-x)(e3-x)
C           - e5^2(e1-x) - e6^2(e2-x) - e4^2(e3-x)
C           + 2e4 e5 e6
C         e1 + e2 + e3 = 0 => terme en x^2 = 0
C         y = x^3 + c x + d
c         yp= 3 x^2 + c
        e42 = e4*e4
        e52 = e5*e5
        e62 = e6*e6
        c = - e1*e1 - e2*e2 - e3*e3 - e42 - e52 - e62
        d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42
     &      - two*e4*e5*e6 
        cc = c*third
        epst = sqrt(-cc)
        epst2 = epst * epst
        y = (epst2 + c)* epst + d
        IF(abs(y)>em8)THEN
          epst = onep75 * epst
          epst2 = epst * epst
          y = (epst2 + c)* epst + d
          yp = three*epst2 + c
          IF(yp/=zero)epst = epst - y/yp
          epst2 = epst * epst
          y = (epst2 + c)* epst + d
          yp = three*epst2 + c
          IF(yp/=zero)epst = epst - y/yp
          epst2 = epst * epst
          y = (epst2 + c)* epst + d
          yp = three*epst2 + c
          IF(yp/=zero)epst = epst - y/yp
          epst2 = epst * epst
          y = (epst2 + c)* epst + d
          yp = three*epst2 + c
          IF(yp/=zero)epst = epst - y/yp
          epst = epst + dav
        ENDIF
C-------------------
C     tension failure 
C-------------------
        fail(i) = max(zero,min(one,
     .            (epsr2-epst)/(epsr2-epsr1) ))
C
c jbm023
c
c        FAIL(I)=1.
      ENDDO
C-------------------
C     CRITERE
C-------------------
      DO i=1,nel
            jj(i) = 1
      ENDDO
      IF (ivector/=1) THEN
       DO i=1,nel
        DO j=2,nrate-1
          IF(epsp(i)>=uparam(iadbufv+6+j)) jj(i) = j
        ENDDO
       ENDDO
      ELSE
        DO j=2,nratem-1
          DO i=1,nel
            IF(nrate-1>=j) THEN
              IF(epsp(i)>=uparam(iadbufv+6+j)) jj(i) = j
            ENDIF
          ENDDO
        ENDDO
      ENDIF
C
      DO i=1,nel
        rate(i,1)=uparam(iadbufv+6+jj(i))
        rate(i,2)=uparam(iadbufv+6+jj(i)+1)
        yfac(i,1)=uparam(iadbufv+6+nrate+jj(i))
        yfac(i,2)=uparam(iadbufv+6+nrate+jj(i)+1)
      ENDDO
C
      DO i=1,nel
        j1 = jj(i)
        j2 = j1+1
C  037        IPOS1(I) = NINT(UVAR(I,J1))
        ipos1(i) = nint(uvar(i,j1+2))
        iad1(i)  = npf(ifunc(j1)) / 2 + 1
        ilen1(i) = npf(ifunc(j1)+1) / 2 - iad1(i) - ipos1(i)
C  037        IPOS2(I) = NINT(UVAR(I,J2))
        ipos2(i) = nint(uvar(i,j2+2))
        iad2(i)  = npf(ifunc(j2)) / 2 + 1
        ilen2(i) = npf(ifunc(j2)+1) / 2 - iad2(i) - ipos2(i)
      ENDDO
C
      CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
      CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
C
C---  PRESSURE DEPENDENT YIELD FUNCTION FACTOR
      IF (ipflag==nel.AND.(iparit==0)) THEN
C------- to optimize in SPMD & PARIT/ON case add a global flag 
C        in a case of homogeneous element groups (to do if needed)
        DO i=1,nel
          iposp(i) = nint(uvar(i,nrate+5))
          iadp(i)  = npf(ifunc(pfun)) / 2 + 1
          ilenp(i) = npf(ifunc(pfun)+1) / 2 - iadp(i) - iposp(i)
        ENDDO
        CALL vinter2(tf,iadp,iposp,ilenp,nel,p0,dfdp,pfac)
        DO i=1,nel
          uvar(i,nrate+5) = iposp(i)
        ENDDO
      ELSEIF (ipflag>0) THEN
        DO i=1,nel
          IF (pfun>0) THEN
            iposp(i) = nint(uvar(i,nrate+5))
            iadp(i)  = npf(ifunc(pfun)) / 2 + 1
            ilenp(i) = npf(ifunc(pfun)+1) / 2 
     .                           - iadp(i) - iposp(i)
            pfac(i)  = finter2(tf   ,iadp(i),iposp(i),ilenp(i),
     .                         p0(i),dfdp(i))
            uvar(i,nrate+5) = iposp(i)
          ENDIF
        ENDDO
      ENDIF
      DO i=1,nel
        j1 = jj(i)
        j2 = j1+1
        y1(i)=y1(i)*yfac(i,1)
        y2(i)=y2(i)*yfac(i,2)
        fac   = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
        yld(i) = fail(i)*(y1(i)    + fac*(y2(i)-y1(i)))
        dydx1(i)=dydx1(i)*yfac(i,1)
        dydx2(i)=dydx2(i)*yfac(i,2)
        h(i)   = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
        h(i) = h(i) * max(zero, pfac(i))
        yld(i) = yld(i) * max(zero, pfac(i))
 
C  037        UVAR(I,J1) = IPOS1(I)
C  037        UVAR(I,J2) = IPOS2(I)
        uvar(i,j1+2) = ipos1(i)
        uvar(i,j2+2) = ipos2(i)
      ENDDO
C-------------------
C     PROJECTION
C-------------------
      IF(ipla==0)THEN
       DO i=1,nel
        vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
     1          +signxy(i)**2+signyz(i)**2+signzx(i)**2
        vm =sqrt(three*vm)
        r = min(one,yld(i)/ max(vm,em20))
c        P = C1 * (RHO(I)/RHO0(I)- ONE)
        p = c1 * amu(i)
        signxx(i)=signxx(i)*r-p
        signyy(i)=signyy(i)*r-p
        signzz(i)=signzz(i)*r-p
        signxy(i)=signxy(i)*r
        signyz(i)=signyz(i)*r
        signzx(i)=signzx(i)*r
        pla(i)=pla(i) + (one - r)*vm/max(g3+h(i),em20)
        uvar(i,1)=pla(i)
        dpla(i) = (one - r)*vm/max(g3+h(i),em20)        
       ENDDO
      ELSEIF(ipla==2)THEN
       DO i=1,nel
        vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
     1          +signxy(i)**2+signyz(i)**2+signzx(i)**2
        vm =sqrt(three*vm)
        r = min(one,yld(i)/ max(vm,em20))
c        P = C1 * (RHO(I)/RHO0(I)-ONE)
        p = c1 * amu(i)
        signxx(i)=signxx(i)*r-p
        signyy(i)=signyy(i)*r-p
        signzz(i)=signzz(i)*r-p
        signxy(i)=signxy(i)*r
        signyz(i)=signyz(i)*r
        signzx(i)=signzx(i)*r
        pla(i)=pla(i) + (one-r)*vm/max(g3,em20)
        uvar(i,1)=pla(i)
        dpla(i) = (one-r)*vm/max(g3,em20)        
       ENDDO
      ELSEIF(ipla==1)THEN
       DO i=1,nel
        vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
     1          +signxy(i)**2+signyz(i)**2+signzx(i)**2
        vm =sqrt(three*vm)
        r = min(one,yld(i)/ max(vm,em20))
C       plastic strain increment.
        dpla(i)=(one - r)*vm/max(g3+h(i),em20)
C       actual yield stress.
        yld(i)=max(yld(i)+dpla(i)*h(i),zero)
        r = min(one,yld(i)/ max(vm,em20))
C
c        P = C1 * (RHO(I)/RHO0(I)-ONE)
        p = c1 * amu(i)
        signxx(i)=signxx(i)*r-p
        signyy(i)=signyy(i)*r-p
        signzz(i)=signzz(i)*r-p
        signxy(i)=signxy(i)*r
        signyz(i)=signyz(i)*r
        signzx(i)=signzx(i)*r
       pla(i)=pla(i) + dpla(i)
        uvar(i,1)=pla(i)
       ENDDO
      ENDIF
C-------------------
      DO i=1,nel
        IF(off(i)<em01) off(i)=zero
        IF(off(i)<one) off(i)=off(i)*four_over_5
      ENDDO
C
      nindx=0
      DO i=1,nel
        IF(pla(i)>epsmax.AND.off(i)==one) THEN
          off(i)=four_over_5
          nindx=nindx+1
          indx(nindx)=i
        ENDIF
      ENDDO
      IF(nindx>0)THEN
        DO j=1,nindx
#include "lockon.inc"
        WRITE(iout, 1000) ngl(indx(j))
        WRITE(istdo,1100) ngl(indx(j)),tt
#include "lockoff.inc"
        ENDDO
      ENDIF
C
 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
 1100 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10,
     .          ' AT TIME :',g11.4)
C
      RETURN
      END
C