diff --git a/src/vrm.f90 b/src/vrm.f90
index 03323d5de26a1cc4c2499987da05a89ecfc959c4..d6ac9ae9b8332d1de74481efc87aa9b071b2218d 100644
--- a/src/vrm.f90
+++ b/src/vrm.f90
@@ -72,120 +72,245 @@ 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.
+ 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
+ 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
- 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
+ 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')
- !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 ('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 ('MC')
- if (params%init_e2d) then
- theta = 0.d0
- else
- theta = lode_angle(J2d,J3d)
+ 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
- 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
+ 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)/180.d0*3.141592654d0
- if (straintot.gt.strain_soft_in) then
- phi=atan2(plasticity_parameters(8),plasticity_parameters(9))
+ 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
- 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
+ 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 ('Tresca')
- if (params%init_e2d) then
- theta = 0.d0
- else
- theta = lode_angle(J2d,J3d)
+ 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
- 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
+ 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.
- case default
- print *,plasticity_type
- call stop_run('error in vrm$')
- flag_vrm_pb = .true.
end select
!write (*,*) 'vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv'
@@ -217,44 +342,38 @@ 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
+ 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
end do
return
-end
\ No newline at end of file
+end