sigeps48.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/.
!||====================================================================
!||    sigeps48   ../engine/source/materials/mat/mat048/sigeps48.F
!||--- called by ------------------------------------------------------
!||    mulaw      ../engine/source/materials/mat_share/mulaw.F90
!||    mulaw8     ../engine/source/materials/mat_share/mulaw8.F90
!||====================================================================
      SUBROUTINE sigeps48(
     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    ,SIGY    ,DEFP    ,DPLA1  ,
     C     AMU    )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.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      "scr17_c.inc"
C
      INTEGER NEL, NUPARAM, NUVAR,IPT,
     .   IPM(NPROPMI,*),NGL(NEL),MAT(NEL)
      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),sigy(nel),defp(nel),
     .    dpla1(mvsiz)
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)
C-------------------------
C    Variables not used (user functions)
C-------------------------
      INTEGER NPF(*),MFUNC,KFUNC(MFUNC)
      my_real
     . TF(*)
C      EXTERNAL FINTER
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,IADBUF,MX
      my_real 
     .        E,NU,P,DAV,VM,R,DPLA,EPST,E1,E2,E3,E4,E5,E6,C,CD,D,
     .        Y,YP,E42,E52,E62,EPST2,PA,PB,PC,PDA,PDB,YY,
     .        C1(MVSIZ),C2(MVSIZ),C3(MVSIZ),
     .        g(mvsiz),g2(mvsiz),g3(mvsiz),epsgm(mvsiz),
     .        epmax(mvsiz),smax(mvsiz),pla(mvsiz),fail(mvsiz),
     .        epsr1(mvsiz),epsr2(mvsiz),yld(mvsiz),h(mvsiz),
     .        ca(mvsiz),cb(mvsiz),cn(mvsiz),cc(mvsiz),co(mvsiz),
     .        cm(mvsiz),ce(mvsiz),ck(mvsiz),eps0(mvsiz)
C==================================================================|
      IF(time==zero)THEN
        DO i=1,nel
          uvar(i,1)=zero
        ENDDO
      ENDIF
C
C---
      mx = mat(1)
      iadbuf  = ipm(7,mx)-1
      DO i=1,nel
        e       = uparam(iadbuf+1)
        nu      = uparam(iadbuf+2)
        ca(i)   = uparam(iadbuf+3)
        smax(i) = uparam(iadbuf+4)
        epmax(i)= uparam(iadbuf+5)
        epsr1(i)= uparam(iadbuf+6)
        epsr2(i)= uparam(iadbuf+7)
        cb(i)   = uparam(iadbuf+8)
        cn(i)   = uparam(iadbuf+9)
        cc(i)   = uparam(iadbuf+10)
        co(i)   = uparam(iadbuf+11)
        cm(i)   = uparam(iadbuf+12)
        g(i)    = uparam(iadbuf+16)
        g2(i)   = uparam(iadbuf+17)
        g3(i)   = uparam(iadbuf+18)
        c1(i)   = uparam(iadbuf+19)
        c2(i)   = uparam(iadbuf+20)
        c3(i)   = uparam(iadbuf+21)
        epsgm(i)= uparam(iadbuf+22)
        eps0(i) = uparam(iadbuf+23)
        ce(i)   = uparam(iadbuf+24)
        ck(i)   = uparam(iadbuf+25)
      ENDDO
C
C---   DEVIATORIC STRESS ESTIMATE
C
      DO i=1,nel
        pla(i) = uvar(i,1)
        p = -(sigoxx(i)+sigoyy(i)+sigozz(i))* third
        dav = (depsxx(i)+depsyy(i)+depszz(i))* third
        signxx(i)=sigoxx(i) + p + g2(i)*(depsxx(i)-dav)
        signyy(i)=sigoyy(i) + p + g2(i)*(depsyy(i)-dav)
        signzz(i)=sigozz(i) + p + g2(i)*(depszz(i)-dav)
        signxy(i)=sigoxy(i) + g(i)*depsxy(i)
        signyz(i)=sigoyz(i) + g(i)*depsyz(i)
        signzx(i)=sigozx(i) + g(i)*depszx(i)
C
        soundsp(i) = sqrt((c1(i)+four*g(i)/three)/rho0(i))
        viscmax(i) = zero
         dpla1(i) = zero                
      ENDDO
C__________
C
C---    STRAIN principal 1, 4 newton iterations 
C
      DO i=1,nel
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 
        cd = c*third
        epst = sqrt(-cd)
        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 = 3*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(i)-epst)/(epsr2(i)-epsr1(i)) ))
C
      ENDDO
C
C--- STRAIN RATE EFFECT, CURRENT YIELD & HARDENING
C---------------------
C Law formulation:
C sig = ca*SigY + cb*epla^cn + (cc-co*epla^cm)*ln(epsp/eps0) + ce*epsp^ck
C---------------------
      DO i=1,nel
        IF(pla(i)<=zero) THEN
         pa=ca(i)
        ELSEIF(pla(i)>epmax(i)) THEN
         pa=ca(i)+cb(i)*epmax(i)**cn(i)
        ELSE
         pa=ca(i)+cb(i)*pla(i)**cn(i)
        ENDIF
        IF(epsp(i)<=eps0(i)) THEN
         pb=zero
        ELSEIF(pla(i)<=zero) THEN
         pb=cc(i)*log(epsp(i)/eps0(i))
        ELSE
         pb=(cc(i)-co(i)*pla(i)**cm(i))*log(epsp(i)/eps0(i))
        ENDIF
        IF(epsp(i)<=zero) THEN
         pc=zero
        ELSE
         pc=ce(i)*epsp(i)**ck(i)
        ENDIF
C-----
        IF(pla(i)>zero. and .cn(i)>=one) THEN
         pda = cb(i)*cn(i)*pla(i)**(cn(i)- one)
        ELSEIF(pla(i)>zero. and .cn(i)<one)THEN
         pda = cb(i)*cn(i)*pla(i)**(one-cn(i))
        ELSE
         pda = e
        ENDIF
        IF(pla(i)<=zero. or .epsp(i)<=eps0(i)) THEN
         pdb = zero
        ELSEIF(cm(i)>=one) THEN
         pdb = co(i)*cm(i)*pla(i)**(cm(i)- one)*log(epsp(i)/eps0(i))
        ELSE
         pdb = co(i)*cm(i)*pla(i)**(one - cm(i))*log(epsp(i)/eps0(i))
        ENDIF
C-----
        yy    = pa + pb + pc
        yld(i)= min(smax(i)+pc, yy)
        h(i) = pda + pdb
        IF (yld(i)<yy) h(i)  = zero
        yld(i)= fail(i)*yld(i)
        h(i)  = fail(i)*h(i)
        sigy(i)=yld(i)
        IF(pla(i)>epmax(i)) yld(i)=zero
      ENDDO        
C
C-------------------
C     VON MISES & RADIAL RETURN
C-------------------
      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=(one -r)*vm/max(g3(i)+h(i),em20)
C       actual yield stress.
        yld(i)=yld(i)+dpla*h(i)
        IF(pla(i)>epmax(i)) yld(i)=zero
        r = min(one,yld(i)/ max(vm,em20))
C
c        P = C1(I) * (RHO(I)/RHO0(I)- ONE)
        p = c1(i) * 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
        uvar(i,1)=pla(i)
        dpla1(i) = dpla        
      ENDDO
C
      DO i=1,nel
        sigy(i)=max(sigy(i),yld(i))
        defp(i)=uvar(i,1)
      ENDDO
C________________________________________________________
      RETURN
      END
C