diff --git a/src/vrm.f90 b/src/vrm.f90 index 98011dc27f86136f19e543e8885b9dcb5969a2a8..9e67b2cc50f85188c51e578ee92067c162c96725 100644 --- a/src/vrm.f90 +++ b/src/vrm.f90 @@ -18,7 +18,7 @@ subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure, & plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb, & - straintot,call_cnt) + straintot) use constants use definitions @@ -42,7 +42,6 @@ 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 ))))))))))))) @@ -57,14 +56,6 @@ 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 @@ -111,17 +102,6 @@ select case (trim(plasticity_type)) 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. @@ -287,16 +267,11 @@ select case (trim(plasticity_type)) 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 + zeta = cos_theta-sin_phi*sin_theta/sqrt3 + yield = max(plasticity_parameters(2)/compressibility,(-pressure*sin_phi+c*cos_phi)/zeta) + fail = 2.d0*viscosity*e2d-yield 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 + viscosity=0.5d0*yield/e2d is_plastic=.true. else is_plastic=.false. @@ -325,14 +300,6 @@ select case (trim(plasticity_type)) 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 !-------------------------------------------------------------------------------