mat104_nodam_newton.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/.
!||====================================================================
!||    mat104_nodam_newton   ../engine/source/materials/mat/mat104/mat104_nodam_newton.F
!||--- called by ------------------------------------------------------
!||    sigeps104             ../engine/source/materials/mat/mat104/sigeps104.f
!||====================================================================
      SUBROUTINE mat104_nodam_newton(
     1     NEL     ,NGL     ,NUPARAM ,NUVAR   ,VOLUME  ,FHEAT   , 
     2     TIME    ,TIMESTEP,UPARAM  ,UVAR    ,JTHE    ,OFF     ,
     3     RHO0    ,RHO     ,PLA     ,DPLA    ,EPSD    ,SOUNDSP ,
     4     DEPSXX  ,DEPSYY  ,DEPSZZ  ,DEPSXY  ,DEPSYZ  ,DEPSZX  ,
     5     SIGOXX  ,SIGOYY  ,SIGOZZ  ,SIGOXY  ,SIGOYZ  ,SIGOZX  ,
     6     SIGNXX  ,SIGNYY  ,SIGNZZ  ,SIGNXY  ,SIGNYZ  ,SIGNZX  ,
     7     SIGY    ,ET      ,TEMP    ,SEQ     ,INLOC   ,DPLANL  ,
     8     JLAG    )
      !=======================================================================
      !      Implicit types
      !=======================================================================
#include      "implicit_f.inc"
      !=======================================================================
      !      Common
      !=======================================================================
      !=======================================================================
      !      Dummy arguments
      !=======================================================================
      INTEGER NEL,NUPARAM,NUVAR,JTHE,INLOC
      INTEGER ,INTENT(IN) :: JLAG
      INTEGER ,DIMENSION(NEL), INTENT(IN)    :: NGL
      my_real 
     .   TIME,TIMESTEP
      my_real,DIMENSION(NUPARAM), INTENT(IN) :: 
     .   UPARAM
      my_real ,DIMENSION(NEL), INTENT(IN)    :: VOLUME
      my_real ,DIMENSION(NEL), INTENT(INOUT) :: fheat
      my_real ,DIMENSION(NEL), INTENT(IN)    :: rho0,rho,dplanl,
     .   depsxx,depsyy,depszz,depsxy,depsyz,depszx,
     .   sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
c
      my_real ,DIMENSION(NEL), INTENT(OUT)   :: 
     .   soundsp,sigy,et,
     .   signxx,signyy,signzz,signxy,signyz,signzx
c
      my_real ,DIMENSION(NEL), INTENT(INOUT)       :: 
     .   pla,dpla,epsd,off,temp,seq
      my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) :: 
     .   uvar
      !=======================================================================
      !      Local Variables
      !=======================================================================
      INTEGER I,II,ITER,NITER,NINDX,INDEX(NEL)
c
      my_real 
     .   young,bulk,lam,g,g2,nu,cdr,kdr,hard,yld0,qvoce,bvoce,jcc,
     .   epsp0,mtemp,tini,tref,eta,cp,dpis,dpad,asrate,afiltr
      my_real 
     .   h,ldav,trdep,sigvm,omega,
     .   dtemp,fcosh,fsinh,dpdt,dlam,ddep
      my_real 
     .   dsdrdj2,dsdrdj3,
     .   dj3dsxx,dj3dsyy,dj3dszz,dj3dsxy,dj3dsyz,dj3dszx,
     .   dj2dsxx,dj2dsyy,dj2dszz,dj2dsxy,dj2dsyz,dj2dszx,
     .   dfdsxx,dfdsyy,dfdszz,dfdsxy,dfdsyz,dfdszx,
     .   normxx,normyy,normzz,normxy,normyz,normzx,
     .   sdpla,dphi_dtrsig,sig_dfdsig,dfdsig2,
     .   dphi_dsig,dphi_dyld,dphi_dfdr,df_dfs,dfs_dft,dphi_dft,
     .   dphi_dfs,dfn_dlam,dfsh_dlam,dfg_dlam,dft_dlam,dpdt_nl,
     .   dfn,dfsh,dfg,dft,dyld_dpla,dyld_dtemp,dtemp_dlam,normsig
c
      my_real, DIMENSION(NEL) ::
     .   dsigxx,dsigyy,dsigzz,dsigxy,dsigyz,dsigzx,trsig,
     .   sxx,syy,szz,sxy,syz,szx,sigm,j2,j3,sigdr,yld,weitemp,
     .   hardp,fhard,frate,ftherm,dtherm,fdr,depxx,depyy,depzz,depxy,depyz,depzx,
     .   phi,ft,fs,fg,fn,fsh,dpla_dlam,dphi_dlam,phi0,
     .   dpxx,dpyy,dpxy,dpyz,dpzx,dpzz,sigdr2,yld2i
      !=======================================================================
      !             DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW
      !=======================================================================
c      
      !=======================================================================
      ! - INITIALISATION OF COMPUTATION ON TIME STEP
      !=======================================================================
      ! Recovering model parameters
      ! Elastic parameters      
      young   = uparam(1)  ! Young modulus
      bulk    = uparam(2)  ! Bulk modulus
      g       = uparam(3)  ! Shear modulus
      g2      = uparam(4)  ! 2*Shear modulus
      lam     = uparam(5)  ! Lambda Hooke parameter
      nu      = uparam(6)  ! Poisson ration
      ! Plastic criterion and hardening parameters [Drucker, 1948]
      cdr     = uparam(12) ! Drucker coefficient
      kdr     = uparam(13) ! Drucker 1/K coefficient
      tini    = uparam(14) ! Initial temperature
      hard    = uparam(15) ! Linear hardening
      yld0    = uparam(16) ! Initial yield stress
      qvoce   = uparam(17) ! 1st Voce parameter
      bvoce   = uparam(18) ! 2nd Voce parameter
      ! Strain-rate dependence parameters
      jcc     = uparam(20) ! Johnson-Cook strain rate coefficient
      epsp0   = uparam(21) ! Johnson-Cook reference strain rate
      ! Thermal softening and self-heating parameters
      mtemp   = uparam(22) ! Thermal softening slope
      tref    = uparam(23) ! Reference temperature
      eta     = uparam(24) ! Taylor-Quinney coefficient
      cp      = uparam(25) ! Thermal mass capacity
      dpis    = uparam(26) ! Isothermal plastic strain rate
      dpad    = uparam(27) ! Adiabatic plastic strain rate
      ! Plastic strain-rate filtering parameters
      asrate  = uparam(28) ! Plastic strain rate filtering frequency
      afiltr  = asrate*timestep/(asrate*timestep + one)
c      
      ! Recovering internal variables
      DO i=1,nel
        IF (off(i) < em01) off(i) = zero
        IF (off(i) <  one) off(i) = off(i)*four_over_5
        ! Standard inputs
        dpla(i)  = zero    ! Initialization of the plastic strain increment
        et(i)    = one     ! Initialization of hourglass coefficient
        hardp(i) = zero    ! Initialization of hardening modulus
      ENDDO
c      
      ! Initialization of temperature and self-heating weight factor
      IF (time == zero .and. jthe == 0) THEN   
        temp(1:nel) = tini 
      ENDIF 
      IF (cp > zero .OR. jthe /= 0) THEN     
        IF (inloc == 0) THEN       
          DO i=1,nel
            IF (epsd(i) < dpis) THEN
              weitemp(i) = zero
            ELSEIF (epsd(i) > dpad) THEN
              weitemp(i) = one
            ELSE
              weitemp(i) = ((epsd(i)-dpis)**2 )
     .                * (three*dpad - two*epsd(i) - dpis)
     .                / ((dpad-dpis)**3)
            ENDIF
          ENDDO
        ENDIF
      ENDIF
c      
      ! Computation of the initial yield stress
      DO i = 1,nel
        ! a) - Hardening law
        fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
        ! b) - Correction factor for strain-rate dependence (Johnson-Cook)
        frate(i)  = one
        IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)   
        ! c) - Correction factor for thermal softening      
        ftherm(i) = one - mtemp * (temp(i) - tref)
        ! d) - Computation of the yield function and value check
        yld(i) = fhard(i)*frate(i)*ftherm(i)
        ! e) - Checking values
        yld(i) = max(em10, yld(i))
      ENDDO
c      
      !========================================================================
      ! - COMPUTATION OF TRIAL VALUES
      !========================================================================
      DO i=1,nel
        ! Computation of the trial stress tensor
        ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lam
        signxx(i) = sigoxx(i) + depsxx(i)*g2 + ldav
        signyy(i) = sigoyy(i) + depsyy(i)*g2 + ldav
        signzz(i) = sigozz(i) + depszz(i)*g2 + ldav
        signxy(i) = sigoxy(i) + depsxy(i)*g
        signyz(i) = sigoyz(i) + depsyz(i)*g
        signzx(i) = sigozx(i) + depszx(i)*g
        ! Computation of the trace of the trial stress tensor
        trsig(i)  = signxx(i) + signyy(i) + signzz(i)
        sigm(i)   =-trsig(i) * third
        ! Computation of the deviatoric trial stress tensor
        sxx(i) = signxx(i) + sigm(i)
        syy(i) = signyy(i) + sigm(i)
        szz(i) = signzz(i) + sigm(i)
        sxy(i) = signxy(i)
        syz(i) = signyz(i)
        szx(i) = signzx(i)
        ! Second deviatoric invariant
        j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
     .        +       sxy(i)**2 + syz(i)**2 + szx(i)**2
        ! Third deviatoric invariant
        j3(i) = sxx(i) * syy(i) * szz(i)  
     .        + sxy(i) * syz(i) * szx(i) * two
     .        - sxx(i) * syz(i)**2
     .        - syy(i) * szx(i)**2
     .        - szz(i) * sxy(i)**2 
        ! Drucker equivalent stress
        fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
        ! Checking equivalent stress values
        IF (fdr(i) > zero) THEN
          sigdr(i) = kdr * exp((one/six)*log(fdr(i)))  ! FDR(I)**(1/6)
        ELSE
          sigdr(i) = zero
        ENDIF
      ENDDO
c
      !========================================================================
      ! - COMPUTATION OF YIELD FONCTION
      !========================================================================
      phi(1:nel) = (sigdr(1:nel) / yld(1:nel))**2 - one
      ! Checking plastic behavior for all elements
      nindx = 0
      DO i=1,nel         
        IF ((phi(i) >= zero).AND.(off(i) == one)) THEN
          nindx=nindx+1
          index(nindx)=i
        ENDIF
      ENDDO
c      
      !====================================================================
      ! - PLASTIC CORRECTION WITH CUTTING PLANE METHOD (SEMI-IMPLICIT)
      !====================================================================         
c      
      ! Number of iterations
      niter = 3
c
      ! Loop over yielding elements
#include "vectorize.inc" 
      DO ii=1,nindx  
c        
        ! Number of the element with plastic behaviour                                          
        i = index(ii) 
c        
        ! Initialization of the iterative Newton procedure
        sigdr2(i) = sigdr(i)**2
        yld2i(i)  = one / yld(i)**2
        dpxx(i)   = zero
        dpyy(i)   = zero
        dpzz(i)   = zero
        dpxy(i)   = zero
        dpyz(i)   = zero
        dpzx(i)   = zero
      ENDDO
c
      ! Loop over the iterations     
      DO iter = 1, niter
#include "vectorize.inc" 
        ! Loop over yielding elements
        DO ii=1,nindx 
          i = index(ii)    
c        
          ! Note: in this part, the purpose is to compute for each iteration
          ! a plastic multiplier allowing to update internal variables to satisfy
          ! the consistency condition using the cutting plane semi-implicit 
          ! iterative procedure.
          ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
          ! -> PHI          : current value of yield function (known)
          ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
          !                   into account of internal variables kinetic : 
          !                   plasticity, temperature and damage (to compute)
c        
          ! 1 - Computation of DPHI_DSIG the normal to the yield surface
          !-------------------------------------------------------------
c        
          ! Derivative with respect to the equivalent stress
          yld2i(i)  = one/(yld(i)**2) 
          dphi_dsig = two*sigdr(i)*yld2i(i)
c          
          ! Computation of the Eulerian norm of the stress tensor
          normsig = sqrt(signxx(i)*signxx(i)
     .            + signyy(i)*signyy(i)
     .            + signzz(i)*signzz(i)  
     .            + two*signxy(i)*signxy(i)
     .            + two*signyz(i)*signyz(i)
     .            + two*signzx(i)*signzx(i))   
          normsig = max(normsig,one)
c        
          ! Derivative with respect to Fdr 
          fdr(i)     =  (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)       
          dphi_dfdr  =  dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))             
          dsdrdj2    =  dphi_dfdr*three*(j2(i)/(normsig**2))**2                             
          dsdrdj3    = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))       
          ! dJ3/dSig                                                        
          dj3dsxx =  two_third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
     .             -  third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
     .             -  third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)                 
          dj3dsyy = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
     .             + two_third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
     .             -  third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)             
          dj3dszz = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
     .              - third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
     .             + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
          dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i) + szx(i)*syz(i))/(normsig**2)                   
          dj3dsyz = two*(sxy(i)*szx(i) + syy(i)*syz(i) + syz(i)*szz(i))/(normsig**2)                     
          dj3dszx = two*(sxx(i)*szx(i) + sxy(i)*syz(i) + szx(i)*szz(i))/(normsig**2)                      
          ! dPhi/dSig
          normxx  = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx
          normyy  = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy
          normzz  = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz
          normxy  = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
          normyz  = two*dsdrdj2*syz(i)/normsig + dsdrdj3*dj3dsyz
          normzx  = two*dsdrdj2*szx(i)/normsig + dsdrdj3*dj3dszx  
c          
          ! 2 - Computation of DPHI_DLAMBDA
          !---------------------------------------------------------
c        
          !   a) Derivative with respect stress increments tensor DSIG
          !   --------------------------------------------------------
          dfdsig2 = normxx * normxx * g2
     .            + normyy * normyy * g2
     .            + normzz * normzz * g2
     .            + normxy * normxy * g
     .            + normyz * normyz * g
     .            + normzx * normzx * g   
c     
          !   b) Derivatives with respect to plastic strain P 
          !   ------------------------------------------------  
c          
          !     i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
          !     ----------------------------------------------------------------------------
          hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
          dyld_dpla  = hardp(i)*frate(i)*ftherm(i)
c          
          !     ii) Derivative of dPLA with respect to DLAM
          !     -------------------------------------------                    
          sig_dfdsig = signxx(i) * normxx
     .               + signyy(i) * normyy
     .               + signzz(i) * normzz
     .               + signxy(i) * normxy
     .               + signyz(i) * normyz
     .               + signzx(i) * normzx  
          dpla_dlam(i) = sig_dfdsig / yld(i)
c          
          !  c) Derivatives with respect to the temperature TEMP 
          !  ---------------------------------------------------
          IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
          !    i)  Derivative of the yield stress with respect to temperature dYLD / dTEMP
          !    ---------------------------------------------------------------------------  
            dyld_dtemp = -fhard(i)*frate(i)*mtemp
          !    ii) Derivative of the temperature TEMP with respect to DLAM
          !    -----------------------------------------------------------
            dtemp_dlam = weitemp(i)*eta/(rho0(i)*cp)*sig_dfdsig
          ELSE
            dyld_dtemp = zero
            dtemp_dlam = zero
          ENDIF
c          
          !  d) Derivative with respect to the yield stress
          !  ----------------------------------------------
          dphi_dyld = -two*sigdr2(i)/(yld(i)**3)
c          
          ! 3 - Computation of plastic multiplier and variables update
          !----------------------------------------------------------
c          
          ! Derivative of PHI with respect to DLAM
          dphi_dlam(i) = - dfdsig2 + (dphi_dyld*dyld_dpla*dpla_dlam(i))
          IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
            dphi_dlam(i) = dphi_dlam(i) + dphi_dyld*dyld_dtemp*dtemp_dlam
          ENDIF
          dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
c          
          ! Computation of the plastic multiplier
          dlam = -phi(i)/dphi_dlam(i)
c          
          ! Plastic strains tensor update
          dpxx(i) = dlam * normxx
          dpyy(i) = dlam * normyy
          dpzz(i) = dlam * normzz
          dpxy(i) = dlam * normxy
          dpyz(i) = dlam * normyz
          dpzx(i) = dlam * normzx  
c          
          ! Elasto-plastic stresses update
          signxx(i) = signxx(i) - dpxx(i)*g2
          signyy(i) = signyy(i) - dpyy(i)*g2
          signzz(i) = signzz(i) - dpzz(i)*g2
          signxy(i) = signxy(i) - dpxy(i)*g
          signyz(i) = signyz(i) - dpyz(i)*g
          signzx(i) = signzx(i) - dpzx(i)*g
          trsig(i)  = signxx(i) + signyy(i) + signzz(i)
          sigm(i)   =-trsig(i)  * third
          sxx(i)    = signxx(i) + sigm(i)
          syy(i)    = signyy(i) + sigm(i)
          szz(i)    = signzz(i) + sigm(i)
          sxy(i)    = signxy(i)
          syz(i)    = signyz(i)
          szx(i)    = signzx(i)          
c          
          ! Cumulated plastic strain and strain rate update           
          ddep    = (dlam/yld(i))*sig_dfdsig
          dpla(i) = max(zero, dpla(i) + ddep)
          pla(i)  = pla(i) + ddep
c
          ! Drucker equivalent stress update
          j2(i)    = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
     .             +         sxy(i)**2 + syz(i)**2 + szx(i)**2
          j3(i)    = sxx(i) * syy(i) * szz(i) + two   * sxy(i) * syz(i) * szx(i)
     .             - sxx(i) * syz(i)**2 - syy(i) * szx(i)**2 - szz(i) * sxy(i)**2 
          fdr(i)   = j2(i)**3 - cdr*(j3(i)**2)
          sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
c          
          ! Temperature update
          IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN   
            dtemp     = weitemp(i)*yld(i)*ddep*eta/(rho0(i)*cp)
            temp(i)   = temp(i) + dtemp
          ENDIF
          ftherm(i) = one - mtemp*(temp(i) - tref)
c
          ! Hardening law update
          fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
c
          ! Yield stress update
          yld(i) = fhard(i)*frate(i)*ftherm(i)
          yld(i) = max(yld(i), em10)
c        
          ! Yield function value update
          sigdr2(i) = sigdr(i)**2
          yld2i(i)  = one / yld(i)**2
          phi(i)    = sigdr2(i) * yld2i(i) - one
c
        ENDDO
        ! End of the loop over the iterations  
      ENDDO 
      ! End of the loop over the iterations 
      !===================================================================
      ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ITERATIVE METHOD
      !===================================================================
c      
      ! Storing new values
      DO i=1,nel  
        ! USR Outputs
        seq(i)     = sigdr(i)
        ! Plastic strain-rate (filtered)
        dpdt       = dpla(i) / max(em20,timestep)
        epsd(i)    = afiltr * dpdt + (one - afiltr) * epsd(i)
        ! Coefficient for hourglass
        IF (dpla(i) > zero) THEN 
          et(i)    = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
        ELSE
          et(i)    = one
        ENDIF
        ! Computation of the sound speed
        soundsp(i) = sqrt((bulk + four_over_3*g)/rho0(i))
        ! Storing the yield stress
        sigy(i)    = yld(i)
        ! Non-local temperature
        IF (inloc > 0 .AND. off(i) == one) THEN 
          IF (cp > zero .AND. jthe == 0) THEN     
            ! Computation of the weight factor
            dpdt_nl = max(dplanl(i),zero)/max(timestep,em20)
            IF (dpdt_nl < dpis) THEN
              weitemp(i) = zero
            ELSEIF (dpdt_nl > dpad) THEN
              weitemp(i) = one
            ELSE
              weitemp(i) = ((dpdt_nl-dpis)**2 )
     .                * (three*dpad - two*dpdt_nl - dpis)
     .                / ((dpad-dpis)**3)
            ENDIF
            dtemp   = weitemp(i)*dplanl(i)*yld(i)*eta/(rho0(i)*cp)
            temp(i) = temp(i) + dtemp 
          ENDIF
        ENDIF
      ENDDO
c      
C-----------------------------------------------
C     plastic work dissipated as heat - in case of /heat/mat
C-----------------------------------------------
      IF (jthe < 0 .AND. jlag /= 0) THEN
        DO i=1,nel
          fheat(i) = fheat(i) + eta*(one-ft(i))*weitemp(i)*yld(i)*dpla(i)*volume(i)
        ENDDO
      END IF
!-----------
      RETURN      
      END