sigeps91gpinch.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/.
!||====================================================================
!||    sigeps91gpinch   ../engine/source/materials/mat/mat091/sigeps91gpinch.F
!||--- called by ------------------------------------------------------
!||    mulawglcpinch    ../engine/source/elements/shell/coqueba/mulawglcpinch.F
!||====================================================================
      SUBROUTINE sigeps91gpinch(
     1                     JFT     ,JLT     ,NUVAR   ,UPARAM  ,RHO0    ,
     2                     THK     ,THK0    ,NEL     ,SSP     ,AREA    ,
     3                     DEPSXX  ,DEPSYY  ,DEPSZZ  ,
     4                     DEPSXY  ,DEPSYZ  ,DEPSZX  ,
     5                     DEPBXX  ,DEPBYY  ,DEPBXY  ,
     6                     DEPPXZ  ,DEPPYZ  ,
     7                     SIGOXX  ,SIGOYY  ,SIGOZZ  ,
     8                     SIGOXY  ,SIGOYZ  ,SIGOZX  ,
     9                     MOMOXX  ,MOMOYY  ,MOMOXY  ,
     A                     MOMOPXZ ,MOMOPYZ ,
     B                     SIGNXX  ,SIGNYY  ,SIGNZZ  ,
     C                     SIGNXY  ,SIGNYZ  ,SIGNZX  ,
     D                     MOMNXX  ,MOMNYY  ,MOMNXY  ,
     E                     MOMNPXZ ,MOMNPYZ ,TIME    ,UVAR    ,DT_INV   ,
     F                     VISCMX  ,ALDT    ,VOL0    ,IPM     ,MAT      ,
     G                     PLA     ,DEGMB   ,DEGFX   ,
     H                     NGL     ,EZZAVG  ,AREAPINCH)
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   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "param_c.inc"
C-----------------------------------------------
C   I N P U T   A R G U M E N T S
C-----------------------------------------------
      INTEGER JFT,JLT,NEL,ISMSTR,NUVAR, IPM(NPROPMI,*),MAT(NEL),NGL(*)
C     REAL
      my_real
     .  UPARAM(*),UVAR(NEL,NUVAR)
      my_real
     .  TIME
      my_real
     .   DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
     .   DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
     .   DEPBXX(NEL),DEPBYY(NEL),DEPBXY(NEL),
     .   DEPPXZ(NEL),DEPPYZ(NEL),
     .   SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
     .   SIGOXY(NEL),SIGOYZ(NEL),SIGOZX(NEL),
     .   MOMOXX(NEL),MOMOYY(NEL),MOMOXY(NEL),
     .   MOMOPXZ(NEL),MOMOPYZ(NEL),
     .   AREA(NEL),THK(NEL),THK0(NEL),RHO0(NEL),DT_INV(NEL),ALDT(MVSIZ),
     .   VOL0(MVSIZ),PLA(NEL),DEGMB(MVSIZ),DEGFX(MVSIZ),EZZAVG(MVSIZ),
     .   areapinch(mvsiz)
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),
     .    momnxx(nel),momnyy(nel),momnxy(nel),
     .    momnpxz(nel),momnpyz(nel),
     .    ssp(nel),rho(nel)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,MX,J,IADBUF
C     REAL
      my_real
     .   PA1,PA2,PA3,PA4,PA5,
     .   E,NU,G,GS,SIGY,HM,SIGY0,
     .   MS,FS,D1,D2,D3,C1,C2,C3, 
     .   AAA,BBB,CCC,P,BULK
      my_real
     .   vol(mvsiz),spheps(mvsiz), 
     .   devdepsxx(mvsiz),devdepsyy(mvsiz),devdepszz(mvsiz),
     .   devdepsxy(mvsiz),devdepsyz(mvsiz),devdepsxz(mvsiz),
     .   b1(mvsiz),b2(mvsiz),b3(mvsiz),b4(mvsiz),b5(mvsiz),
     .   thk08(mvsiz),viscmx(mvsiz),thkx(mvsiz),mu(mvsiz),
     .   svonm(mvsiz),crit(mvsiz),sphsig(mvsiz),magdev(mvsiz),invmagdev(mvsiz),
     .   devsxx(mvsiz),devsyy(mvsiz),devszz(mvsiz),
     .   devsxy(mvsiz),devsyz(mvsiz),devszx(mvsiz),
     .   n1(mvsiz),n2(mvsiz),n3(mvsiz),n4(mvsiz),n5(mvsiz),n6(mvsiz),
     .   rr(mvsiz),unsyeq(mvsiz),hh(mvsiz),etse(mvsiz),
     .   degsh_loc(mvsiz),degmb_loc(mvsiz),degfx_loc(mvsiz),dwelm(mvsiz),
     .   dwelf(mvsiz),dwpla(mvsiz),dpla(mvsiz),
     .   sdevxx(mvsiz),sdevyy(mvsiz),sdevzz(mvsiz),
     .   signdevxx(mvsiz),signdevyy(mvsiz),signdevzz(mvsiz),pnew(mvsiz),ptrial(mvsiz),
     .   pold(mvsiz),sigodevxx(mvsiz),sigodevyy(mvsiz),sigodevzz(mvsiz),dd(mvsiz)
C-----------------------------------------------
C
C initialize
      IF (time == zero) THEN
        DO i=1,nel
          uvar(i,1) = areapinch(i)*thk(i) 
          uvar(i,2) = thk(i)
        ENDDO
      ENDIF
C retrieve material parameters from the buffer       
        mx = mat(1)
        iadbuf = ipm(7,mx)
        e = uparam(iadbuf)
        nu  = uparam(iadbuf+1)
        sigy0    = uparam(iadbuf+2)
        hm      = uparam(iadbuf+3)
        g    = half*e/(one+nu)
        gs   = 5.d0/6.d0*g
        bulk = e/(three*(one-two*nu))
C other derived material parameters
        pa1  = e*(one-nu)/(one+nu)/(one-two*nu)
        pa2  = e*nu/(one+nu)/(one-two*nu)
        pa3  = g
        pa4  = (one+nu)*(one-two*nu)/(one-nu**2)/(one-nu)*pa1
        pa5  = (one+nu)*(one-two*nu)/(one-nu**2)*pa2 
C C1/C2/C3 are PM(28:30) [ref:sigeps02g]
        c1 = one/e
        c2 = -nu*c1
        c3 = one/g
C
      DO i=1,nel
        thk08(i)= thk0(i)*one_over_12
        b1(i)   = pa1*thk08(i)
        b2(i)   = pa2*thk08(i)
        b3(i)   = pa3*thk08(i)
        b4(i)   = pa4*thk08(i)
        b5(i)   = pa5*thk08(i)
      ENDDO
C
      DO i=1,nel
C        THKX(I) = UVAR(I,2)*(1+HALF*DEPSZZ(I))/(1-HALF*DEPSZZ(I))
         thkx(i) = uvar(i,2)*(1+half*ezzavg(i))/(1-half*ezzavg(i))
 
        uvar(i,2) = thkx(i)
        vol(i) = areapinch(i)*thkx(i)
        rho(i) = uvar(i,1)*rho0(i)/vol(i)
        mu(i) = rho(i)/rho0(i)-one
        pnew(i) = bulk*mu(i)
C        print*, NGL(I),PNEW(I)
      ENDDO
 
        
      DO i=jft,jlt
C       energy of transverse shear + pinching
        degsh_loc(i) = sigoyz(i)*depsyz(i)+sigozx(i)*depszx(i)
C       energy of only in plane components, xx,yy,xy
        degmb_loc(i) = degmb(i) - degsh_loc(i) 
C       bending energy
        degfx_loc(i) = degfx(i) 
C       calculate old pressure and old derivatoric stress
        pold(i) = -third*(sigoxx(i)+sigoyy(i)+sigozz(i))
        sigodevxx(i) = sigoxx(i)+pold(i)
        sigodevyy(i) = sigoyy(i)+pold(i)
        sigodevzz(i) = sigozz(i)+pold(i)
C       calculate change in spherical strain = DD
        dd(i) = third*(depsxx(i)+depsyy(i)+depszz(i))
C       update trial deviatoric stress to calculate trial deviatoric stress
        signdevxx(i)=sigodevxx(i)+two*g*(depsxx(i)-dd(i))
        signdevyy(i)=sigodevyy(i)+two*g*(depsyy(i)-dd(i))
        signdevzz(i)=sigodevzz(i)+two*g*(depszz(i)-dd(i))
C       shear stresses
        signxy(i)=sigoxy(i)+g*depsxy(i)
        signyz(i)=sigoyz(i)+gs*(depsyz(i)+zero*deppyz(i))
        signzx(i)=sigozx(i)+gs*(depszx(i)+zero*deppxz(i))
C       moments
        momnxx(i)=momoxx(i)+b4(i)*depbxx(i)+b5(i)*depbyy(i)
        momnyy(i)=momoyy(i)+b5(i)*depbxx(i)+b4(i)*depbyy(i)
        momnxy(i)=momoxy(i)+b3(i)*depbxy(i)
        momnpxz(i)=momopxz(i)+b3(i)*deppxz(i)
        momnpyz(i)=momopyz(i)+b3(i)*deppyz(i)
C       Sound speed like solids
        ssp(i) = sqrt(pa1/rho0(i))
 
 
 
 
 
 
      ENDDO
C     calculate equivalent stress 2D
      DO i=jft,jlt
        ms = momnxx(i)+momnyy(i)
C     theory manual (global plasticity algorithm)
        unsyeq(i) = one/
     .    sqrt(max(sixteen*(ms*ms + three*(momnxy(i)*momnxy(i) - momnxx(i)*momnyy(i)))
     .         +three*half*(signdevxx(i)**2+signdevyy(i)**2+signdevzz(i)**2+two*signxy(i)**2),em20))
        sigy = sigy0 + hm*pla(i)
        rr(i)  = min(one,sigy*unsyeq(i))
      ENDDO
C
      DO i=jft,jlt
        IF (rr(i) < one) THEN
          hh(i) = zero
          etse(i) = zero
        ENDIF
      ENDDO
C     energies calculation
      DO i=jft,jlt
C
          signxx(i) = signdevxx(i)*rr(i)-pnew(i)
          signyy(i) = signdevyy(i)*rr(i)-pnew(i)
          signxy(i) = signxy(i)*rr(i)
          signzz(i) = signdevzz(i)*rr(i)-pnew(i)
C
          d1 = signxx(i)-sigoxx(i)
          d2 = signyy(i)-sigoyy(i)
          d3 = signzz(i)-sigozz(i)
C
          dwelm(i) = (signxx(i)+sigoxx(i))*(c1*d1+c2*d2+c2*d3)+
     .               (signyy(i)+sigoyy(i))*(c2*d1+c1*d2+c2*d3)+
     .               (signzz(i)+sigozz(i))*(c2*d1+c2*d2+c1*d3)+
     .               (signxy(i)+sigoxy(i))*(c3*(signxy(i)-sigoxy(i)))
          degmb_loc(i) = degmb_loc(i)+signxx(i)*depsxx(i)+signyy(i)*depsyy(i)
     .                               +signzz(i)*depszz(i)
     .                               +signxy(i)*depsxy(i)
!
          momnxx(i) = momnxx(i)*rr(i)
          momnyy(i) = momnyy(i)*rr(i)
          momnxy(i) = momnxy(i)*rr(i)
          d1 = momnxx(i)-momoxx(i)
          d2 = momnyy(i)-momoyy(i)
          dwelf(i) = twelve*(
     .              (momnxx(i)+momoxx(i))*(c1*d1+c2*d2)
     .             +(momnyy(i)+momoyy(i))*(c2*d1+c1*d2)
     .             +(momnxy(i)+momoxy(i))*(c3*(momnxy(i)-momoxy(i))) )
          degfx_loc(i) = degfx_loc(i)+ momnxx(i)*depbxx(i)+momnyy(i)*depbyy(i)
     .                                +momnxy(i)*depbxy(i)
 
      ENDDO
C     delta of plastic work  
      DO i=jft,jlt
          dwpla(i) = degmb_loc(i)+degfx_loc(i)*thk0(i)-dwelm(i)-dwelf(i)
      ENDDO
C     equivalent plastic strain 
      DO i=jft,jlt
          dpla(i) = max(zero,half*dwpla(i)/sigy)
          pla(i) = pla(i) + dpla(i)
      ENDDO  
      RETURN
      END