Updated Drucker-Prager yield options

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