Fixed problems with Drucker-Prager yield criteria owing to a missing factor of 3 in pressure terms.

src/vrm.f90

@@ -100,7 +100,7 @@ select case (trim(plasticity_type))
                           flag_vrm_pb)
     endif
     friction_angle = phi*180.d0/pi
-    alpha = (2.d0*sin(phi))  /(sqrt(3.d0)*(3.d0-sin(phi)))
+    alpha = (6.d0*sin(phi))  /(sqrt(3.d0)*(3.d0-sin(phi)))
     k     = (6.d0*c*cos(phi))/(sqrt(3.d0)*(3.d0-sin(phi)))
     if (params%plastic_stress_min.lt.eps) then
       yield = max(c/compressibility,k-alpha*pressure)
@@ -127,7 +127,7 @@ select case (trim(plasticity_type))
                           flag_vrm_pb)
     endif
     friction_angle = phi*180.d0/pi
-    alpha = (2.d0*sin(phi))  /(sqrt(3.d0)*(3.d0+sin(phi)))
+    alpha = (6.d0*sin(phi))  /(sqrt(3.d0)*(3.d0+sin(phi)))
     k     = (6.d0*c*cos(phi))/(sqrt(3.d0)*(3.d0+sin(phi)))
     if (params%plastic_stress_min.lt.eps) then
       yield = max(c/compressibility,k-alpha*pressure)
@@ -154,7 +154,7 @@ select case (trim(plasticity_type))
                           flag_vrm_pb)
     endif
     friction_angle = phi*180.d0/pi
-    alpha = (2.d0*sin(phi))  /(sqrt(3.d0)*3.d0)
+    alpha = (6.d0*sin(phi))  /(sqrt(3.d0)*3.d0)
     k     = (6.d0*c*cos(phi))/(sqrt(3.d0)*3.d0)
     if (params%plastic_stress_min.lt.eps) then
 	  yield = max(c/compressibility,k-alpha*pressure)
@@ -181,7 +181,7 @@ select case (trim(plasticity_type))
                           flag_vrm_pb)
     endif
     friction_angle = phi*180.d0/pi
-    alpha = tan(phi)/(sqrt(9.d0+12.d0*(tan(phi))**2.d0))
+    alpha = (3.d0*tan(phi))/(sqrt(9.d0+12.d0*(tan(phi))**2.d0))
     k     = (3.d0*c)/(sqrt(9.d0+12.d0*(tan(phi))**2.d0))
     if (params%plastic_stress_min.lt.eps) then
       yield = max(c/compressibility,k-alpha*pressure)
@@ -296,19 +296,6 @@ select case (trim(plasticity_type))
     fail  = 2.d0*viscosity*e2d-yield
     if (fail.gt.0.d0) then
       viscosity=0.5d0*yield/e2d
-      !if (viscosity.le.params%viscositymin) then
-      !  write(*,*) 'VISCOSITY BELOW MINIMUM'
-      !  write(*,*) 'viscosity: ',viscosity
-      !  write(*,*) 'yield: ',yield
-      !  write(*,*) 'fail: ',fail
-      !  write(*,*) '-pressure: ',-pressure
-      !  write(*,*) 'phi: ',phi
-      !  write(*,*) 'sin_phi: ',sin_phi
-      !  write(*,*) 'c: ',c
-      !  write(*,*) 'cos_phi: ',cos_phi
-      !  write(*,*) 'zeta: ',zeta
-      !  write(*,*) 'e2d: ',e2d
-      !endif
       is_plastic=.true.
     else
       is_plastic=.false.