Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! VRM Feb. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure, &
plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb, &
straintot,call_cnt)
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
integer call_cnt
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
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
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=plasticity_parameters(3)
if (strain_soft_in.gt.0.d0) then
strain_soft_out=plasticity_parameters(4)
strain_soft_phi=plasticity_parameters(5)*pi/180.d0
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
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
!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.
else
is_plastic=.false.
end if
case ('DPII')
! Strain softening/hardening
phi = plasticity_parameters(1)*pi/180.d0
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)*pi/180.d0
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
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
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=plasticity_parameters(3)
if (strain_soft_in.gt.0.d0) then
strain_soft_out=plasticity_parameters(4)
strain_soft_phi=plasticity_parameters(5)*pi/180.d0
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
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)
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=plasticity_parameters(3)
if (strain_soft_in.gt.0.d0) then
strain_soft_out=plasticity_parameters(4)
strain_soft_phi=plasticity_parameters(5)*pi/180.d0
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
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)
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)
yield = max(k/compressibility,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
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)*pi/180.d0
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
c = plasticity_parameters(2)
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
is_plastic=.true.
else
is_plastic=.false.
end if
case ('DPVII')
Dave Whipp
committed
! Strain softening/hardening
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)*pi/180.d0
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
c = plasticity_parameters(2)
sin_phi=sin(phi)
cos_phi=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
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
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
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)*pi/180.d0
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
c = plasticity_parameters(2)
sin_phi = sin(phi)
cos_phi = cos(phi)
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.
!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')
Dave Whipp
committed
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
Dave Whipp
committed
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
Dave Whipp
committed
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)
k = (6.d0*c*cos(phi))/(sqrt(3.d0)*3.d0)
mat(i)%plasticity_parameters(8) = alpha
mat(i)%plasticity_parameters(9) = k
Dave Whipp
committed
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.d0))
k = (3.d0*c)/(sqrt(9.d0+12.d0*(tan(phi))**2.d0))
mat(i)%plasticity_parameters(8) = alpha
mat(i)%plasticity_parameters(9) = k
case ('DPV')
case ('DPVI')
Dave Whipp
committed
phi = mat(i)%plasticity_parameters(1)*pi/180.d0
mat(i)%plasticity_parameters(8) = tan(phi)
Dave Whipp
committed
phi = mat(i)%plasticity_parameters(1)*pi/180.d0
mat(i)%plasticity_parameters(8) = sin(phi)
mat(i)%plasticity_parameters(9) = cos(phi)
Dave Whipp
committed
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 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