sigeps102.F File Reference

#include "implicit_f.inc"

Go to the source code of this file.

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

Function/Subroutine Documentation

sigeps102()

subroutine sigeps102	(	integer, intent(in)	nel,
		integer, intent(in)	nuparam,
		dimension(nuparam), intent(in)	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,
		dimension(nel)	signxx,
		dimension(nel)	signyy,
		dimension(nel)	signzz,
		dimension(nel)	signxy,
		dimension(nel)	signyz,
		dimension(nel)	signzx,
		dimension(nel)	off,
		dimension(nel)	psh,
		dimension(nel)	pnew,
		dimension(nel)	dpdm,
		dimension(nel)	ssp,
		dimension(nel)	pla,
		intent(in)	pmin )

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,INTENT(IN)                   :: NEL,NUPARAM
      my_real,INTENT(IN)                   :: uparam(nuparam)
      my_real,INTENT(IN),DIMENSION(NEL)    :: rho, rho0
      my_real,INTENT(IN),DIMENSION(NEL)    :: depsxx, depsyy, depszz, depsxy, depsyz, depszx
      my_real,INTENT(IN),DIMENSION(NEL)    :: sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx
      my_real,INTENT(INOUT),DIMENSION(NEL) :: signxx, signyy, signzz, signxy, signyz, signzx
      my_real,INTENT(IN),DIMENSION(NEL)    :: pnew, psh
      my_real,INTENT(IN),DIMENSION(NEL)    :: off
      my_real,INTENT(IN)                   :: pmin
      my_real,INTENT(INOUT),DIMENSION(NEL) :: ssp, dpdm, pla
C----------------------------------------------------------------
C  L O C A L  V A R I B L E S
C----------------------------------------------------------------
      my_real :: a0,a1,a2,amax
      my_real :: 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
      my_real :: mu(nel),mu2(nel)
      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))
      !----------------------------------------------------------------!
      !  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    
      !----------------------------------------------------------------!
      !  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 (EOS + SHEAR)                                     !
      !----------------------------------------------------------------!      
      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))
          IF(ptot <= pmin) g0(i) = zero !Drucker-Prager cut-off
          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 <= pmin) g0(i) = zero !Drucker-Prager cut-off
          IF(ptot <= pstar)g0(i)=zero !tensile cut-off (PSTAR is the root of the yield function)
          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