m2iter_imp.F File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	m2iter_imp (sig, epxe, aj2, g, ca, cb, cn, epd, sigmx, epmx, dpla1, ak, qh, sigy, fisokin, nel)

Function/Subroutine Documentation

m2iter_imp()

subroutine m2iter_imp	(	dimension(nel,6)	sig,
		dimension(nel)	epxe,
		dimension(nel)	aj2,
		dimension(nel)	g,
		dimension(nel)	ca,
			cb,
			cn,
		dimension(nel)	epd,
		dimension(nel)	sigmx,
			epmx,
		dimension(nel)	dpla1,
		dimension(nel)	ak,
		dimension(nel)	qh,
		dimension(nel)	sigy,
			fisokin,
		integer	nel )

Definition at line 28 of file m2iter_imp.F.

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