dprag51.F File Reference

#include "implicit_f.inc"
#include "com08_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	dprag51 (nel, sigd, vol, epseq, vfrac, deps, uparam, volume, eint, plas, pfrac, pold, pext, de)

Function/Subroutine Documentation

dprag51()

subroutine dprag51	(	integer, intent(in)	nel,
		dimension(6,nel), intent(inout)	sigd,
		dimension(nel), intent(in)	vol,
		dimension(nel), intent(inout)	epseq,
		dimension(nel), intent(in)	vfrac,
		dimension(6,nel), intent(in)	deps,
		dimension(*), intent(in)	uparam,
		dimension(nel), intent(in)	volume,
		dimension(nel), intent(inout)	eint,
		dimension(nel), intent(inout)	plas,
		intent(in)	pfrac,
		dimension(nel), intent(inout)	pold,
		intent(in)	pext,
		dimension(nel), intent(in)	de )

Definition at line 28 of file dprag51.F.

 
 
!common points with M10LAW()
        !SIGD   = deviatoric stress
        !T      = SIGD copy
        !G*DEPSXX = G * DT1*D1
 
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      "com08_c.inc"
C-----------------------------------------------
C   I N P U T   O U T P U T   A r g u m e n t s
C-----------------------------------------------
      INTEGER,INTENT(IN) :: NEL
      my_real, INTENT(IN) :: pfrac,deps(6,nel),vol(nel),uparam(*),volume(nel),pext,de(nel)
      my_real,INTENT(INOUT) :: sigd(6,nel),pold(nel),eint(nel),epseq(nel),plas(nel)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I,KK
      my_real :: gg,a0,a1,a2,amax,tmelt,thetl,amumx,sigmx
      my_real :: t0,fact
      my_real :: ptot
      my_real, intent(in) :: vfrac(nel)
      my_real :: young, nu
      my_real :: t1(nel), t2(nel), t3(nel), t4(nel), t5(nel), t6(nel)
      my_real :: r
      my_real :: bulk(nel), bulk2(nel), g(nel), g43(nel), g0(nel), g2(nel)
      my_real :: j2(nel),yield2(nel)
      my_real :: pstar,sigdo(6,nel)
      my_real :: einc, vm, vol_avg
      my_real :: dpla
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      a0           = uparam(16)
      a1           = uparam(17)
      a2           = uparam(18)
      amax         = uparam(19)
      young        = uparam(02)
      nu           = uparam(22)
      pstar        = uparam(23)
      tmelt        = uparam(08)
      thetl        = uparam(09)
      amumx        = uparam(20)
      sigmx        = uparam(11)
      t0           = uparam(13)
      gg           = uparam(01)
      bulk(1:nel)  = uparam(21)
      bulk2(1:nel) = uparam(21)
 
      DO i=1,nel
        g(i)    = gg !* DT1
        g2(i)   = two * g(i)
        g43(i)  = onep333*gg
 
        t1(i)=sigd(1,i)
        t2(i)=sigd(2,i)
        t3(i)=sigd(3,i)
        t4(i)=sigd(4,i)
        t5(i)=sigd(5,i)
        t6(i)=sigd(6,i)
 
        sigdo(1:6,i) = sigd(1:6,i)
        fact = vfrac(i)
        IF(vfrac(i) > two*em02)THEN
          t1(i) = t1(i) + g2(i)* (deps(1,i)-de(i))*fact
          t2(i) = t2(i) + g2(i)* (deps(2,i)-de(i))*fact
          t3(i) = t3(i) + g2(i)* (deps(3,i)-de(i))*fact
          t4(i) = t4(i) + g(i) * deps(4,i)*fact
          t5(i) = t5(i) + g(i) * deps(5,i)*fact
          t6(i) = t6(i) + g(i) * deps(6,i)*fact
        ENDIF
 
        j2(i)=half*(t1(i)**2+t2(i)**2+t3(i)**2)+t4(i)**2+t5(i)**2+t6(i)**2
        vm=sqrt(three*j2(i))
        ptot = pold(i)+pext
        g0(i) =a0+a1*ptot+a2*ptot*ptot
        g0(i)= min(amax,g0(i))
        g0(i)= max(zero,g0(i))
        IF(ptot <= pstar)g0(i)=zero
        IF(pold(i) < pfrac)g0(i)=zero
        yield2(i)=j2(i)-g0(i)
 
 
        if(vfrac(i) > two*em02) then
        
          r = zero
          IF(g0(i) > zero)THEN
            IF(yield2(i) > zero)THEN
              r = sqrt(three*g0(i))/(vm+ em14)
              dpla = (one - r)*vm  /max(three*g(i),em15)
              plas(i)  = plas(i) + dpla
              epseq(i) = epseq(i)+ dpla
            ELSE
              r = one-em02 ! 1-epsilon
            ENDIF
          ENDIF
                
          ! deviatoric stress
          sigd(1,i) = t1(i) * r
          sigd(2,i) = t2(i) * r
          sigd(3,i) = t3(i) * r
          sigd(4,i) = t4(i) * r
          sigd(5,i) = t5(i) * r
          sigd(6,i) = t6(i) * r
          !plastic work
          vol_avg = half*(vfrac(i)*volume(i)+vol(i))
          einc    = half*vol_avg*
     .                       ( (sigdo(1,i)+sigd(1,i)) * deps(1,i)
     .                       + (sigdo(2,i)+sigd(2,i)) * deps(2,i)
     .                       + (sigdo(3,i)+sigd(3,i)) * deps(3,i)
     .                       + (sigdo(4,i)+sigd(4,i)) * deps(4,i)
     .                       + (sigdo(5,i)+sigd(5,i)) * deps(5,i)
     .                       + (sigdo(6,i)+sigd(6,i)) * deps(6,i))
          eint(i) = eint(i) + einc
        elseif(vfrac(i) < em02)then
          plas(i) = zero
          sigd(1,i) = zero
          sigd(2,i) = zero
          sigd(3,i) = zero
          sigd(4,i) = zero
          sigd(5,i) = zero
          sigd(6,i) = zero
        else
         !smooth transition vfrac \in [0.01 0.02]
          r = (vfrac(i)-em02)*hundred
          plas(i) = r * plas(i)
          if(yield2(i) > zero)then
            r = r * sqrt(three*g0(i))/(vm+ em14)
          end if
          if(g0(i)  == zero) r = zero
          !--------------
          ! projection
          !--------------
          sigd(1,i) = sigd(1,i) * r
          sigd(2,i) = sigd(2,i) * r
          sigd(3,i) = sigd(3,i) * r
          sigd(4,i) = sigd(4,i) * r
          sigd(5,i) = sigd(5,i) * r
          sigd(6,i) = sigd(6,i) * r
        end if
 
      ENDDO !next I
 
 
      RETURN