Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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 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=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.
Dave Whipp
committed
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)
Dave Whipp
committed
! 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
Dave Whipp
committed
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)
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)
Dave Whipp
committed
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
Dave Whipp
committed
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)
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)
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
Dave Whipp
committed
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)
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)
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
Dave Whipp
committed
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)
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)
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)
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
Dave Whipp
committed
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)
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)
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')
Dave Whipp
committed
! Strain softening/hardening
Dave Whipp
committed
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)
friction_angle = phi*180.d0/pi
alpha = sin(phi)
k = c*cos(phi)
if (params%plastic_stress_min.lt.eps) then
Dave Whipp
committed
yield = max(c/compressibility,k-alpha*pressure)
yield = max(params%plastic_stress_min,k-alpha*pressure)
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)
Dave Whipp
committed
! Strain softening/hardening
phi = plasticity_parameters(1)*pi/180.d0
Dave Whipp
committed
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)
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
Dave Whipp
committed
yield = max(c/compressibility,(-pressure*sin_phi+c*cos_phi)/zeta)
yield = max(params%plastic_stress_min,(-pressure*sin_phi+c*cos_phi)/zeta)
fail = 2.d0*viscosity*e2d-yield
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
Dave Whipp
committed
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.
Dave Whipp
committed
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
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
Dave Whipp
committed
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
! 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
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------