sigeps28.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/.
!||====================================================================
!||    sigeps28    ../engine/source/materials/mat/mat028/sigeps28.F
!||--- called by ------------------------------------------------------
!||    mulaw       ../engine/source/materials/mat_share/mulaw.F90
!||    mulaw8      ../engine/source/materials/mat_share/mulaw8.F90
!||--- calls      -----------------------------------------------------
!||    finter      ../engine/source/tools/curve/finter.F
!||    vinter2     ../engine/source/tools/curve/vinter.F
!||    vinter2dp   ../engine/source/tools/curve/vinter.F
!||====================================================================
      SUBROUTINE sigeps28 (
     1     NEL    ,NUPARAM,NUVAR   ,NFUNC   ,IFUNC   ,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    ,AMU    )
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 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 NFUNC   |  1      | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
C IFUNC   | NFUNC   | I | R | FUNCTION INDEX 
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---------+---------+---+---+--------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
#include      "param_c.inc"
#include      "scr17_c.inc"
#include      "scr05_c.inc"
#include      "com08_c.inc"
#include      "units_c.inc"
C-----------------------------------------------
C   I N P U T   A r g u m e n t s
C-----------------------------------------------
C
      INTEGER NEL, NUPARAM, NUVAR,IPT,
     .   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),
     .   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)
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 FOR FUNCTION INTERPOLATION 
C-----------------------------------------------
      INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
      my_real
     .        FINTER ,TF(*)
      EXTERNAL FINTER
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,IADBUF,IF1,IF2,AUX,IC,II,K,
     .        IAD1(MVSIZ),IPOS1(MVSIZ),ILEN1(MVSIZ),
     .        INDEX(MVSIZ),
     .        NINDX,INDX(MVSIZ)
      my_real
     .        E11,E22,E33,G12,G23,G31,
     .        Y11,Y22,Y33,Y12,Y23,Y31,EP1,EP2,EP3,EP4,EP5,EP6,
     .        YC(MVSIZ),FAC1(MVSIZ),FAC2(MVSIZ),
     .        FAC3(MVSIZ),FAC4(MVSIZ),FAC5(MVSIZ),FAC6(MVSIZ),
     .        DYDX,DYDXV(MVSIZ),EP(MVSIZ,6),EPC(MVSIZ),
     .        EMX11,EMX22,EMX33,EMX12,EMX23,EMX31,AMUV
C=======================================================================
C     USER VARIABLES INITIALIZATION
C-----------------------------------------------
      IF(time==zero)THEN
       IF (nuvar>0) THEN
        DO 100 j=1,nuvar
        DO 100 i=1,nel
           uvar(i,j)=zero
 100    CONTINUE
       ENDIF
      ENDIF
C-----------------------------------------------
      DO i=1,nel
C
        iadbuf = ipm(7,mat(i))
        e11 = uparam(iadbuf)
        e22 = uparam(iadbuf+1)
        e33 = uparam(iadbuf+2)
        g12 = uparam(iadbuf+3)
        g23 = uparam(iadbuf+4)
        g31 = uparam(iadbuf+5)
C
        signxx(i) = sigoxx(i) + e11 * depsxx(i)
        signyy(i) = sigoyy(i) + e22 * depsyy(i)
        signzz(i) = sigozz(i) + e33 * depszz(i)
        signxy(i) = sigoxy(i) + g12 * depsxy(i)
        signyz(i) = sigoyz(i) + g23 * depsyz(i)
        signzx(i) = sigozx(i) + g31 * depszx(i)
C
        soundsp(i) = sqrt(max(e11,e22,e33,g12,g23,g31)/rho0(i))
        viscmax(i) = zero
C
      ENDDO
 
      nindx=0
      DO i=1,nel
        iadbuf = ipm(7,mat(i)) - 1  
        emx11 = uparam(iadbuf+9)    
        emx22 = uparam(iadbuf+10)   
        emx33 = uparam(iadbuf+11)   
        emx12 = uparam(iadbuf+12)   
        emx23 = uparam(iadbuf+13)   
        emx31 = uparam(iadbuf+14)
        fac1(i) = uparam(iadbuf+15) 
        fac2(i) = uparam(iadbuf+16) 
        fac3(i) = uparam(iadbuf+17) 
        fac4(i) = uparam(iadbuf+18) 
        fac5(i) = uparam(iadbuf+19) 
        fac6(i) = uparam(iadbuf+20) 
        IF((epsxx(i)>emx11.OR.                                  
     .      epsyy(i)>emx22.OR.                                  
     .      epszz(i)>emx33.OR.                                  
     .      abs(epsxy(i)/two)>emx12.OR.                             
     .      abs(epsyz(i)/two)>emx23.OR.                             
     .      abs(epszx(i)/two)>emx31).AND.off(i)/=zero) THEN       
          off(i) = zero                                            
          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
      DO i=1,nel
C
c        AMU = RHO(I)/RHO0(I) - ONE ! Change using AMU(I) to take into account THERM_STRESS
        ep(i,1) = amu(i)
        ep(i,2) = amu(i)
        ep(i,3) = amu(i)
        ep(i,4) = amu(i)
        ep(i,5) = amu(i)
        ep(i,6) = amu(i)
        iadbuf = ipm(7,mat(i))     
        if1=nint(uparam(iadbuf+6)) 
        if2=nint(uparam(iadbuf+7)) 
        IF(if1==1)THEN           
          ep(i,1) = epsxx(i)       
          ep(i,2) = epsyy(i)       
          ep(i,3) = epszz(i)       
        ELSEIF(if1==-1)THEN      
          ep(i,1) = -epsxx(i)      
          ep(i,2) = -epsyy(i)      
          ep(i,3) = -epszz(i)      
        ENDIF                      
        IF(if2==1)THEN           
          ep(i,4) = epsxy(i)       
          ep(i,5) = epsyz(i)       
          ep(i,6) = epszx(i)       
        ELSEIF(if2==-1)THEN      
          ep(i,4) = -epsxy(i)      
          ep(i,5) = -epsyz(i)      
          ep(i,6) = -epszx(i)      
        ENDIF                      
       ENDDO
C---------------------------
      DO j = 1, 6
        ic = 0
        DO i = 1, nel
          nfunc  = ipm(10,mat(i))
          IF (nfunc>=j) THEN
           aux = ipm(10+j,mat(i))
           IF (aux/=0) THEN
            ic = ic + 1
            index(ic) = i
            ipos1(ic) = nint(uvar(i,j))
            iad1(ic)  = npf(aux) / 2 + 1
            ilen1(ic) = npf(aux+1) / 2 - iad1(ic) - ipos1(ic)
            epc(ic) = ep(i,j)
           ENDIF
          ENDIF
        ENDDO
 
        IF (iresp==1) THEN
          CALL vinter2dp(tf,iad1,ipos1,ilen1,ic,epc,dydxv,yc)
        ELSE
          CALL vinter2(tf,iad1,ipos1,ilen1,ic,epc,dydxv,yc)
        ENDIF
C
        IF (j==1) THEN
#include "vectorize.inc"
          DO ii = 1, ic
            i = index(ii)
            uvar(i,j)=ipos1(ii)
            signxx(i)=sign(min(abs(signxx(i)),yc(ii)*fac1(i)),signxx(i))
          ENDDO
        ELSEIF (j==2) THEN
#include "vectorize.inc"
          DO ii = 1, ic
            i = index(ii)
            uvar(i,j)=ipos1(ii) 
            signyy(i)=sign(min(abs(signyy(i)),yc(ii)*fac2(i)),signyy(i))
          ENDDO
        ELSEIF (j==3) THEN
#include "vectorize.inc"
          DO ii = 1, ic
            i = index(ii)
            uvar(i,j)=ipos1(ii) 
            signzz(i)=sign(min(abs(signzz(i)),yc(ii)*fac3(i)),signzz(i))
          ENDDO
        ELSEIF (j==4) THEN
#include "vectorize.inc"
          DO ii = 1, ic
            i = index(ii)
            uvar(i,j)=ipos1(ii) 
            signxy(i)=sign(min(abs(signxy(i)),yc(ii)*fac4(i)),signxy(i))
          ENDDO
        ELSEIF (j==5) THEN
#include "vectorize.inc"
          DO ii = 1, ic
            i = index(ii)
            uvar(i,j)=ipos1(ii) 
            signyz(i)=sign(min(abs(signyz(i)),yc(ii)*fac5(i)),signyz(i))
          ENDDO
        ELSEIF (j==6) THEN
#include "vectorize.inc"
          DO ii = 1, ic
            i = index(ii)
            uvar(i,j)=ipos1(ii) 
            signzx(i)=sign(min(abs(signzx(i)),yc(ii)*fac6(i)),signzx(i))
          ENDDO
        ENDIF
C
      ENDDO
 
 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
 1100 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10,
     .          ' AT TIME :',g11.4)
      RETURN
      END