sigeps120_connect_tab_vm.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/.
!||====================================================================
!||    sigeps120_connect_tab_vm   ../engine/source/materials/mat/mat120/sigeps120_connect_tab_vm.F
!||--- called by ------------------------------------------------------
!||    sigeps120_connect_main     ../engine/source/materials/mat/mat120/sigeps120_connect_main.F
!||--- calls      -----------------------------------------------------
!||    table_vinterp              ../engine/source/tools/curve/table_tools.F
!||--- uses       -----------------------------------------------------
!||    interface_table_mod        ../engine/share/modules/table_mod.F
!||    table_mod                  ../engine/share/modules/table_mod.F
!||====================================================================
      SUBROUTINE sigeps120_connect_tab_vm(
     1     NEL     ,NGL     ,TIME    ,TIMESTEP,UPARAM  ,OFF     ,
     2     EPSD    ,STIFM   ,THICK   ,JTHE    ,
     3     AREA    ,DEPSZZ  ,DEPSYZ  ,DEPSZX  ,NUPARAM ,
     4     SIGOZZ  ,SIGOYZ  ,SIGOZX  ,SIGNZZ  ,SIGNYZ  ,SIGNZX  ,
     5     PLA     ,JSMS    ,DMELS   ,UVAR    ,NUVAR   ,
     6     NUMTABL ,ITABLE  ,TABLE   ,NVARTMP ,VARTMP  ,TEMP    ,DMG)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE table_mod
      USE interface_table_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com04_c.inc"
#include      "units_c.inc"
#include      "comlock.inc"
#include      "sms_c.inc"
C----------------------------------------------------------------
C  D u m m y   A R G U M E N T S
C----------------------------------------------------------------
      INTEGER :: NEL,NUPARAM,NUVAR,JSMS,NUMTABL,NVARTMP,JTHE
      my_real :: TIME,TIMESTEP
      INTEGER ,DIMENSION(NEL)    ,INTENT(IN) :: NGL
      INTEGER ,DIMENSION(NUMTABL),INTENT(IN) :: ITABLE
      my_real ,DIMENSION(NUPARAM),INTENT(IN) :: UPARAM
      my_real ,DIMENSION(NEL)    ,INTENT(IN) :: THICK,TEMP,
     .         depszz,depsyz,depszx,
     .         sigozz,sigoyz,sigozx,area
c
      my_real ,DIMENSION(NEL) ,INTENT(OUT) :: epsd
      my_real ,DIMENSION(NEL) ,INTENT(OUT) :: 
     .                   signzz,signyz,signzx
      my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: off,pla,stifm,dmels,dmg
      my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
      INTEGER ,DIMENSION(NEL,NVARTMP) ,INTENT(INOUT) :: VARTMP
      TYPE(ttable) ,DIMENSION(NTABLE) ,INTENT(IN)    :: TABLE
C----------------------------------------------------------------
C  L O C A L  V A R I B L E S
C----------------------------------------------------------------
      INTEGER  :: I,II,ITER,NITER,NINDX,NINDF,ITRX,IDAM
      my_real :: xscale,yscale
      INTEGER :: IPOS1(NEL,1), IPOS2(NEL,2), IPOS3(NEL,3)
      INTEGER :: NDIM_YLD,FUNC_YLD
      my_real :: XVEC1(NEL,1), XVEC2(NEL,2), XVEC3(NEL,3)
      INTEGER  ,DIMENSION(NEL) :: INDX,INDF
c     material parameters
      my_real :: young,nu,g,g2,bulk,yld0,qq,beta,hh,c11,c12,
     .   a1f,a2f,as,d1c,d2c,d1f,d2f,d_trx,d_jc,dm,expn
c     Johnson & Cook rate-dependency
      my_real cc,epsdref,epsdmax,frate,fcut,alpha,alphai,epsdot
c     Local variables
      my_real facr,triax,rvoce,i1p,ip,jp,tp,fact,facd,dfac,dtinv,
     .   dlam,ddep,dfdlam,dyld_dpla,dpla_dlam,dphi_dlam,dtb,
     .   gamc,gamf,np_xx,np_yy,np_zz,np_xy,np_yz,np_zx,
     .   nf_xx,nf_yy,nf_zz,nf_xy,nf_yz,nf_zx,dpxx,dpyy,dpzz,dpxy,dpyz,dpzx
      my_real ,DIMENSION(NEL) :: i1,j2,yld,phi,dam,fdam,jcrate,stf,hardp,hardp_i,fvm,
     .   sxx,syy,szz,sxy,syz,szx,dpla,jcr_log,gamma,dgamm,pla_nl,dpdt_nl,yld_i
C---------------------------------------------
 
      young = uparam(1)
      nu    = uparam(2)
      g     = uparam(3)
      bulk  = uparam(4)
c     Parameters of yield function/hardening law and flow potential
      yld0  = uparam(5)
      qq    = uparam(6)
      beta  = uparam(7)
      hh    = uparam(8)
      a1f   = uparam(9) 
      a2f   = uparam(10)
c      A1H   = UPARAM(11)
c      A2H   = UPARAM(12)
      as    = uparam(13)
c     Johnson-Cook failure parameters
      d1c   = uparam(14)
      d2c   = uparam(15)
      d1f   = uparam(16)
      d2f   = uparam(17)
      d_trx = uparam(18)
      d_jc  = uparam(19)
      expn  = uparam(20) 
c     Johnson-Cook strain rate-dependency and distortional hardening
      cc      = uparam(21) 
      epsdref = uparam(22) 
      epsdmax = uparam(23) 
      fcut    = uparam(24)
c
      itrx    = nint(uparam(26)) 
      idam    = nint(uparam(27))
c
      xscale  = uparam(30)
      yscale  = uparam(31)
c
      func_yld = itable(1)
      ndim_yld = table(func_yld)%NDIM
c
      !YOUNG  = YOUNG/THICK
      !G      = G/THICK
      !BULK   = BULK/THICK
      g2     = g * two
      alpha  = 0.005
      alphai = one - alpha
      dtinv   = one / (max(timestep, em20))
C---------------------------------------------
      stf(1:nel)     = young *  area(1:nel)                                            
      stifm(1:nel)   = stifm(1:nel)  + stf(1:nel)*off(1:nel)                                              
C---------------------------------------------
      dam(1:nel) = dmg(1:nel)
      dpla(1:nel) = zero
      fdam(1:nel) = one - dam(1:nel)
      DO i = 1,nel
        IF (off(i) < 0.001) off(i) = zero
        IF (off(i) < one)    off(i) = off(i)*four_over_5
        IF (off(i) == one) THEN
          signzz(i) = sigozz(i)/fdam(i) + (depszz(i) /thick(i) )*young
          signyz(i) = sigoyz(i)/fdam(i) + (depsyz(i) /thick(i) )*g
          signzx(i) = sigozx(i)/fdam(i) + (depszx(i) /thick(i) )*g
          i1(i)     = signzz(i)
c
          sxx(i) = - i1(i) * third
          syy(i) = - i1(i) * third
          szz(i) = signzz(i) - i1(i) * third
          syz(i) = signyz(i)
          szx(i) = signzx(i)
          j2(i)   = (sxx(i)**2 + syy(i)**2 + szz(i)**2)*half
     .            +  syz(i)**2 + szx(i)**2
        END IF
      ENDDO
 
c     Johnson & Cook rate dependency  (total strain rate)
      jcrate(1:nel)  = one
      DO i=1,nel   
        IF (off(i) == one) THEN
          epsdot  = sqrt(depsyz(i)**2 + depszx(i)**2 + depszz(i)**2)
          epsdot  = epsdot * dtinv / thick(i)
          epsd(i) = min(epsdot ,epsdmax)
          epsd(i) = alpha *epsd(i) + alphai * uvar(i,3)
          uvar(i,3) = epsd(i)
          IF (epsd(i) > epsdref) THEN
                jcr_log(i) = log(epsd(i) / epsdref)
                jcrate(i)  = one + cc * jcr_log(i)
          END IF
        END IF
      ENDDO
c----------------------------------------------------      
c     Computation of the initial yield stress
c----------------------------------------------------      
      IF (ndim_yld == 1) THEN
        xvec1(1:nel,1) = pla(1:nel)
        ipos1(1:nel,1) = vartmp(1:nel,1)
c
        CALL table_vinterp(table(func_yld),nel,nel,ipos1,xvec1,yld,hardp)
c
        yld(1:nel)   = yld(1:nel)   * yscale * jcrate(1:nel)
        hardp(1:nel) = hardp(1:nel) * yscale * jcrate(1:nel)
        vartmp(1:nel,1) = ipos1(1:nel,1)
        
      ELSE IF (ndim_yld == 2) THEN
        xvec2(1:nel,1) = pla(1:nel)
        xvec2(1:nel,2) = epsd(1:nel)
        ipos2(1:nel,1) = vartmp(1:nel,1)
        ipos2(1:nel,2) = 1
c
        CALL table_vinterp(table(func_yld),nel,nel,ipos2,xvec2,yld,hardp)
c
        yld(1:nel)   = yld(1:nel)   * yscale
        hardp(1:nel) = hardp(1:nel) * yscale
        vartmp(1:nel,1) = ipos2(1:nel,1)
        
      ELSE
        xvec3(1:nel,1) = pla(1:nel)
        xvec3(1:nel,2) = epsd(1:nel)
        IF (jthe > 0) THEN
          xvec3(1:nel,3) = temp(1:nel)
        ELSE
          xvec3(1:nel,3) = zero
        END IF
        ipos3(1:nel,1) = vartmp(1:nel,1)
        ipos3(1:nel,2) = 1
        ipos3(1:nel,3) = 1
c
        CALL table_vinterp(table(func_yld),nel,nel,ipos3,xvec3,yld,hardp)
c
        yld(1:nel)   = yld(1:nel)   * yscale
        hardp(1:nel) = hardp(1:nel) * yscale
        vartmp(1:nel,1) = ipos3(1:nel,1)
        
      END IF
c
      DO i = 1,nel
        fvm(i) = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f )
        phi(i) = j2(i) + third*a2f * fvm(i)**2 
     .         - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
      ENDDO
c------------------------------------------------------------------
c     Check plasticity
c------------------------------------------------------------------
      nindx = 0
      nindf = 0
      DO i = 1,nel
        IF (phi(i) > zero .and. off(i) == one) THEN
          nindx = nindx + 1
          indx(nindx) = i
        END IF
      ENDDO
c
      IF (nindx > 0) THEN
      
        niter = 4      
        dpla(1:nel) = zero
      
        DO iter = 1,niter
c
          DO ii = 1,nindx
            i = indx(ii)
            
            ! normal to non associated plastic flow surface = d_psi/d_sig
            i1p   = max(zero, i1(i))
            jp    = third * as * i1p
            tp    = two  * jp
c
            np_xx = sxx(i) 
            np_yy = syy(i) 
            np_zz = szz(i) + tp 
            np_yz = syz(i)
            np_zx = szx(i)
c
            ! normal to plastic yield surface = d_phi/d_sig
            ip = two_third * a2f * fvm(i)
            nf_xx = sxx(i)
            nf_yy = syy(i) 
            nf_zz = szz(i) + ip
            nf_yz = syz(i)
            nf_zx = szx(i)
c         
            ! derivative of yld criterion : dphi/dlam = d_phi/d_sig * dsig/dlam
            dfdlam = - nf_zz * (young * np_zz)
     .             -(  nf_yz*np_yz + nf_zx*np_zx)*two*g2
c
            dpla_dlam = two*sqrt((j2(i) + tp*as*i1(i)))
            dphi_dlam = dfdlam - two*yld(i)*hardp(i)*dpla_dlam
c----------------------------------------------------------------
            dlam = -phi(i) / dphi_dlam
c----------------------------------------------------------------
            ! Plastic strains tensor update
                         
            dpzz = dlam*np_zz                                                   
            dpyz = dlam*np_yz                          
            dpzx = dlam*np_zx                 
            ! Elasto-plastic stresses update   
 
            signzz(i) = signzz(i) - young*dpzz     
            signyz(i) = signyz(i) - g2*dpyz  
            signzx(i) = signzx(i) - g2*dpzx
c
            ! Plastic strain increments update
            ddep   = dlam * dpla_dlam
            dpla(i)= max(zero, dpla(i) + ddep)
            pla(i) = pla(i) + ddep
          ENDDO  ! II = 1,NINDX
c-----------------------------------------        
c         Update yield from tabulated data   
c-----------------------------------------        
c
          IF (ndim_yld == 1) THEN
            DO ii = 1, nindx 
              i = indx(ii)
              xvec1(ii,1) = pla(i)
              ipos1(ii,1) = vartmp(i,1)
            ENDDO
c      
            CALL table_vinterp(table(func_yld),nel,nel,ipos1,xvec1,yld_i,hardp_i)
c      
            DO ii = 1, nindx 
              i = indx(ii)
              yld(i)   = yld_i(ii)*yscale
              hardp(i) = hardp_i(ii)*yscale
              vartmp(i,1) = ipos1(ii,1)
            ENDDO
            
          ELSE IF (ndim_yld == 2) THEN
            DO ii = 1, nindx 
              i = indx(ii)
              xvec2(ii,1) = pla(i)
              xvec2(ii,2) = epsd(i)
              ipos2(ii,1) = vartmp(i,1)
            ENDDO
c      
            CALL table_vinterp(table(func_yld),nel,nindx,ipos2,xvec2,yld_i,hardp_i)
c      
            DO ii = 1, nindx 
              i = indx(ii)
              yld(i)   = yld_i(ii)*yscale
              hardp(i) = hardp_i(ii)*yscale
              vartmp(i,1) = ipos2(ii,1)
            ENDDO
            
          ELSE
            DO ii = 1, nindx 
              i = indx(ii)
              xvec3(ii,1) = pla(i)
              xvec3(ii,2) = epsd(i)
              xvec3(ii,3) = temp(i)
              ipos3(ii,1) = vartmp(i,1)
            ENDDO
c      
            CALL table_vinterp(table(func_yld),nel,nel,ipos3,xvec3,yld_i,hardp_i)
c      
            DO ii = 1, nindx 
              i = indx(ii)
              yld(i)   = yld_i(ii)*yscale
              hardp(i) = hardp_i(ii)*yscale
              vartmp(i,1) = ipos3(ii,1)
            ENDDO            
          END IF
c----------------------        
c           Update stress invariants and criterion function value    
c----------------------        
          DO ii = 1,nindx
            i = indx(ii)
            i1(i)  = signzz(i) 
            sxx(i) = - i1(i)*third
            syy(i) = - i1(i)*third
            szz(i) = signzz(i) - i1(i)*third
            syz(i) = signyz(i)
            szx(i) = signzx(i)
            j2(i)  =(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) * half
     .           +   syz(i)**2 + szx(i)**2
            fvm(i) = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f )
            phi(i) = j2(i) + third*a2f * fvm(i)**2 
     .             - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
            uvar(i,4) = phi(i) / yld(i)**2
          END DO 
c
        END DO  ! Newton iterations
c
      END IF   ! NINDX > 0                                  
c---------------------------------------------------------------------
c     Update damage
c---------------------------------------------------------------------
      IF (idam > 0) THEN
        IF (idam == 1) THEN
            dgamm(1:nel) = dpla(1:nel)
            gamma(1:nel) = pla(1:nel)
        ELSE
            dgamm(1:nel)  = dpla(1:nel) * fdam(1:nel)
            uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
            gamma(1:nel) = uvar(1:nel,1)
        END IF
c
        DO ii = 1,nindx
          i = indx(ii)
c
          triax = zero
          fact = one
          facr = one + d_jc * jcr_log(i)  ! total strain rate factor on damage
          IF ( j2(i)/= zero)triax = third * i1(i) / sqrt(three*j2(i))
          IF (itrx == 1 ) THEN
            fact  = exp(-d_trx*triax)
          ELSE
           IF (triax > zero )fact  = exp(-d_trx*triax)
          ENDIF
         gamc = (d1c + d2c * fact) * facr
          gamf = (d1f + d2f * fact) * facr
          IF (gamma(i) > gamc) THEN
            IF (dam(i) == zero) THEN
              nindf = nindf + 1
              indf(nindf) = i
            END IF
            IF (expn == one) THEN
              dam(i) = dam(i) + dgamm(i) / (gamf - gamc)
            ELSE
              dfac   = (gamma(i) - gamc) / (gamf - gamc)
              dam(i) = dam(i) + expn * dfac**(expn-one) * dgamm(i) / (gamf - gamc)
            END IF
            IF (dam(i) >= one) THEN
              dam(i) = one
              off(i) = four_over_5
            END IF
            dmg(i) = dam(i)
          END IF ! GAMMA > GAMC
        ENDDO
c       Update damaged stress
        DO i = 1, nel
          facd = one - dam(i)
          signzz(i) = signzz(i) * facd
          signyz(i) = signyz(i) * facd
          signzx(i) = signzx(i) * facd        
        END DO
      END IF
c-------------------------
      IF (nindf > 0) THEN
        DO ii=1,nindf
#include "lockon.inc"
          WRITE(iout, 1000) ngl(indf(ii))
          WRITE(istdo,1100) ngl(indf(ii)),time
#include "lockoff.inc"
        ENDDO
      END IF
c-----------------------------------------------------------------
 1000 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10)
 1100 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10,1x,' AT TIME :',g11.4) 
c-----------------------------------------------------------------
c-----------------------------------------------------      
c omega = sqrt(2k/2*dmels), dt=2/omega, 2*dmels=dt**2 * 2k / 4
      IF (idtmins==2 .AND. jsms/=0 ) THEN
        dtb = (dtmins/dtfacs)**2
        DO i=1,nel                                                 
          dmels(i)=max(dmels(i),half*dtb*stf(i)*off(i))
        ENDDO                                                        
      ENDIF
c-----------------------------------------------------------------
      RETURN
      END