Skip to content
Snippets Groups Projects
Commit beded7fa authored by Dave Whipp's avatar Dave Whipp
Browse files

Fixed typos in the DP plasticity options and added some debug output

parent ab96bc63
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 21 additions and 9 deletions
src/vrm.f90
+ 21
− 9
View file @ beded7fa
...@@ -18,7 +18,7 @@ ...@@ -18,7 +18,7 @@
subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure, & subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure, &
plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb, & plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb, &
straintot) straintot,call_cnt)
use constants use constants
use definitions use definitions
...@@ -42,6 +42,7 @@ double precision :: J2d,J3d,viscosity,pressure,compressibility,straintot ...@@ -42,6 +42,7 @@ double precision :: J2d,J3d,viscosity,pressure,compressibility,straintot
character (len=8) plasticity_type character (len=8) plasticity_type
double precision plasticity_parameters(9) double precision plasticity_parameters(9)
logical :: is_plastic,flag_vrm_pb logical :: is_plastic,flag_vrm_pb
integer call_cnt
!------------------------------------------------------------------------------| !------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables ))))))))))))) !(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
...@@ -108,8 +109,19 @@ select case (trim(plasticity_type)) ...@@ -108,8 +109,19 @@ select case (trim(plasticity_type))
c = plasticity_parameters(2) c = plasticity_parameters(2)
alpha = (2.d0*sin(phi)) /(sqrt(3.d0)*(3.d0-sin(phi))) 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))) 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 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 if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d viscosity=0.5d0*yield/e2d
is_plastic=.true. is_plastic=.true.
...@@ -133,7 +145,7 @@ select case (trim(plasticity_type)) ...@@ -133,7 +145,7 @@ select case (trim(plasticity_type))
c = plasticity_parameters(2) c = plasticity_parameters(2)
alpha = (2.d0*sin(phi)) /(sqrt(3.d0)*(3.d0+sin(phi))) 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))) 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 fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d viscosity=0.5d0*yield/e2d
...@@ -158,7 +170,7 @@ select case (trim(plasticity_type)) ...@@ -158,7 +170,7 @@ select case (trim(plasticity_type))
c = plasticity_parameters(2) c = plasticity_parameters(2)
alpha = (2.d0*sin(phi)) /(sqrt(3.d0)*3.d0) alpha = (2.d0*sin(phi)) /(sqrt(3.d0)*3.d0)
k = (6.d0*c*cos(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 fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d viscosity=0.5d0*yield/e2d
...@@ -183,7 +195,7 @@ select case (trim(plasticity_type)) ...@@ -183,7 +195,7 @@ select case (trim(plasticity_type))
c = plasticity_parameters(2) c = plasticity_parameters(2)
alpha = tan(phi)/(sqrt(9.d0+12.d0*(tan(phi))**2.d0)) 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)) 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 fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d viscosity=0.5d0*yield/e2d
...@@ -195,7 +207,7 @@ select case (trim(plasticity_type)) ...@@ -195,7 +207,7 @@ select case (trim(plasticity_type))
case ('DPV') case ('DPV')
alpha = plasticity_parameters(1) alpha = plasticity_parameters(1)
k = plasticity_parameters(2) k = plasticity_parameters(2)
yield = max(k/compressibility,k+alpha*pressure) yield = max(k/compressibility,k-alpha*pressure)
fail = 2.d0*viscosity*e2d-yield fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d viscosity=0.5d0*yield/e2d
...@@ -218,7 +230,7 @@ select case (trim(plasticity_type)) ...@@ -218,7 +230,7 @@ select case (trim(plasticity_type))
endif endif
endif endif
c = plasticity_parameters(2) 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 fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d viscosity=0.5d0*yield/e2d
...@@ -243,7 +255,7 @@ select case (trim(plasticity_type)) ...@@ -243,7 +255,7 @@ select case (trim(plasticity_type))
c = plasticity_parameters(2) c = plasticity_parameters(2)
sin_phi=sin(phi) sin_phi=sin(phi)
cos_phi=cos(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 fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d viscosity=0.5d0*yield/e2d
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment