Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! VRM Feb. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure, &
plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb, &
straintot)
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 :: 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 (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
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
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
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')
! 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)*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)
! 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)*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
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
373
374
375
select case (trim(mat(i)%plasticity_type))
case ('No')
case ('vM')
case ('DPI')
case ('DPII')
case ('DPIII')
case ('DPIV')
case ('DPV')
case ('DPVI')
case ('DPVII')
case ('MC')
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