From beded7fa9369ba3fccc11436f9772f3178dc7cc1 Mon Sep 17 00:00:00 2001
From: Dave Whipp <dwhipp@dal.ca>
Date: Fri, 21 Jan 2011 18:02:58 +0000
Subject: [PATCH] Fixed typos in the DP plasticity options and added some debug
output
---
src/vrm.f90 | 30 +++++++++++++++++++++---------
1 file changed, 21 insertions(+), 9 deletions(-)
diff --git a/src/vrm.f90 b/src/vrm.f90
index 7414aefd..98011dc2 100644
--- a/src/vrm.f90
+++ b/src/vrm.f90
@@ -18,7 +18,7 @@
subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure, &
plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb, &
- straintot)
+ straintot,call_cnt)
use constants
use definitions
@@ -42,6 +42,7 @@ double precision :: J2d,J3d,viscosity,pressure,compressibility,straintot
character (len=8) plasticity_type
double precision plasticity_parameters(9)
logical :: is_plastic,flag_vrm_pb
+integer call_cnt
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
@@ -108,8 +109,19 @@ select case (trim(plasticity_type))
c = plasticity_parameters(2)
alpha = (2.d0*sin(phi)) /(sqrt(3.d0)*(3.d0-sin(phi)))
k = (6.d0*c*cos(phi))/(sqrt(3.d0)*(3.d0-sin(phi)))
- yield = max(c/compressibility,k+alpha*pressure)
- fail = 2.d0*viscosity*e2d-yield
+ yield = max(c/compressibility,k-alpha*pressure)
+ fail = 2.d0*viscosity*e2d-yield
+ !if (call_cnt.le.10) then
+ ! write(*,*) '******************************************'
+ ! write(*,*) 'phi: ',phi
+ ! write(*,*) 'c: ',c
+ ! write(*,*) 'alpha: ',alpha
+ ! write(*,*) 'k: ',k
+ ! write(*,*) 'pressure: ',pressure
+ ! write(*,*) 'yield: ',yield
+ ! write(*,*) 'fail: ',fail
+ ! write(*,*) '******************************************'
+ !endif
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
is_plastic=.true.
@@ -133,7 +145,7 @@ select case (trim(plasticity_type))
c = plasticity_parameters(2)
alpha = (2.d0*sin(phi)) /(sqrt(3.d0)*(3.d0+sin(phi)))
k = (6.d0*c*cos(phi))/(sqrt(3.d0)*(3.d0+sin(phi)))
- yield = max(c/compressibility,k+alpha*pressure)
+ yield = max(c/compressibility,k-alpha*pressure)
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
@@ -158,7 +170,7 @@ select case (trim(plasticity_type))
c = plasticity_parameters(2)
alpha = (2.d0*sin(phi)) /(sqrt(3.d0)*3.d0)
k = (6.d0*c*cos(phi))/(sqrt(3.d0)*3.d0)
- yield = max(c/compressibility,k+alpha*pressure)
+ yield = max(c/compressibility,k-alpha*pressure)
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
@@ -183,7 +195,7 @@ select case (trim(plasticity_type))
c = plasticity_parameters(2)
alpha = tan(phi)/(sqrt(9.d0+12.d0*(tan(phi))**2.d0))
k = (3.d0*c)/(sqrt(9.d0+12.d0*(tan(phi))**2.d0))
- yield = max(c/compressibility,k+alpha*pressure)
+ yield = max(c/compressibility,k-alpha*pressure)
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
@@ -195,7 +207,7 @@ select case (trim(plasticity_type))
case ('DPV')
alpha = plasticity_parameters(1)
k = plasticity_parameters(2)
- yield = max(k/compressibility,k+alpha*pressure)
+ yield = max(k/compressibility,k-alpha*pressure)
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
@@ -218,7 +230,7 @@ select case (trim(plasticity_type))
endif
endif
c = plasticity_parameters(2)
- yield = max(c/compressibility,c+tan(phi)*pressure)
+ yield = max(c/compressibility,c-tan(phi)*pressure)
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
@@ -243,7 +255,7 @@ select case (trim(plasticity_type))
c = plasticity_parameters(2)
sin_phi=sin(phi)
cos_phi=cos(phi)
- yield = max(plasticity_parameters(2)/compressibility,pressure*sin_phi+c*cos_phi)
+ yield = max(plasticity_parameters(2)/compressibility,-pressure*sin_phi+c*cos_phi)
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
--
GitLab