sigeps102.F File Reference

#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	sigeps102 (nel, nuparam, nuvar, uparam, rho0, rho, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, soundsp, uvar, off, et, psh, pnew, dpdm, ssp, pla)

Function/Subroutine Documentation

sigeps102()

subroutine sigeps102	(	integer	nel,
		integer	nuparam,
		integer	nuvar,
			uparam,
		dimension(nel)	rho0,
		dimension(nel)	rho,
		dimension(nel)	depsxx,
		dimension(nel)	depsyy,
		dimension(nel)	depszz,
		dimension(nel)	depsxy,
		dimension(nel)	depsyz,
		dimension(nel)	depszx,
		dimension(nel)	sigoxx,
		dimension(nel)	sigoyy,
		dimension(nel)	sigozz,
		dimension(nel)	sigoxy,
		dimension(nel)	sigoyz,
		dimension(nel)	sigozx,
			signxx,
			signyy,
			signzz,
			signxy,
			signyz,
			signzx,
			soundsp,
			uvar,
			off,
			et,
		dimension(nel)	psh,
		dimension(nel)	pnew,
		dimension(nel)	dpdm,
		dimension(nel)	ssp,
		dimension(nel)	pla )

Definition at line 28 of file sigeps102.F.

C-----------------------------------------------
C   I M P L I C I T   T Y P E S
C-----------------------------------------------
#include "implicit_f.inc"
C----------------------------------------------------------------
C  I N P U T   A R G U M E N T S
C----------------------------------------------------------------
      INTEGER NEL,NVARF,NUPARAM,NUVAR
      my_real :: time,timestep   
      my_real :: uparam(nuparam)
      my_real ,DIMENSION(NEL) :: rho, rho0,
     .   depsxx, depsyy, depszz, depsxy, depsyz, depszx,
     .   sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx,
     .   ssp   , dpdm  , pnew, psh, pla
C----------------------------------------------------------------
C  O U T P U T   A R G U M E N T S
C----------------------------------------------------------------
      my_real
     .      signxx(nel), signyy(nel), signzz(nel),
     .      signxy(nel), signyz(nel), signzx(nel),
     .      soundsp(nel), et(nel)
C----------------------------------------------------------------
C  I N P U T  O U T P U T   A R G U M E N T S
C----------------------------------------------------------------
      my_real
     .      uvar(nel,nuvar), off(nel) ,mu(nel),mu2(nel)
C----------------------------------------------------------------
C  L O C A L  V A R I B L E S
C----------------------------------------------------------------
      my_real :: a0,a1,a2,amax
      my_real :: dav,pold(nel)
      my_real :: t1(nel),t2(nel),t3(nel),t4(nel),t5(nel),t6(nel)
      my_real :: ptot,g0(nel),ratio(nel),yield2(nel)
      my_real :: pstar,g,gg,scrt(nel),aj2(nel),dpla(nel)
      my_real :: i3(nel),cos3t(nel),sqrt_j2,theta,c,phi,k
      integer :: I,IFORM
C----------------------------------------------------------------
C  S o u r c e   L i n e s
C----------------------------------------------------------------
      c         = uparam(1)
      phi       = uparam(2)
      pstar     = uparam(3)
      a0        = uparam(4)
      a1        = uparam(5)
      a2        = uparam(6)
      amax      = uparam(7) 
      g         = uparam(8)
      gg        = two*g
      iform     = nint(uparam(9))
      !----------------------------------------------------------------!
      !  state init.                                                   !
      !----------------------------------------------------------------!        
      DO i=1,nel
        pold(i) = -(sigoxx(i)+sigoyy(i)+sigozz(i))*third 
        scrt(i) =  (depsxx(i)+depsyy(i)+depszz(i))*third 
        mu(i)   = rho(i)/rho0(i) - one
        mu2(i)  = mu(i) * max(zero,mu(i))
      ENDDO !next I    
      !----------------------------------------------------------------!
      !  TEMPORARY DEVIATORIC STRESS TENSOR : T(1:6)                   !
      !----------------------------------------------------------------!  
      DO i=1,nel
        t1(i)=sigoxx(i)+pold(i)+gg*(depsxx(i)-scrt(i))
        t2(i)=sigoyy(i)+pold(i)+gg*(depsyy(i)-scrt(i))
        t3(i)=sigozz(i)+pold(i)+gg*(depszz(i)-scrt(i))
        t4(i)=sigoxy(i)        + g*depsxy(i)
        t5(i)=sigoyz(i)        + g*depsyz(i)
        t6(i)=sigozx(i)        + g*depszx(i)
      ENDDO !next I  
      !----------------------------------------------------------------!
      !  SOUND SPEED                                                   !
      !----------------------------------------------------------------!      
      DO i=1,nel
        dpdm(i) = dpdm(i) + onep333*g
        ssp(i)  = sqrt(abs(dpdm(i))/rho0(i))
      ENDDO !next I      
      !----------------------------------------------------------------!
      !  YIELD SURFACE                                                 !
      !----------------------------------------------------------------!     
      DO i=1,nel
        aj2(i)= half*(t1(i)**2+t2(i)**2+t3(i)**2)+t4(i)**2+t5(i)**2+t6(i)**2
      ENDDO
      !----SUBCASE --- ORIGINAL MOHR COULOMB
      IF(iform==4)THEN
        k = one/sqrt(three)
        DO i=1,nel
          i3(i)    = t2(i)*t3(i)*t1(i)-t2(i)*t6(i)*t6(i)-t3(i)*t4(i)*t4(i)-t5(i)*t5(i)*t1(i)+2*t5(i)*t4(i)*t6(i)
          sqrt_j2  = sqrt(aj2(i))
          cos3t(i) = nine*i3(i)/two/sqrt(three)/sqrt_j2/sqrt_j2/sqrt_j2
          theta    = acos(max(zero,min(one,cos3t(i))))
          ptot     = pnew(i)+psh(i)
          g0(i)    = -ptot*sin(phi)+sqrt_j2*(cos(theta)-k*sin(theta)*sin(phi))-c*cos(phi)
          g0(i)    = max(zero,g0(i))
          yield2(i)= aj2(i)-g0(i)
        ENDDO
      !----SUBCASE --- FITTED DRUCKER PRAGER FROM MOHR COULOMB PARAMETERS (A0,A1,A2 CALCULATED DURING STARTER)
      ELSE
        DO i=1,nel
          ptot     = pnew(i)+psh(i)
          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
          yield2(i)=aj2(i)-g0(i)
        ENDDO !next I           
      ENDIF 
 
      !----------------------------------------------------------------!
      !  PROJECTION FACTOR ON YIELD SURFACE                            !
      !----------------------------------------------------------------!      
      DO  i=1,nel
        ratio(i)=zero
        IF(yield2(i)<=zero .AND. g0(i)>zero)THEN
          ratio(i)=one
        ELSE
          ratio(i)=sqrt(g0(i)/(aj2(i)+ em14))
        ENDIF
      ENDDO !next I 
      !----------------------------------------------------------------!
      !  UPDATE DEVIATORIC STRESS TENSOR IN SIG(:,:)                   !
      !----------------------------------------------------------------!      
      DO i=1,nel
        signxx(i)=ratio(i)*t1(i)*off(i) - pnew(i)
        signyy(i)=ratio(i)*t2(i)*off(i) - pnew(i)
        signzz(i)=ratio(i)*t3(i)*off(i) - pnew(i)
        signxy(i)=ratio(i)*t4(i)*off(i)
        signyz(i)=ratio(i)*t5(i)*off(i)
        signzx(i)=ratio(i)*t6(i)*off(i)
        dpla(i)  =(one -ratio(i))*sqrt(aj2(i)) / max(em20,three*g)
      ENDDO !next I     
      pla(1:nel) = pla(1:nel) + dpla(1:nel)
c-----------         
      RETURN