!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! VRM Feb. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure, &
plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb, &
straintot)
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
character (len=8) plasticity_type
double precision plasticity_parameters(9)
logical :: is_plastic,flag_vrm_pb
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
double precision :: yield,fail,theta,c,e2d,zeta
double precision :: sin_theta,cos_theta,sin_phi,cos_phi
double precision :: strain_soft_in,strain_soft_out,strain_soft_phi,phi,fact
double precision :: strain_soft_yield
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
flag_vrm_pb = .false.
!write (*,*) 'vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv'
!write (*,*) 'init_e2d: ',params%init_e2d
!write (*,*) 'real e2d: ',sqrt(J2d)
!write (*,*) 'viscosity: ',viscosity
!write (*,*) 'is_plastic: ',is_plastic
!write (*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
!write (*,*) ''
if (params%init_e2d) then
e2d=params%e2d0
else
e2d=sqrt(J2d)
endif
select case (trim(plasticity_type))
case ('vM')
yield = plasticity_parameters(1)
! strain softening
strain_soft_in=plasticity_parameters(3)
if (strain_soft_in.gt.0.d0) then
strain_soft_out=plasticity_parameters(4)
strain_soft_yield=plasticity_parameters(2)
if (straintot.gt.strain_soft_in) then
yield = (strain_soft_yield-yield)/(strain_soft_out-strain_soft_in)*(straintot-strain_soft_out) + strain_soft_yield
yield = max(yield, strain_soft_yield)
end if
end if
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','DPII')
!yield = abs(plasticity_parameters(9)-plasticity_parameters(8)*pressure)
! strain softening
phi=plasticity_parameters(1)
strain_soft_in=plasticity_parameters(3)
if (strain_soft_in.gt.0.d0) then
strain_soft_out=plasticity_parameters(4)
strain_soft_phi=plasticity_parameters(5)/180.d0*3.141592654d0
if (straintot.gt.strain_soft_in) then
fact=(straintot-strain_soft_in)/(strain_soft_out-strain_soft_in)
fact=min(fact,1.d0)
phi=phi+(strain_soft_phi-phi)*fact
endif
endif
yield = max(plasticity_parameters(2)/compressibility,plasticity_parameters(2)+tan(phi)*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')
!yield = plasticity_parameters(2)+plasticity_parameters(1)*pressure
yield = max(plasticity_parameters(2)/compressibility,plasticity_parameters(2)+plasticity_parameters(1)*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)
sin_phi = plasticity_parameters(8)
cos_phi = plasticity_parameters(9)
c = plasticity_parameters(2)
! strain softening
strain_soft_in=plasticity_parameters(3)
if (strain_soft_in.gt.0.d0) then
strain_soft_out=plasticity_parameters(4)
strain_soft_phi=plasticity_parameters(5)/180.d0*3.141592654d0
if (straintot.gt.strain_soft_in) then
phi=atan2(plasticity_parameters(8),plasticity_parameters(9))
fact=(straintot-strain_soft_in)/(strain_soft_out-strain_soft_in)
fact=min(fact,1.d0)
phi=phi+(strain_soft_phi-phi)*fact
sin_phi=sin(phi)
cos_phi=cos(phi)
endif
endif
zeta = cos_theta-sin_phi*sin_theta/sqrt3
! modified by Jean
!fail = -pressure*sin_phi+e2d*zeta-c*cos_phi
fail = pressure*sin_phi+2.d0*viscosity*e2d*zeta-c*cos_phi
if (fail.gt.0.d0) then
viscosity=0.5d0/e2d*(c*cos_phi-pressure*sin_phi)/zeta
!if (viscosity<0.d0) then
! viscosity=0.d0
! !print *,viscosity,c*cos_phi,pressure*sin_phi
!end if
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)
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
!write (*,*) 'vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv'
!write (*,*) 'init_e2d: ',params%init_e2d
!write (*,*) 'e2d out: ',sqrt(J2d)
!write (*,*) 'viscosity out: ',viscosity
!write (*,*) 'is_plastic out: ',is_plastic
!write (*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
!write (*,*) ''
return
end
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
subroutine compute_plastic_params (params,mat)
use definitions
use constants
implicit none
type (parameters) params
type (material) :: mat(0:params%nmat)
double precision c,k,alpha,phi
integer i,iproc,nproc,ierr
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
do i=0,params%nmat
select case (trim(mat(i)%plasticity_type))
case ('No')
case ('vM')
case ('DPI')
phi = mat(i)%plasticity_parameters(1)*pi/180.d0
c = mat(i)%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))
mat(i)%plasticity_parameters(8)=alpha
mat(i)%plasticity_parameters(9)=k
if (iproc.eq.0) then
write(8,*) 'material ',i
write(8,*) ' phi =',phi
write(8,*) ' c =',c
write(8,*) ' alpha=',alpha
write(8,*) ' k =',k
end if
case ('DPII')
phi = mat(i)%plasticity_parameters(1)*pi/180.d0
c = mat(i)%plasticity_parameters(2)
alpha = tan(phi)/sqrt(9.d0+12.d0*(tan(phi))**2)
k = 3.d0*c /sqrt(9.d0+12.d0*(tan(phi))**2)
mat(i)%plasticity_parameters(8)=alpha
mat(i)%plasticity_parameters(9)=k
case ('MC')
phi = mat(i)%plasticity_parameters(1)*pi/180.d0
mat(i)%plasticity_parameters(8)=sin(phi)
mat(i)%plasticity_parameters(9)=cos(phi)
case ('Tresca')
case ('DPIII')
case default
call stop_run('error in compute_plastic_params$')
end select
if (iproc.eq.0) then
write(8,*) 'material ',i,'plasticity_type ',mat(i)%plasticity_type
write(8,*) 'plasticity_parameters',mat(i)%plasticity_parameters(:)
end if
end do
return
end