Modfied to use an initial e2d value for the first plasticity calculation

src/vrm.f90

@@ -16,9 +16,12 @@
 !------------------------------------------------------------------------------|
 !------------------------------------------------------------------------------|
 
-subroutine vrm(J2d,J3d,compressibility,viscosity,pressure,plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb,straintot)
+subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure,             &
+               plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb,   &
+               straintot)
 
 use constants
+use definitions
 use invariants
 
 !------------------------------------------------------------------------------|
@@ -34,7 +37,7 @@ use invariants
 
 implicit none
 
-
+type (parameters) params
 double precision :: J2d,J3d,viscosity,pressure,compressibility,straintot
 character (len=8) plasticity_type
 double precision plasticity_parameters(9)
@@ -53,12 +56,24 @@ double precision :: strain_soft_yield
 
 flag_vrm_pb = .false.
 
+!write (*,*) 'vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv'
+!write (*,*) 'init_e2d: ',params%init_e2d
+!write (*,*) 'real e2d: ',sqrt(J2d)
+!write (*,*) 'viscosity: ',viscosity
+!write (*,*) 'is_plastic: ',is_plastic
+!write (*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
+!write (*,*) ''
+
+if (params%init_e2d) then
+  e2d=params%e2d0
+else
+  e2d=sqrt(J2d)
+endif
+
 select case (trim(plasticity_type))
 
    case ('vM')
-
       yield = plasticity_parameters(1)
-
       ! strain softening
       strain_soft_in=plasticity_parameters(3)
       if (strain_soft_in.gt.0.d0) then
@@ -69,10 +84,7 @@ select case (trim(plasticity_type))
             yield = max(yield, strain_soft_yield)
          end if
       end if
-
-    
-      e2d   = sqrt(J2d)
-      fail  = 2.d0*viscosity*e2d-yield 
+      fail = 2.d0*viscosity*e2d-yield
       if (fail.gt.0.d0) then
          viscosity=0.5d0*yield/e2d
          is_plastic=.true.
@@ -82,7 +94,7 @@ select case (trim(plasticity_type))
 
    case ('DPI','DPII')
       !yield = abs(plasticity_parameters(9)-plasticity_parameters(8)*pressure)
-! strain softening
+      ! strain softening
       phi=plasticity_parameters(1)
       strain_soft_in=plasticity_parameters(3)
       if (strain_soft_in.gt.0.d0) then
@@ -95,7 +107,6 @@ select case (trim(plasticity_type))
          endif
       endif
       yield = max(plasticity_parameters(2)/compressibility,plasticity_parameters(2)+tan(phi)*pressure)
-      e2d   = sqrt(J2d)
       fail  = 2.d0*viscosity*e2d-yield 
       if (fail.gt.0.d0) then
          viscosity=0.5d0*yield/e2d
@@ -105,7 +116,6 @@ select case (trim(plasticity_type))
       end if
 
    case ('DPIII') 
-      e2d   = sqrt(J2d)
       !yield = plasticity_parameters(2)+plasticity_parameters(1)*pressure
       yield = max(plasticity_parameters(2)/compressibility,plasticity_parameters(2)+plasticity_parameters(1)*pressure)
       fail  = 2.d0*viscosity*e2d-yield 
@@ -117,14 +127,17 @@ select case (trim(plasticity_type))
       end if
 
    case ('MC')
-      e2d       = sqrt(J2d)
-      theta     = lode_angle(J2d,J3d)
+      if (params%init_e2d) then
+        theta = 0.d0
+      else
+        theta = lode_angle(J2d,J3d)
+      endif
       sin_theta = sin(theta)
       cos_theta = cos(theta)
       sin_phi   = plasticity_parameters(8)
       cos_phi   = plasticity_parameters(9)
       c         = plasticity_parameters(2)
-! strain softening
+      ! strain softening
       strain_soft_in=plasticity_parameters(3)
       if (strain_soft_in.gt.0.d0) then
       strain_soft_out=plasticity_parameters(4)
@@ -139,23 +152,26 @@ select case (trim(plasticity_type))
         endif
       endif
       zeta      = cos_theta-sin_phi*sin_theta/sqrt3
-! modified by Jean
-!      fail      = -pressure*sin_phi+e2d*zeta-c*cos_phi
+      ! modified by Jean
+      !fail      = -pressure*sin_phi+e2d*zeta-c*cos_phi
       fail      = pressure*sin_phi+2.d0*viscosity*e2d*zeta-c*cos_phi
       if (fail.gt.0.d0) then
          viscosity=0.5d0/e2d*(c*cos_phi-pressure*sin_phi)/zeta
-!if (viscosity<0.d0) then
-!   viscosity=0.d0
-!   !print *,viscosity,c*cos_phi,pressure*sin_phi
-!end if
+         !if (viscosity<0.d0) then
+         !  viscosity=0.d0
+         !  !print *,viscosity,c*cos_phi,pressure*sin_phi
+         !end if
          is_plastic=.true.
       else
          is_plastic=.false.
       end if
 
    case ('Tresca')
-      e2d       = sqrt(J2d)
-      theta     = lode_angle(J2d,J3d)
+      if (params%init_e2d) then
+        theta = 0.d0
+      else
+        theta = lode_angle(J2d,J3d)
+      endif
       yield     = plasticity_parameters(1)
       cos_theta = cos(theta)
       fail  = 2.d0*viscosity*e2d*cos_theta-yield 
@@ -169,15 +185,22 @@ select case (trim(plasticity_type))
    case default
       print *,plasticity_type
       call stop_run('error in vrm$')
-       flag_vrm_pb = .true.
+      flag_vrm_pb = .true.
 end select
- 
+
+!write (*,*) 'vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv'
+!write (*,*) 'init_e2d: ',params%init_e2d
+!write (*,*) 'e2d out: ',sqrt(J2d)
+!write (*,*) 'viscosity out: ',viscosity
+!write (*,*) 'is_plastic out: ',is_plastic
+!write (*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
+!write (*,*) ''
+
 return
 end
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 
-
 subroutine compute_plastic_params (params,mat)
 use definitions
 use constants
@@ -234,7 +257,4 @@ do i=0,params%nmat
 end do
 
 return
-end
-
-
-
+end
\ No newline at end of file