sigeps53.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/.
!||====================================================================
!||    sigeps53      ../engine/source/materials/mat/mat053/sigeps53.F
!||--- called by ------------------------------------------------------
!||    mulaw         ../engine/source/materials/mat_share/mulaw.F90
!||    mulaw8        ../engine/source/materials/mat_share/mulaw8.F90
!||--- calls      -----------------------------------------------------
!||    ancmsg        ../engine/source/output/message/message.F
!||    arret         ../engine/source/system/arret.F
!||    finter2       ../engine/source/tools/curve/vinter.F
!||    vinter        ../engine/source/tools/curve/vinter.F
!||--- uses       -----------------------------------------------------
!||    message_mod   ../engine/share/message_module/message_mod.F
!||====================================================================
       SUBROUTINE sigeps53(
     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    ,SEQ_OUTPUT)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
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"
       INTEGER NEL, NUPARAM, NUVAR,IPT,IPLA
       INTEGER 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),seq_output(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(*), MFUNC, KFUNC(MFUNC)
       my_real
     .   FINTER2, TF(*)
        EXTERNAL FINTER2
C
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(MVSIZ),J1,J2,JJ(MVSIZ),NFUNC,
     .        NRATE(MVSIZ),IPOS1(MVSIZ),IPOS2(MVSIZ),IAD1(MVSIZ),
     .        ILEN1(MVSIZ),IAD2(MVSIZ),ILEN2(MVSIZ),
     .        IFUNC(MVSIZ,100),NFUNCV(MVSIZ),PFUN(MVSIZ),
     .        IPOSP(MVSIZ),IADP(MVSIZ),ILENP(MVSIZ),NFUNCM,NRATEM,
     .        IPFLAG,IPARAM,NPAR, IFLAG(MVSIZ)
       my_real   
     .        RATE(MVSIZ,2),YFAC(MVSIZ,2),
     .        Y1(MVSIZ),Y2(MVSIZ),H(MVSIZ),DYDX1(MVSIZ),
     .        DYDX2(MVSIZ),PLA(MVSIZ),FAIL(MVSIZ),EPSR1(MVSIZ),
     .        P0(MVSIZ),PFAC(MVSIZ),
     .        dfdp(mvsiz),dpla,
     .        evl(mvsiz),
     .        e11(mvsiz), e22(mvsiz), g12(mvsiz),g23(mvsiz),
     .        f1(mvsiz), f2(mvsiz), f11(mvsiz), f22(mvsiz), f12(mvsiz), 
     .        f44(mvsiz), f55(mvsiz), f23(mvsiz),f(mvsiz), s1c(mvsiz),   
     .        s2c(mvsiz),s2t(mvsiz), s3c(mvsiz),s3t(mvsiz), s4t(mvsiz),
     .        s4c(mvsiz), s45c(mvsiz),s45t(mvsiz),scale,bb,aa,dd,cc,ss1,
     .        ss2,s1t(mvsiz), evol(mvsiz)
          
C-----------------------------------------------
C     USER VARIABLES INITIALIZATION
C-----------------------------------------------
C
       DO i=1,nel
         iadbufv(i) = ipm(7,mat(i))-1
         e11(i)   =  uparam(iadbufv(i)+1)
         e22(i)   =  uparam(iadbufv(i)+2)
         g12(i)   =  uparam(iadbufv(i)+3)
         g23(i)   =  uparam(iadbufv(i)+4)
         iflag(i) = nint(uparam(iadbufv(i)+5))
C
         nfunc  = ipm(10,mat(i))
         DO j=1,nfunc
           ifunc(i,j) = ipm(10+j,mat(i))
         ENDDO    
       ENDDO
C
       IF(time==zero)THEN
         DO i=1,nel
            uvar(i,1)=zero
          DO j=1,nfunc
           uvar(i,j+1)=zero
          ENDDO
         ENDDO
       ENDIF
C-----------------------------------------------
C
       DO i=1,nel
          pla(i) = uvar(i,1)
        pla(i) = pla(i) +  depsxx(i) + depsyy(i) + depszz(i)
         signxx(i)=sigoxx(i)+ e11(i) * depsxx(i)
         signyy(i)=sigoyy(i)+ e22(i) * depsyy(i)
         signzz(i)=sigozz(i)+ e22(i) * depszz(i)
         signxy(i)=sigoxy(i)+ g12(i) * depsxy(i)
         signyz(i)=sigoyz(i)+ g23(i) * depsyz(i)
        signzx(i)=sigozx(i)+ g12(i) * depszx(i)
C
      soundsp(i) = sqrt(max(e11(i),e22(i),g12(i),g23(i))/rho0(i))
         evol(i) =one - exp(pla(i))
        viscmax(i) = zero
        uvar(i,1)= pla(i)
       ENDDO
C-------------------
C     CRITERE
C-------------------
        DO i=1,nel
            jj(i) = 1
             evl(i)=zero
        ENDDO
C
        DO j=1, nfunc
        DO i=1,nel
C---first direction--
         ipos1(i) = nint(uvar(i,1+j))
         iad1(i)  = npf(ifunc(i,j)) / 2 + 1
         ilen1(i) = npf(ifunc(i,j)+1) / 2 - iad1(i) - ipos1(i)
C--- second direction
C        IPOS2(I) = NINT(UVAR(I,3))
C        IAD2(I)  = 0.5*NPF(IFUNC(I,2)) + 1
C        ILEN2(I) = 0.5*NPF(IFUNC(I,3)+1) - IAD2(I) - IPOS2(I)
C--- third direction
C       IPOS3(I) = NINT(UVAR(I,4))
C        IAD3(I)  = 0.5*NPF(IFUNC(I,3)) + 1
C        ILEN3(I) = 0.5*NPF(IFUNC(I,4)+1) - IAD3(I) - IPOS3(I) 
C....fourth direction
C       IPOS4(I) = NINT(UVAR(I,5))
C        IAD4(I)  = 0.5*NPF(IFUNC(I,4)) + 1
C        ILE4(I)  = 0.5*NPF(IFUNC(I,5)+1) - IAD4(I) - IPOS4(I)
C... 45 direction curve
C       IPOS5(I) = NINT(UVAR(I,6))
C        IAD5(I)  = 0.5*NPF(IFUNC(I,5)) + 1
C        ILE5(I)  = 0.5*NPF(IFUNC(I,6)+1) - IAD5(I) - IPOS5(I)
        ENDDO
C  extrapollation dans chaque direction
C---first direction--
        CALL vinter(tf,iad1,ipos1,ilen1,nel,evol,dydx1,y1)
       CALL vinter(tf,iad1,ipos1,ilen1,nel,evl,dydx1,y2)
C
        IF(j==1) THEN
          DO i=1, nel
          IF(iflag(i)/=1)THEN
           IF(evol(i)>zero)THEN
            s1c(i)= y1(i) 
            s1t(i)= y2(i)
            ELSE
            s1c(i)= y2(i) 
            s1t(i)= y1(i)
           ENDIF
           ELSE
          s1c(i)= y1(i) 
            s1t(i)= y1(i)
         ENDIF
          uvar(i,1+j)= ipos1(i)
          ENDDO
C second direction
         ELSEIF(j==2)THEN
C          CALL vINTER(TF,IAD2,IPOS2,ILEN2,NEL,PLA,DYDX1,Y1)
C           CALL VINTER(TF,IAD2,IPOS2,ILEN2,NEL,EVL,DYDX1,Y2)
C
        DO i=1, nel
          IF(iflag(i)/=1)THEN
          IF(evol(i)>zero)THEN
            s2c(i)= y1(i) 
           s2t(i)= y2(i)
            ELSE
            s2c(i)= y2(i) 
           s2t(i)= y1(i)
          ENDIF
         ELSE
         s2c(i)= y1(i) 
          s2t(i)= y1(i)
          ENDIF
        uvar(i,1+j)= ipos1(i)
        ENDDO
CC third direction
       ELSEIF(j==3)THEN
C          CALL VINTER(TF,IAD3,IPOS3,ILEN3,NEL,PLA,DYDX1,Y1)
C           CALL VINTER(TF,IAD3,IPOS3,ILEN3,NEL,EVL,DYDX1,Y2)
C 
       DO i=1, nel
          IF(iflag(i)/=1)THEN
          IF(evol(i)>zero)THEN
            s3c(i)= y1(i) 
            s3t(i)= y2(i)
          ELSE
           s3c(i)= y2(i) 
            s3t(i)= y1(i)
           ENDIF
         ELSE
           s3c(i)= y1(i) 
           s3t(i)= y1(i)
          ENDIF
         uvar(i,1+j)= ipos1(i)
         ENDDO
CCC fourth direction
        ELSEIF(j==4)THEN
C          CALL VINTER(TF,IAD4,IPOS4,ILEN4,NEL,PLA,DYDX1,Y1)
C           CALL VINTER(TF,IAD4,IPOS4,ILEN4,NEL,EVL,DYDX1,Y2)
C
      DO i=1, nel
        IF(iflag(i)/=1)THEN
         IF(evol(i)>zero)THEN
          s4c(i)= y1(i) 
           s4t(i)= y2(i)
             ELSE
            s4c(i)= y2(i) 
            s4t(i)= y1(i)
           ENDIF
          ELSE
           s4c(i)= y1(i) 
           s4t(i)= y1(i)
         ENDIF
           uvar(i,1+j)= ipos1(i)
         ENDDO      
CCCC 45 direction
        ELSEIF(j==5)THEN
C          CALL VINTER(TF,IAD5,IPOS5,ILEN5,NEL,PLA,DYDX1,Y1)
C           CALL VINTER(TF,IAD5,IPOS5,ILEN5,NEL,EVL,DYDX1,Y2)
C
          DO i=1, nel
            IF(iflag(i)/=1)THEN
              IF(evol(i)>zero)THEN
               s45c(i)= y1(i) 
               s45t(i)= y2(i)
              ELSE
               s45c(i)= y2(i) 
               s45t(i)= y1(i)
              ENDIF
            ELSE
              s45c(i)= y1(i) 
              s45t(i)= y1(i)
            ENDIF   
            uvar(i,1+j)= ipos1(i)   
          ENDDO
        ENDIF
      ENDDO
C--YEILD SURFACE COEFFICIENTS
        DO i=1,nel         
          f1(i) = -one/s1c(i) + one/s1t(i)
          f2(i) = -one/s2c(i) + one/s2t(i)
          f11(i)= one/(s1c(i)*s1t(i))
          f22(i)= one/(s2c(i)*s2t(i))
          f44(i)= one/(s3c(i)*s3t(i))
          f55(i)= one/(s4c(i)*s4t(i))
          f12(i)=(-half)*sqrt(f11(i)*f22(i))
          f23(i)=(-half)*f22(i)     
C... CORRECTION
           IF(ifunc(i,5)>0)THEN
           f12(i)= two/(s45c(i)*s45c(i)) - 
     .        (half)*(f11(i)+ f22(i) + f44(i)) +
     .         (f1(i) + f2(i))/s45c(i)     
           ENDIF
C... Check the yield condition according to Tsay Wu.
       f(i) = f1(i)*signxx(i) + f2(i)* signyy(i) + f2(i)*signzz(i)+
     . f11(i)*signxx(i)**2 + f22(i)*signyy(i)**2 + f22(i)*signzz(i)**2
     .+f44(i)*signxy(i)**2 + f55(i)*signyz(i)**2 + f44(i)*signzx(i)**2
     .+two*f12(i)*signxx(i)*signyy(i) + two*f23(i)*signyy(i)*signzz(i)
     .+two*f12(i)*signzz(i)*signxx(i)
C
!!        SEQ_OUTPUT(I) = F(I)
C             
        IF(f(i)<=one) THEN
          scale=one
        ELSE !IF(F(I)>ONE)THEN
          bb= f1(i)*signxx(i) + f2(i)*signyy(i)+ f2(i)*signzz(i)
          cc = -one
          aa =
     . f11(i)*signxx(i)**2 + f22(i)*signyy(i)**2 + f22(i)*signzz(i)**2
     .+f44(i)*signxy(i)**2 + f55(i)*signyz(i)**2 + f44(i)*signzx(i)**2
     .+two*f12(i)*signxx(i)*signyy(i) + two*f23(i)*signyy(i)*signzz(i)
     .+two*f12(i)*signzz(i)*signxx(i)
C
          dd= bb**2 - four*aa*cc
          IF(dd<zero) THEN
            CALL ancmsg(msgid=136,anmode=aninfo)
            CALL arret(2)
          ENDIF
C...
          ss1 = (-bb + sqrt(dd))/(two*aa)
          ss2 = (-bb - sqrt(dd))/(two*aa)
C
          IF(ss1<=zero) scale = ss2 
          IF(ss2<=zero) scale = ss1 
C
          IF(ss1>0.AND.ss2>0) THEN
            WRITE(*,*)' TWO POSITIVE ROOTS IN STRANDFOAM'
            IF(ss1<ss2) scale = ss1
            IF(ss2<ss1) scale = ss2
          ENDIF
        ENDIF
CCC STRESS SCALING
C radial in stress space
        signxx(i) = signxx(i) * scale   
        signyy(i) = signyy(i) * scale   
        signzz(i) = signzz(i) * scale   
        signxy(i) = signxy(i) * scale   
        signyz(i) = signyz(i) * scale   
        signzx(i) = signzx(i) * scale   
 
      ENDDO
C-----
      RETURN
      END