m2iter_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/.
!||====================================================================
!||    m2iter_imp   ../engine/source/materials/mat/mat002/m2iter_imp.F
!||--- called by ------------------------------------------------------
!||    m2law        ../engine/source/materials/mat/mat002/m2law.F
!||====================================================================
      SUBROUTINE m2iter_imp(
     1   SIG,     EPXE,    AJ2,     G,
     2   CA,      CB,      CN,      EPD,
     3   SIGMX,   EPMX,    DPLA1,   AK,
     4   QH,      SIGY,    FISOKIN, NEL)
c--------------------------------------------------------------
c
c    -------------------
c    -- parameters in --
c    -------------------
c
c    SIG   = Elastic Predictor (Deviatoric Stress Tensor)
c    EPXE  = previous time-step equivalent plastic strain  (on input)
c    AJ2   = von Mises stress at elastic predictor
c    G     = Kirhhoff's modulus (elastic shear modulus))
c
c    CA, CB, CN = coefficients of the hardening rule
c                 Sigma_Yield =   CA + CB * (equiv_pl_strn)**CN
c    EPD   = scalling fac due to velocity effect
c    SIGMX = max allowable von Mises stress   
c    EPMX  = max allowable equivalent plastic strain
c    FISOKIN= Mixed hardening ratio
c    --------------------
c    -- parameters out --      
c    --------------------
c
c    DPLA1 = equivalent plastic strain increment
c    EPXE  = updated  equivalent plastic strain (on output)
c    AK    = current yield stress
c    SIGY  = current yield stress (same as AK)
c    QH    = current hardening parameter
c
c--------------------------------------------------------------
c qz !! this subroutine will be optimized after
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER :: NEL
      my_real :: FISOKIN,EPMX,CB,CN
      my_real ,DIMENSION(NEL)   :: EPXE,EPD,G,AK,QH,AJ2,CA,SIGMX,DPLA1,SIGY
      my_real ,DIMENSION(NEL,6) :: SIG
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I, ITER, NITER
      MY_REAL :: XSI,DXSI,LHS,RHS,ALPHA_RADIAL,BETA,G3(MVSIZ),HTOT
C=======================================================================
c ---- max number of iterations ----
       niter = 10
C       
       DO 100 i=1,nel
C
       xsi = zero
C
C      -- PREVIOUS TIME-STEP EQUIVALENT PLASTIC STRAIN
       epxe(i) = max(zero,epxe(i))
       g3(i) =max(three*g(i),em15)
C
C      -- PREVIOUS TIME-STEP YIELD STRESS !!!!! if EPXE=0!!!!!
C       AK(I)=CA(I)(I)+CB*EPXE(I)**CN !..SIGMA0(EPS_P)
C
C      -- COMPARE IT WITH ELASTIC PREDICTOR
       IF(aj2(i)<ak(i)) GO TO 90 
C
       beta = one-fisokin
       DO 80 iter = 1,niter ! .. NR ITERATIONS ..
C
       IF(epxe(i)>zero) THEN !..TESTING FOR ZERO EQUIV.PL.STRN
         ak(i)=ca(i)+beta*cb*epxe(i)**cn !..SIGMA0(EPS_P)
C
C      --COMPUTE THE HARDENING PARAMETER H
C      --(JUST THE DERIVATIVE OF THE SIG-EPS_PL FUNCTION)
C
         IF(cn>one) THEN
           qh(i) = (cb*cn*epxe(i)**(cn - one))*epd(i)
         ELSEIF(cn==one) THEN !.. LINEAR HARDENING
           qh(i)= cb*epd(i)
         ELSE !.. POWER LESS THAN UNITY
           qh(i) = (cb*cn/epxe(i)**(one -cn))*epd(i)
         ENDIF
C
       ELSEIF(epxe(i)==zero) THEN
C
         ak(i)=ca(i)
C
         IF(cn>one )THEN
               qh(i) = zero
         ELSEIF(cn==one) THEN 
               qh(i) = cb*epd(i)
         ELSE
               qh(i) = ep10 *epd(i) !.. SHOULD BE +INFINITY
         ENDIF
C
       ELSE !--ERROR
CCC             WRITE(*,*)'M2LAW-NEGATIVE EQUIVALENT PLASTIC STRN',EPXE(I)
CCC             CALL ARRET(2)
C
       ENDIF !--EPXE(I)>ZERO ..END TESTING FOR ZERO EQUIV.PL.STRN
C
       htot = g3(i) + fisokin*qh(i)
       rhs = aj2(i) -htot * xsi - ak(i)
       lhs = g3(i) +  qh(i)
       dxsi = rhs/lhs
       xsi = xsi + dxsi
C
       epxe(i) = epxe(i) + dxsi 
       epxe(i) = max(zero,epxe(i))
          IF( abs(dxsi)<em10.AND.
     $        abs(rhs )<em10) GO TO 90
C
 80    CONTINUE
CCC      WRITE(*,*)'M2LAW--NON-CONVERGED ITERATION', ABS(DXSI),ABS(RHS)
 90    CONTINUE
C
       IF(xsi<zero) xsi = zero
       dpla1(i) = xsi !.. PLASTIC STRAIN MULTIPLIER "delta lambda"
CCC       IF(I==1)WRITE(*,*)'D_LAMBDA =',XSI
C
C     .. ENFORCING SOME CUT-OFFS, LEFT UNCHANGED
C  ----!!!cette partie devrait etre   interieure de boucle ---  tester
       ak(i)=ak(i)*epd(i)
       IF(sigmx(i)<ak(i))THEN
        ak(i)=sigmx(i)
        qh(i)=zero
       ENDIF
C       SIGY(I) = AK(I)
       IF(epxe(i)>epmx)THEN
         ak(i)=zero
         qh(i)=zero
       ENDIF
C
C     -- RADIAL RETURN ------
C     -- ALPHA_RADIAL := PARAMETER ALPHA OF THE RADIAL RETURN
C
       alpha_radial= min(one,ak(i)/ max(aj2(i),em15))
       sig(i,1)=alpha_radial*sig(i,1)
       sig(i,2)=alpha_radial*sig(i,2)
       sig(i,3)=alpha_radial*sig(i,3)
       sig(i,4)=alpha_radial*sig(i,4)
       sig(i,5)=alpha_radial*sig(i,5)
       sig(i,6)=alpha_radial*sig(i,6)
C
 100  CONTINUE                  !..LOOP ON ELEMENTS
C
      RETURN
      END