constitutive_plastic_isotropic.f90

Go to the documentation of this file.

# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_isotropic.f90"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_isotropic.f90"
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_isotropic
 
  type :: tparameters
    real(pReal) :: &
      M, &                                                                                          !< Taylor factor
      dot_gamma_0, &                                                                                !< reference strain rate
      n, &                                                                                          !< stress exponent
      h0, &
      h_ln, &
      xi_inf, &                                                                                     !< maximum critical stress
      a, &
      c_1, &
      c_4, &
      c_3, &
      c_2
    integer :: &
      of_debug = 0
    logical :: &
      dilatation
    character(len=pStringLen), allocatable, dimension(:) :: &
      output
  end type tparameters
 
  type :: tisotropicstate
    real(pReal), pointer, dimension(:) :: &
      xi, &
      gamma
  end type tisotropicstate
 
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
 type(tParameters),     allocatable, dimension(:) :: param
 type(tIsotropicState), allocatable, dimension(:) :: &
   dotState, &
   state
 
contains
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_init
 
  integer :: &
    Ninstance, &
    p, &
    NipcMyPhase, &
    sizeState, sizeDotState
  real(pReal) :: &
    xi_0
  character(len=pStringLen) :: &
    extmsg = ''
 
  write(6,'(/,a)') ' <<<+-  plastic_'//plasticity_isotropic_label//' init  -+>>>'; flush(6)
 
  write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
  write(6,'(a)')   ' https://doi.org/10.1016/j.scriptamat.2017.09.047'
 
  ninstance = count(phase_plasticity == plasticity_isotropic_id)
  if (iand(debug_level(debug_constitutive),debug_levelbasic) /= 0) &
    write(6,'(a16,1x,i5,/)') '# instances:',ninstance
 
  allocate(param(ninstance))
  allocate(state(ninstance))
  allocate(dotstate(ninstance))
 
  do p = 1, size(phase_plasticity)
    if (phase_plasticity(p) /= plasticity_isotropic_id) cycle
    associate(prm => param(phase_plasticityinstance(p)), &
              dot => dotstate(phase_plasticityinstance(p)), &
              stt => state(phase_plasticityinstance(p)), &
              config => config_phase(p))
 
    prm%output = config%getStrings('(output)',defaultval=emptystringarray)
 
 
 
 
 
 
    xi_0            = config%getFloat('tau0')
    prm%xi_inf      = config%getFloat('tausat')
    prm%dot_gamma_0 = config%getFloat('gdot0')
    prm%n           = config%getFloat('n')
    prm%h0          = config%getFloat('h0')
    prm%M           = config%getFloat('m')
    prm%h_ln        = config%getFloat('h0_slopelnrate', defaultval=0.0_preal)
    prm%c_1         = config%getFloat('tausat_sinhfita',defaultval=0.0_preal)
    prm%c_4         = config%getFloat('tausat_sinhfitb',defaultval=0.0_preal)
    prm%c_3         = config%getFloat('tausat_sinhfitc',defaultval=0.0_preal)
    prm%c_2         = config%getFloat('tausat_sinhfitd',defaultval=0.0_preal)
    prm%a           = config%getFloat('a')
 
    prm%dilatation  = config%keyExists('/dilatation/')
 
!--------------------------------------------------------------------------------------------------
!  sanity checks
    if (xi_0            <  0.0_preal) extmsg = trim(extmsg)//' xi_0'
    if (prm%dot_gamma_0 <= 0.0_preal) extmsg = trim(extmsg)//' dot_gamma_0'
    if (prm%n           <= 0.0_preal) extmsg = trim(extmsg)//' n'
    if (prm%a           <= 0.0_preal) extmsg = trim(extmsg)//' a'
    if (prm%M           <= 0.0_preal) extmsg = trim(extmsg)//' M'
 
!--------------------------------------------------------------------------------------------------
! allocate state arrays
    nipcmyphase = count(material_phaseat == p) * discretization_nip
    sizedotstate = size(['xi               ','accumulated_shear'])
    sizestate = sizedotstate
 
    call material_allocateplasticstate(p,nipcmyphase,sizestate,sizedotstate,0)
 
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
    stt%xi  => plasticstate(p)%state   (1,:)
    stt%xi  =  xi_0
    dot%xi  => plasticstate(p)%dotState(1,:)
    plasticstate(p)%atol(1) = config%getFloat('atol_xi',defaultval=1.0_preal)
    if (plasticstate(p)%atol(1) < 0.0_preal) extmsg = trim(extmsg)//' atol_xi'
 
    stt%gamma  => plasticstate(p)%state   (2,:)
    dot%gamma  => plasticstate(p)%dotState(2,:)
    plasticstate(p)%atol(2) = config%getFloat('atol_gamma',defaultval=1.0e-6_preal)
    if (plasticstate(p)%atol(2) < 0.0_preal) extmsg = trim(extmsg)//' atol_gamma'
    ! global alias
    plasticstate(p)%slipRate        => plasticstate(p)%dotState(2:2,:)
 
    plasticstate(p)%state0 = plasticstate(p)%state                                                  ! ToDo: this could be done centrally
 
    end associate
 
!--------------------------------------------------------------------------------------------------
!  exit if any parameter is out of range
    if (extmsg /= '') call io_error(211,ext_msg=trim(extmsg)//'('//plasticity_isotropic_label//')')
 
  enddo
 
end subroutine plastic_isotropic_init
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_lpanditstangent(lp,dlp_dmp,mp,instance,of)
 
  real(pReal), dimension(3,3),     intent(out) :: &
    Lp
  real(pReal), dimension(3,3,3,3), intent(out) :: &
    dLp_dMp
 
  real(pReal), dimension(3,3), intent(in) :: &
    Mp
  integer,                     intent(in) :: &
    instance, &
    of
 
  real(pReal), dimension(3,3) :: &
    Mp_dev
  real(pReal) :: &
    dot_gamma, &                                                                                    !< strainrate
    norm_Mp_dev, &                                                                                  !< norm of the deviatoric part of the Mandel stress
    squarenorm_Mp_dev
  integer :: &
    k, l, m, n
 
  associate(prm => param(instance), stt => state(instance))
 
  mp_dev = math_deviatoric33(mp)
  squarenorm_mp_dev = math_tensordot(mp_dev,mp_dev)
  norm_mp_dev = sqrt(squarenorm_mp_dev)
 
  if (norm_mp_dev > 0.0_preal) then
    dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_preal) * norm_mp_dev/(prm%M*stt%xi(of))) **prm%n
 
    lp = dot_gamma/prm%M * mp_dev/norm_mp_dev
# 194 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_isotropic.f90"
    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
      dlp_dmp(k,l,m,n) = (prm%n-1.0_preal) * mp_dev(k,l)*mp_dev(m,n) / squarenorm_mp_dev
    forall (k=1:3,l=1:3) &
      dlp_dmp(k,l,k,l) = dlp_dmp(k,l,k,l) + 1.0_preal
    forall (k=1:3,m=1:3) &
      dlp_dmp(k,k,m,m) = dlp_dmp(k,k,m,m) - 1.0_preal/3.0_preal
    dlp_dmp = dot_gamma / prm%M * dlp_dmp / norm_mp_dev
  else
    lp = 0.0_preal
    dlp_dmp = 0.0_preal
  end if
 
  end associate
 
end subroutine plastic_isotropic_lpanditstangent
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_lianditstangent(li,dli_dmi,mi,instance,of)
 
  real(pReal), dimension(3,3),     intent(out) :: &
    Li
  real(pReal), dimension(3,3,3,3), intent(out)  :: &
    dLi_dMi
 
  real(pReal), dimension(3,3), intent(in) :: &
    Mi
  integer,                     intent(in) :: &
    instance, &
    of
 
  real(pReal) :: &
    tr
  integer :: &
    k, l, m, n
 
  associate(prm => param(instance), stt => state(instance))
 
  tr=math_trace33(math_spherical33(mi))
 
  if (prm%dilatation .and. abs(tr) > 0.0_preal) then                                                 ! no stress or J2 plasticity --> Li and its derivative are zero
    li = math_i3 &
       * prm%dot_gamma_0/prm%M * (3.0_preal*prm%M*stt%xi(of))**(-prm%n) &
       * tr * abs(tr)**(prm%n-1.0_preal)
 
# 249 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_isotropic.f90"
 
    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
      dli_dmi(k,l,m,n) = prm%n / tr * li(k,l) * math_i3(m,n)
 
  else
    li      = 0.0_preal
    dli_dmi = 0.0_preal
  endif
 
  end associate
 
 end subroutine plastic_isotropic_lianditstangent
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_dotstate(mp,instance,of)
 
  real(pReal), dimension(3,3),  intent(in) :: &
    Mp
  integer,                      intent(in) :: &
    instance, &
    of
 
  real(pReal) :: &
    dot_gamma, &                                                                                    !< strainrate
    xi_inf_star, &                                                                                  !< saturation xi
    norm_Mp
 
  associate(prm => param(instance), stt => state(instance), dot => dotstate(instance))
 
  if (prm%dilatation) then
    norm_mp = sqrt(math_tensordot(mp,mp))
  else
    norm_mp = sqrt(math_tensordot(math_deviatoric33(mp),math_deviatoric33(mp)))
  endif
 
  dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_preal) * norm_mp /(prm%M*stt%xi(of))) **prm%n
 
  if (dot_gamma > 1e-12_preal) then
    if (deq0(prm%c_1)) then
      xi_inf_star = prm%xi_inf
    else
      xi_inf_star = prm%xi_inf &
                  + asinh( (dot_gamma / prm%c_1)**(1.0_preal / prm%c_2))**(1.0_preal / prm%c_3) &
                  / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_preal / prm%n)
    endif
    dot%xi(of) = dot_gamma &
               * ( prm%h0 + prm%h_ln * log(dot_gamma) ) &
               * abs( 1.0_preal - stt%xi(of)/xi_inf_star )**prm%a &
               * sign(1.0_preal, 1.0_preal - stt%xi(of)/xi_inf_star)
  else
    dot%xi(of) = 0.0_preal
  endif
 
  dot%gamma(of) = dot_gamma                                                                         ! ToDo: not really used
 
  end associate
 
end subroutine plastic_isotropic_dotstate
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_results(instance,group)
 
  integer,          intent(in) :: instance
  character(len=*), intent(in) :: group
 
  integer :: o
 
  associate(prm => param(instance), stt => state(instance))
  outputsloop: do o = 1,size(prm%output)
    select case(trim(prm%output(o)))
      case ('flowstress')                                                                           ! ToDo: should be 'xi'
        call results_writedataset(group,stt%xi,'xi','resistance against plastic flow','Pa')
    end select
  enddo outputsloop
  end associate
 
end subroutine plastic_isotropic_results
 
 
end submodule plastic_isotropic