!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! VRM Feb. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure, &
friction_angle,plasticity_parameters,plasticity_type, &
plasticity_ss_type_coh,plasticity_ss_type_phi,is_plastic, &
flag_vrm_pb,straintot,yield)
use constants
use definitions
use invariants
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
!This routine computes the rescaled viscosity in the element in the case it is
!above the failure criterion
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
implicit none
type (parameters) params
double precision :: J2d,J3d,viscosity,pressure,compressibility,straintot
double precision :: friction_angle
character (len=8) plasticity_type
character (len=12) plasticity_ss_type
double precision plasticity_parameters(9)
logical :: is_plastic,flag_vrm_pb
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
double precision :: yield,fail,theta,c,phi,e2d,zeta,alpha,k,eps
double precision :: sin_theta,cos_theta,sin_phi,cos_phi
double precision :: strain_soft_in_c,strain_soft_in_phi
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
flag_vrm_pb = .false.
eps = 1.d-10
c = 0.d0
phi = 0.d0
if (params%init_e2d) then
e2d=params%e2d0
else
e2d=sqrt(J2d)
endif
select case (trim(plasticity_type))
case ('vM')
yield = plasticity_parameters(1)
! Set initial cohesion (yield) value for cohesion (yield) weakening
strain_soft_in_c = plasticity_parameters(6)
! Call strain softening routine if initial cohesion (yield) is positive
if (strain_soft_in_c.gt.0.d0) then
call strain_soften (yield,phi,plasticity_ss_type_coh, &
plasticity_ss_type_phi,plasticity_parameters, &
straintot,e2d,flag_vrm_pb)
endif
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
is_plastic=.true.
else
is_plastic=.false.
end if
case ('DPI')
! Strain softening/hardening
phi = plasticity_parameters(1)*pi/180.d0
strain_soft_in_phi = plasticity_parameters(3)
c = plasticity_parameters(2)
strain_soft_in_c = plasticity_parameters(6)
if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
plasticity_parameters,straintot,e2d,flag_vrm_pb)
endif
friction_angle = phi*180.d0/pi
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)))
if (params%plastic_stress_min.lt.eps) then
yield = max(c/compressibility,k-alpha*pressure)
else
yield = max(params%plastic_stress_min,k-alpha*pressure)
endif
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
is_plastic=.true.
else
is_plastic=.false.
end if
case ('DPII')
! Strain softening/hardening
phi = plasticity_parameters(1)*pi/180.d0
strain_soft_in_phi = plasticity_parameters(3)
c = plasticity_parameters(2)
strain_soft_in_c = plasticity_parameters(6)
if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
plasticity_parameters,straintot,e2d,flag_vrm_pb)
endif
friction_angle = phi*180.d0/pi
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)))
if (params%plastic_stress_min.lt.eps) then
yield = max(c/compressibility,k-alpha*pressure)
else
yield = max(params%plastic_stress_min,k-alpha*pressure)
endif
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
is_plastic=.true.
else
is_plastic=.false.
end if
case ('DPIII')
! Strain softening/hardening
phi = plasticity_parameters(1)*pi/180.d0
strain_soft_in_phi = plasticity_parameters(3)
c = plasticity_parameters(2)
strain_soft_in_c = plasticity_parameters(6)
if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
plasticity_parameters,straintot,e2d,flag_vrm_pb)
endif
friction_angle = phi*180.d0/pi
alpha = (2.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)
else
yield = max(params%plastic_stress_min,k-alpha*pressure)
endif
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
is_plastic=.true.
else
is_plastic=.false.
end if
case ('DPIV')
! Strain softening/hardening
phi = plasticity_parameters(1)*pi/180.d0
strain_soft_in_phi = plasticity_parameters(3)
c = plasticity_parameters(2)
strain_soft_in_c = plasticity_parameters(6)
if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
plasticity_parameters,straintot,e2d,flag_vrm_pb)
endif
friction_angle = phi*180.d0/pi
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))
if (params%plastic_stress_min.lt.eps) then
yield = max(c/compressibility,k-alpha*pressure)
else
yield = max(params%plastic_stress_min,k-alpha*pressure)
endif
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
is_plastic=.true.
else
is_plastic=.false.
end if
case ('DPV')
alpha = plasticity_parameters(1)
k = plasticity_parameters(2)
if (params%plastic_stress_min.lt.eps) then
yield = max(k/compressibility,k-alpha*pressure)
else
yield = max(params%plastic_stress_min,k-alpha*pressure)
endif
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
is_plastic=.true.
else
is_plastic=.false.
end if
case ('DPVI')
! Strain softening/hardening
phi = plasticity_parameters(1)*pi/180.d0
strain_soft_in_phi = plasticity_parameters(3)
c = plasticity_parameters(2)
strain_soft_in_c = plasticity_parameters(6)
if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
plasticity_parameters,straintot,e2d,flag_vrm_pb)
endif
friction_angle = phi*180.d0/pi
alpha = tan(phi)
k = c
if (params%plastic_stress_min.lt.eps) then
yield = max(c/compressibility,k-alpha*pressure)
else
yield = max(params%plastic_stress_min,k-alpha*pressure)
endif
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
is_plastic=.true.
else
is_plastic=.false.
end if
case ('DPVII')
! Strain softening/hardening
phi = plasticity_parameters(1)*pi/180.d0
strain_soft_in_phi = plasticity_parameters(3)
c = plasticity_parameters(2)
strain_soft_in_c = plasticity_parameters(6)
if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
plasticity_parameters,straintot,e2d,flag_vrm_pb)
endif
friction_angle = phi*180.d0/pi
alpha = sin(phi)
k = c*cos(phi)
if (params%plastic_stress_min.lt.eps) then
yield = max(c/compressibility,k-alpha*pressure)
else
yield = max(params%plastic_stress_min,k-alpha*pressure)
endif
fail = 2.d0*viscosity*e2d-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d
is_plastic=.true.
else
is_plastic=.false.
end if
case ('MC')
if (params%init_e2d) then
theta = 0.d0
else
theta = lode_angle(J2d,J3d)
endif
sin_theta = sin(theta)
cos_theta = cos(theta)
! Strain softening/hardening
phi = plasticity_parameters(1)*pi/180.d0
strain_soft_in_phi = plasticity_parameters(3)
c = plasticity_parameters(2)
strain_soft_in_c = plasticity_parameters(6)
if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
plasticity_parameters,straintot,e2d,flag_vrm_pb)
endif
friction_angle = phi*180.d0/pi
sin_phi = sin(phi)
cos_phi = cos(phi)
zeta = cos_theta-sin_phi*sin_theta/sqrt3
if (params%plastic_stress_min.lt.eps) then
yield = max(c/compressibility,(-pressure*sin_phi+c*cos_phi)/zeta)
else
yield = max(params%plastic_stress_min,(-pressure*sin_phi+c*cos_phi)/zeta)
endif
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.
end if
case ('Tresca')
if (params%init_e2d) then
theta = 0.d0
else
theta = lode_angle(J2d,J3d)
endif
yield = plasticity_parameters(1)
! Set initial cohesion (yield) value for cohesion (yield) weakening
strain_soft_in_c = plasticity_parameters(6)
! Call strain softening routine if initial cohesion (yield) is positive
if (strain_soft_in_c.gt.0.d0) then
call strain_soften (yield,phi,plasticity_ss_type_coh, &
plasticity_ss_type_phi,plasticity_parameters, &
straintot,e2d,flag_vrm_pb)
endif
cos_theta = cos(theta)
fail = 2.d0*viscosity*e2d*cos_theta-yield
if (fail.gt.0.d0) then
viscosity=0.5d0*yield/e2d/cos_theta
is_plastic=.true.
else
is_plastic=.false.
end if
case default
print *,plasticity_type
call stop_run('error in vrm$')
flag_vrm_pb = .true.
end select
return
end subroutine vrm
!-------------------------------------------------------------------------------
subroutine strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
plasticity_parameters,straintot,e2d,flag_vrm_pb)
use constants
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine calculates strain softened values of either the angle of internal
! friction or cohesion/yield strength as a function of total strain or strain
! rate
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
implicit none
character (len=12) plasticity_ss_type_coh,plasticity_ss_type_phi
double precision plasticity_parameters(9)
double precision :: c,phi,straintot,e2d
logical :: flag_vrm_pb
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
double precision :: strain_soft_phi_in,strain_soft_phi_out,strain_soft_phi
double precision :: strain_soft_c_in,strain_soft_c_out,strain_soft_phi
double precision :: fact
! Set strain softening onset/termination and min values
strain_soft_c_in = plasticity_parameters(6)
strain_soft_c_out = plasticity_parameters(7)
strain_soft_c = plasticity_parameters(8)
strain_soft_phi_in = plasticity_parameters(3)
strain_soft_phi_out = plasticity_parameters(4)
strain_soft_phi = plasticity_parameters(5)*pi/180.d0
! Apply cohesion softening if onset cohesion value is positive
if (strain_soft_c_in.gt.0.d0) then
! Set type of strain softening
ss_type_c: select case (trim(plasticity_ss_type_coh))
case ('tot_strain')
strain_soft_ref=straintot
case ('strain_rate')
strain_soft_ref=e2d
case default
write (*,*) 'plasticity_ss_type_coh: ',plasticity_ss_type_coh
call stop_run('error in vrm$')
flag_vrm_pb = .true.
end select ss_type_c
! Soften cohesion if beyond onset value
if (strain_soft_ref.gt.strain_soft_in_c) then
fact=(strain_soft_ref-strain_soft_in_c)/(strain_soft_out_c-strain_soft_in_c)
fact=min(fact,1.d0)
c=c+(strain_soft_c-c)*fact
endif
endif
! Apply internal angle softening if onset phi value is positive
if (strain_soft_phi_in.gt.0.d0) then
! Set type of strain softening
ss_type_phi: select case (trim(plasticity_ss_type_phi))
case ('tot_strain')
strain_soft_ref=straintot
case ('strain_rate')
strain_soft_ref=e2d
case default
write (*,*) 'plasticity_ss_type_phi: ',plasticity_ss_type_phi
call stop_run('error in vrm$')
flag_vrm_pb = .true.
end select ss_type_phi
! Soften phi if beyond onset value
if (strain_soft_ref.gt.strain_soft_in_phi) then
fact=(strain_soft_ref-strain_soft_in_phi)/(strain_soft_out_phi-strain_soft_in_phi)
fact=min(fact,1.d0)
phi=phi+(strain_soft_phi-phi)*fact
endif
endif
end subroutine strain_soften
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------