m36iter_imp.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/.
!||====================================================================
!||    m36iter_imp   ../engine/source/materials/mat/mat036/m36iter_imp.F
!||--- called by ------------------------------------------------------
!||    sigeps36      ../engine/source/materials/mat/mat036/sigeps36.F
!||--- calls      -----------------------------------------------------
!||    finter2       ../engine/source/tools/curve/vinter.F
!||    imp_stop      ../engine/source/implicit/imp_solv.F
!||    vinter2       ../engine/source/tools/curve/vinter.F
!||    vinter2dp     ../engine/source/tools/curve/vinter.F
!||====================================================================
      SUBROUTINE m36iter_imp( SIGNXX,  SIGNYY,  SIGNZZ,
     $                        SIGNXY,  SIGNYZ,  SIGNZX,
     $                        PLA,        YLD,      G3,
     $                        YFAC,     DPLA1,       H,   
     $                        TF,        IAD1,  
     $                        ILEN1,     NEL,
     $                        FISOKIN,  VARTMP, PLA0, 
     $                        PLAM    ,Y1,  DYDX1,
     $                        IPFLAG, PFUN, IPFUN, IPOSP, 
     $                        NRATE,  NPF,  IADP,  ILENP,
     $                        PFAC, PSCALE1, DFDP, FAIL , 
     A                        NVARTMP ,
     B                        SIGBXX,SIGBYY,SIGBZZ,SIGBXY,SIGBYZ,SIGBZX)
c--------------------------------------------------------------
c
c      This subroutine computes plastic strain rate multiplier
c      "delta_lambda" for the von Mises material characterized by
c      piecewise-linear hardening.
c      In addition, it computes the deviatoric stresses
c      via radial return.
c
c    -------------------
c    -- parameters in --
c    -------------------
c
c    SIGNXX  = Elastic Predictor  (on input)
c              (XX-Component Deviatoric Stress Tensor)
c              (...)
c    SIGNZX  = Elastic Predictor  (on input)
c              (ZX-Component Deviatoric Stress Tensor)
c
c    PLA     = previous time-step equivalent plastic strain  (on input)
c    YLD     = previous time-step yield stress               (on input)
c    G3      = 3 x Kirchhoff's modulus (elastic shear modulus))
c    YFAC    = scaling factor of sigma-eps_plast relation
c
c    TF      = array containing discrete sigma-eps_plast relation
c    IAD1,  ILEN1 = pre-computed data required by interpolation
c                         routine "VINTER"
c    NEL     = size of the element group being processed
c    FISOKIN = indicator for isotropic/kinematic/moixed hardening
c    UVAR    = array containing various state variables 
c    IPFLAG  = flag for existence of pressure scaling (scalar)
c    PFUN    = flag for existence of pressure scaling (array)
c    PFAC    = scaling factor for pressure at time t (not needed)
c    PSCALE1 = pressure at time t + delta_t
c    IPOSP, ILENP, NRATE,  NPF,  IADP = arrays and parameters required
c              by pressure interpolating routine (for pressure scaling)
c
c
c    --------------------
c    -- parameters out --      
c    --------------------
c
c    DPLA1   = equivalent plastic strain increment
c    PLA     = updated equivalent plastic strain (on output)
c    H       = updated hardening parameter
c    SIGNXX  = updated XX-component of deviatoric stress tensor 
c              (...)
c    SIGNZX  = updated ZX-component of deviatoric stress tensor, 
c    UVAR    = array containing various state variables
c
c    -----------------
c    -- work arrays --      
c    -----------------
c
c    Y1     =  yield stress  
c    DYDX1  =  derivative in sig0-eps_plast curve
c
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      "parit_c.inc"
#include      "scr05_c.inc"
#include      "units_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
C
      INTEGER  NEL, NVARTMP
      INTEGER  IAD1(MVSIZ),  ILEN1(MVSIZ),
     .         IPFLAG, PFUN,  IPFUN, IPOSP(MVSIZ),
     .          IADP(MVSIZ), ILENP(MVSIZ), NRATE, NPF(*)
      INTEGER :: VARTMP(NEL,NVARTMP)
C
C     REAL
      my_real
     .   SIGNXX(*),  SIGNYY(*),  SIGNZZ(*),
     .   SIGNXY(*),  SIGNYZ(*),  SIGNZX(*), 
     .   PLA(*),         G3(*),     YLD(*),
     .   DPLA1(*),        H(*),
     .   tf(*),          y1(*),   dydx1(*),
     .   yfac(mvsiz,*),
     .   pla0(mvsiz), plam(mvsiz),
     .   fisokin,
     .   pfac(mvsiz),  pscale1(mvsiz), dfdp(mvsiz),
     .   fail(mvsiz)
      my_real
     .   sigbxx(nel),sigbyy(nel),sigbzz(nel),sigbxy(nel),sigbyz(nel),sigbzx(nel)
CC
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, II, ITER, NITER, I_BILIN, IPOS
C
C     REAL
      MY_REAL
     .   XSI, DXSI, LHS, RHS, ALPHA_RADIAL, VM, HEP, R, DPLA,
     .   YLD_CURR, YLD_PREV, YLD_M, H_M, P_SCAL, Y_SCAL,
     .   total_scal, yld_m_prev, h_m_prev,
     .   sxx, syy, szz, sxy, syz, szx
     
C
C-----------------------------------------------
C   External Functions
C-----------------------------------------------
C
      my_real
     .   finter2
      EXTERNAL finter2
C
C=======================================================================
C
C
c ---- hardening law flag: i_bilin = 1 -- bilinear hardening law
c ----                     i_bilin = 0 -- general nonlinear hardening
c
       i_bilin = 0
c
c ---- max number of iterations ----
       niter = 20   
c
c    -- scaling if pressure dependent yield function --
c    -- compute the scaling factor for pressure + delta-pressure
c
        DO i=1,nel
          pfac(i)  =  one
        ENDDO
c
        IF (ipflag == nel.AND.(iparit == 0)) THEN
          DO i=1,nel
            iposp(i) = vartmp(i,2)
            iadp(i)  = npf(ipfun) / 2 + 1
            ilenp(i) = npf(ipfun + 1) / 2 - iadp(i) - iposp(i)
          ENDDO
          IF (iresp == 1) THEN
            CALL vinter2dp(tf,iadp,iposp,ilenp,nel,pscale1,dfdp,pfac)
          ELSE
            CALL vinter2(tf,iadp,iposp,ilenp,nel,pscale1,dfdp,pfac)
          ENDIF
          pfac(i) = max(zero, pfac(i))
        ELSEIF (ipflag  >  0) THEN
          IF (pfun  >  0) THEN
            DO i=1,nel
              iposp(i) = vartmp(i,2)
              iadp(i)  = npf(ipfun) / 2 + 1
              ilenp(i) = npf(ipfun  + 1) / 2 - iadp(i) - iposp(i)
              pfac(i)  = finter2(tf   ,iadp(i),iposp(i),ilenp(i),
     .                         pscale1(i),dfdp(i))
              pfac(i) = max(zero, pfac(i))
            ENDDO
          ENDIF
        ENDIF ! .. IPFLAG ...
c    -- end of computation of pressure-scaling factor --
c
c
c    -- store equivalent plastic strain at the beginning of time-step
c    -- store "fractional" equivalent plastic strain
c
       DO i=1,nel
         pla0(i) = pla(i)
         plam(i) = (one - fisokin) * pla(i) 
       ENDDO
c
c      -- computes the yield stress --
c
       DO i=1,nel
C  2Andrzej--------!!!IPOS is not initialized-----       
          ipos = 0
          yld(i) = finter2( tf, iad1(i), ipos, ilen1(i),
     $                         plam(i),dydx1(i) )
C------------------same than explicit
          yld(i)  = max(yld(i),em20)
       ENDDO
c
       DO 100 i = 1, nel ! .. loop on a group of elements ..
c
c     ... temporarily stores elastic predictor
          sxx = signxx(i)
          syy = signyy(i)
          szz = signzz(i)
          sxy = signxy(i)
          syz = signyz(i)
          szx = signzx(i)
c
          IF(i_bilin/=1 )THEN !..  general, nonlinear hardening
c
            xsi = zero ! .. plastic strain rate multiplier ..
c
c          -- PREVIOUS TIME-STEP EQUIVALENT PLASTIC STRAIN
            pla(i) = max(zero, pla(i))
c
c          -- von Mises stress at elastic predictor --
            vm =three*(half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
     $                        +signxy(i)**2+signyz(i)**2+signzx(i)**2)
            vm =sqrt(vm)
c
c
c          -- COMPARE VM with PREVIOUS TIME-STEP YIELD STRESS --
c          -- YLD is already multiplied by PFAC (at pressure)
c
            IF( vm > yld(i) )then ! .. plastically active process
c
c
c          -- scaling for sigma-eps_plast function --
              y_scal = yfac(i,1)
c
c          -- scaling if pressure dependent yield function --
              p_scal = pfac(i) * fail(i) 
c
c          -- total scaling
              total_scal = y_scal * p_scal
c
c            ------------------------------------------------------
c            yield stress and hardening parameter at m*(eps)
c            m =: (1 - FISOKIN)
c            ------------------------------------------------------
c          -- compute yield stress and hardening parameter DYDX1
c             at equivalent plastic strain PLAM, via interpolation
c
              ipos = 0 !.. FINTER2 will search the whole table
c
c          -- yield stress at m*(eps)  --
              yld_m_prev = finter2( tf, iad1(i), ipos, ilen1(i),
     $                         plam(i),dydx1(i) )
c
c          -- hardening parameter --
              h_m_prev  = dydx1(i)
c
c          -- apply scaling --
              yld_m_prev =  yld_m_prev * total_scal 
              h_m_prev  = h_m_prev * total_scal
              yld_m_prev  = max(yld_m_prev,em20)
c          -- check if H is nonnegative and y0 positive
              if( .not.( h_m_prev >= zero ).or.
     $            .not.( yld_m_prev  > zero ) )then 
               write(istdo,1000)plam(i),h_m_prev,yld_m_prev
               write(iout,1000)plam(i),h_m_prev,yld_m_prev
               CALL imp_stop(-1)
              endif
c
c            -------------------
c            yield stress at eps
c            -------------------
c
              ipos = 0 !.. FINTER2 will search the whole table
c
c          -- yield stress at  eps --
c
              yld_prev  = finter2( tf, iad1(i), ipos, ilen1(i),
     $                             pla0(i),dydx1(i) )
c
c          -- apply scaling --
              yld_prev =  yld_prev * total_scal
              yld_prev =  max(yld_prev,em20)
c
c          -- check if H is nonnegative and y0 positive
              if( .not.( dydx1(i) >= zero ).or.
     $            .not.( yld_prev  > zero ) )then
               write(istdo,1000)pla0(i),dydx1(i),yld_prev
               write(iout,1000)pla0(i),dydx1(i),yld_prev
               CALL imp_stop(-1)
              endif
c
c
c
c          =================
c            NR iterations
c          =================
c
            DO 80 iter = 1,niter 
c
c            ---------------------------------------------------
c            yield stress and hardening parameter at eps + deps
c            ---------------------------------------------------
C          -- compute yield stress and hardening parameter DYDX1
c             at equivalent plastic strain PLA, via interpolation
c
              ipos = 0 !.. FINTER2 will search the whole table
c
c          -- yield stress at eps + deps --
c
              yld_curr = finter2( tf,iad1(i), ipos, ilen1(i),
     $                            pla(i),dydx1(i) )
c
c          -- hardening parameter --
              h(i)  = dydx1(i)
c
c          -- apply scaling --
              yld_curr =  yld_curr * total_scal
              h(i) = h(i) * total_scal
              yld_curr =  max(yld_curr,em20)
c
c          -- check if H is nonnegative and y0 positive
              if( .not.( h(i) >= zero ).or.
     $            .not.( yld_curr  > zero ) )then
               write(istdo,1000)pla(i),h(i),yld_curr
               write(iout,1000)pla(i),h(i),yld_curr
               CALL imp_stop(-1)
              endif
c
c
c            -----------------------------------------
c            RHS and LHS for Newton at the Gauss point
c            -----------------------------------------
c
              IF(fisokin == zero)then
                 rhs = vm - g3(i) * xsi - yld_curr
                 lhs = g3(i) +  h(i)
              ELSEIF(fisokin == one)then
                 rhs = vm - g3(i)*xsi + yld_prev - yld_curr - yld_m_prev
                 lhs = g3(i) +  h(i)
              ELSE
                 rhs = vm - g3(i)*xsi + yld_prev - yld_curr - yld_m_prev
                 lhs = g3(i) +  h(i)
              ENDIF    
c
c
c            --------------------------------------
c            Solution for Newton at the Gauss point
c            --------------------------------------
              dxsi = rhs/lhs
C          -- update plastic strain rate multiplier
              xsi = xsi + dxsi
C          -- update equivalent plastic strain --      
              pla(i) =  pla(i) + dxsi
              pla(i) =  max(zero,pla(i)) !..to prevent negative value   
c
c
c          -- convergence criterion should be stringent --           
              IF( abs(dxsi)<em10 ) GO TO 90
c
 80         CONTINUE ! -- ITER --
c          --  we need a warning here --
ccc            WRITE(*,*)
ccc     $     'M36ITER_IMP--NON-CONVERGED ITERATION', ABS(DXSI),ABS(RHS)
c
 90         CONTINUE
c
c       ... equivalent plastic strain increment
            dpla = max(zero,xsi) !..to prevent negative value 
c
c
c
c          --------------------------------------
c           update of stresses and back stresses
c          --------------------------------------
c
            IF(fisokin == zero)then
              alpha_radial =  one - ( g3(i)/max(vm,em20) )*dpla
c
              signxx(i) = signxx(i) * alpha_radial
              signyy(i) = signyy(i) * alpha_radial
              signzz(i) = signzz(i) * alpha_radial
              signxy(i) = signxy(i) * alpha_radial
              signyz(i) = signyz(i) * alpha_radial
              signzx(i) = signzx(i) * alpha_radial
c
            ELSEIF(fisokin > zero.and.fisokin<=one)then
c
c          -- update "fractional" equivalent plastic strain --  
              plam(i) = (one - fisokin)*pla(i) !. for mixed hardeneing
c
c
c            ------------------------------------------------------
c            yield stress and hardening parameter at m*(eps+deps)
c            m =: (1 - FISOKIN)
c            ------------------------------------------------------
c          -- compute yield stress and hardening parameter DYDX1
c             at equivalent plastic strain PLAM, via interpolation
c
              ipos = 0 !.. FINTER2 will search the whole table
c
c          -- yield stress at m*(eps+deps)  --
              yld_m = finter2( tf, iad1(i), ipos, ilen1(i),
     $                         plam(i),dydx1(i) )
c
c          -- hardening parameter --
              h_m  = dydx1(i)
c          -- apply scaling --
              yld_m =  yld_m * total_scal 
              h_m  = h_m * total_scal
              yld_m =  max(yld_m,em20)
c
c          -- check if H is nonnegative and y0 positive
              if( .not.( h_m >= zero ).or.
     $            .not.( yld_m  > zero ) )then 
               write(istdo,1000)plam(i),h_m,yld_m
               write(iout,1000)plam(i),h_m,yld_m
               CALL imp_stop(-1)
              endif
c
c
c
c           ..updates (shifted) stresses
c
c          -- ALPHA_RADIAL := PARAMETER ALPHA OF THE RADIAL RETURN
c
              alpha_radial =  one - ( g3(i)/max(vm,em20) )*dpla
     $        - (yld_curr - yld_prev -yld_m  + yld_m_prev )/max(vm,em20)
c
              signxx(i) = signxx(i) * alpha_radial
              signyy(i) = signyy(i) * alpha_radial
              signzz(i) = signzz(i) * alpha_radial
              signxy(i) = signxy(i) * alpha_radial
              signyz(i) = signyz(i) * alpha_radial
              signzx(i) = signzx(i) * alpha_radial
c
c           ..updates back stresses
              hep =  ( yld_curr - yld_prev
     $                            -yld_m  + yld_m_prev  )/max(vm,em20)
c
              sigbxx(i)  = sigbxx(i)  + sxx *hep
              sigbyy(i)  = sigbyy(i)  + syy *hep
              sigbzz(i)  = sigbzz(i)  + szz *hep
              sigbxy(i)  = sigbxy(i) + sxy *hep 
              sigbyz(i)  = sigbyz(i) + syz *hep
              sigbzx(i)  = sigbzx(i) + szx *hep  
C
c           --adds the back stresses  to shifted stresses
             signxx(i) = signxx(i) + sigbxx(i) 
             signyy(i) = signyy(i) + sigbyy(i) 
             signzz(i) = signzz(i) + sigbzz(i) 
             signxy(i) = signxy(i) + sigbxy(i)
             signyz(i) = signyz(i) + sigbyz(i) 
             signzx(i) = signzx(i) + sigbzx(i) 
c
            ENDIF ! ..FISOKIN > ZERO.and.FISOKIN<=ONE .. 
c
c
c        .. updates yield stress.
            IF(fisokin == zero)then
              yld(i) = yld_curr   
            ELSEIF(fisokin == one)then
cccc              YLD(I) =  YLD(I)
            ELSE
              yld(i) =  yld_m
            ENDIF
c
            dpla1(i)  = dpla !--PLASTIC STRAIN MULTIPLIER "delta lambda"
 
c  
c
c           --- end of plastically active process ---
c
            ELSE !.. elasticity
 
c
c        .. gets stresses from shifted stresses and back stresses
            IF(fisokin > zero)then 
             signxx(i) = signxx(i) + sigbxx(i) 
             signyy(i) = signyy(i) + sigbyy(i) 
             signzz(i) = signzz(i) + sigbzz(i) 
             signxy(i) = signxy(i) + sigbxy(i)
             signyz(i) = signyz(i) + sigbyz(i) 
             signzx(i) = signzx(i) + sigbzx(i) 
            ENDIF
 
            ENDIF ! .. VM > YLD(I)
c
c-------------------------------------------------------------------
        ELSE ! -- i_bilin == 1
c   --- bilinear plastic hardening (isotropic, kinematic and mixed)
c-------------------------------------------------------------------
c
c         .. elastic predictor is temporarily stored in sxx, ... , szx
c
c         .. von Mises stress**2  at the elastic predictor ..
            vm =three*(half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
     $               +signxy(i)**2+signyz(i)**2+signzx(i)**2)
c
            IF(vm > yld(i)*yld(i))THEN !..  plastically active process
              vm =sqrt(vm)
              r = yld(i)/ max(vm,em20)
c          .. plastic strain increment.
              dpla=(one - r)*vm/max(g3(i)+h(i),em20)
c
c         ... updates stresses ..
              alpha_radial = one - g3(i)*( (one - r)/( g3(i) + h(i) ) )
c
              signxx(i) = signxx(i) * alpha_radial
              signyy(i) = signyy(i) * alpha_radial
              signzz(i) = signzz(i) * alpha_radial
              signxy(i) = signxy(i) * alpha_radial
              signyz(i) = signyz(i) * alpha_radial
              signzx(i) = signzx(i) * alpha_radial
c
              IF(fisokin > zero)then
c
c             --adds the back stresses at the beginning of the increment
                signxx(i) = signxx(i) + sigbxx(i)
                signyy(i) = signyy(i) + sigbyy(i)
                signzz(i) = signzz(i) + sigbzz(i)
                signxy(i) = signxy(i) + sigbxy(i)
                signyz(i) = signyz(i) + sigbyz(i)
                signzx(i) = signzx(i) + sigbzx(i)
c
c             ..updates back stresses
                hep = fisokin*h(i)* dpla/vm 
c
                sigbxx(i) = sigbxx(i) + sxx *hep
                sigbyy(i) = sigbyy(i) + syy *hep
                sigbzz(i) = sigbzz(i) + szz *hep
                sigbxy(i) = sigbxy(i) + sxy *hep 
                sigbyz(i) = sigbyz(i) + syz *hep
                sigbzx(i) = sigbzx(i) + szx *hep  
c
              ENDIF ! .. FISOKIN > ZERO ..       
c
c          .. updated yield stress.
              yld(i) = max(yld(i)+(one - fisokin)*dpla*h(i),zero)
c
c          .. updates equivalent plastic strain
              pla(i)    = pla(i) + dpla   
c  
              dpla1(i)  = dpla      
c
            ELSE !.. elasticity
c
            IF(fisokin > zero)then
c          .. gets stresses from shifted stresses and back stresses
              signxx(i) = signxx(i) + sigbxx(i)
              signyy(i) = signyy(i) + sigbyy(i)
              signzz(i) = signzz(i) + sigbzz(i)
              signxy(i) = signxy(i) + sigbxy(i)
              signyz(i) = signyz(i) + sigbyz(i)
              signzx(i) = signzx(i) + sigbzx(i)
            ENDIF
c
            ENDIF ! .. VM > YLD(I)*YLD(I)
c
c
        ENDIF ! .. i_bilin ..    
c
 100  CONTINUE  !..LOOP ON A GROUP OF ELEMENTS
c
 1000 FORMAT(3x,'Hardening parameters in law36 is not positive'/,
     $                   2x,'Eps_plast =',e11.4,2x,'H =',e11.4,
     $                   2x,'Sy0 =',e11.4)
c-----------
      RETURN
      END