recasting vacancy dot state as microstructure state

code/constitutive.f90

@@ -567,8 +567,8 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el)
    LOCAL_DAMAGE_anisoBrittle_ID, &
    LOCAL_DAMAGE_anisoDuctile_ID, &
    LOCAL_DAMAGE_gurson_ID, &
-   LOCAL_DAMAGE_phaseField_ID
+   LOCAL_DAMAGE_phaseField_ID, &
-
+   LOCAL_VACANCY_generation_ID
  use plastic_titanmod, only: &
    plastic_titanmod_microstructure
  use plastic_nonlocal, only: &
@@ -589,6 +589,8 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el)
    damage_gurson_microstructure
  use damage_phaseField, only: &
    damage_phaseField_microstructure
+ use vacancy_generation, only: &
+   vacancy_generation_microstructure    
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -641,6 +643,15 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el)
 
  end select
 
+ select case (phase_damage(material_phase(ipc,ip,el)))
+   case (LOCAL_VACANCY_generation_ID)
+     call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,ipc, ip, el)
+     call vacancy_generation_microstructure(constitutive_homogenizedC(ipc,ip,el), Fe, &
+                                            nSlip,accumulatedSlip,constitutive_getTemperature(ipc,ip,el), &
+                                            subdt,ipc,ip,el)
+
+ end select
+
 end subroutine constitutive_microstructure
 
 
@@ -1084,8 +1095,7 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
    PLASTICITY_dislokmc_ID, &
    PLASTICITY_titanmod_ID, &
    PLASTICITY_nonlocal_ID, &
-   LOCAL_DAMAGE_gurson_ID, &
+   LOCAL_DAMAGE_gurson_ID
-   LOCAL_VACANCY_generation_ID
  use plastic_j2, only:  &
    plastic_j2_dotState
  use plastic_phenopowerlaw, only: &
@@ -1105,8 +1115,6 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
    totalNslip
  use damage_gurson, only: &
    damage_gurson_dotState
- use vacancy_generation, only: &
-   vacancy_generation_dotState
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -1128,8 +1136,6 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
    tick, tock, & 
    tickrate, &
    maxticks
- real(pReal), dimension(:), allocatable :: &
-   accumulatedSlip
  integer(pInt) :: &
    nSlip
  
@@ -1165,13 +1171,6 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
      call damage_gurson_dotState(Tstar_v, Lp, ipc, ip, el)
  end select
 
- select case (phase_vacancy(material_phase(ipc,ip,el)))
-   case (LOCAL_VACANCY_generation_ID)
-     call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,ipc,ip,el)
-     call vacancy_generation_dotState(nSlip,accumulatedSlip,Tstar_v,constitutive_getTemperature(ipc,ip,el), &
-                                      ipc, ip, el)
- end select
-
  if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) then
    call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
    !$OMP CRITICAL (debugTimingDotState)

code/vacancy_generation.f90

@@ -49,7 +49,7 @@ module vacancy_generation
    vacancy_generation_init, &
    vacancy_generation_stateInit, &
    vacancy_generation_aTolState, &
-   vacancy_generation_dotState, &
+   vacancy_generation_microstructure, &
    vacancy_generation_getLocalConcentration, &
    vacancy_generation_putLocalConcentration, &
    vacancy_generation_getConcentration, &
@@ -289,14 +289,20 @@ end subroutine vacancy_generation_aTolState
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-subroutine vacancy_generation_dotState(nSlip, accumulatedSlip, Tstar_v, Temperature, ipc, ip, el)
+subroutine vacancy_generation_microstructure(C, Fe, nSlip, accumulatedSlip, Temperature, subdt, &
+                                             ipc, ip, el)
  use material, only: &
    mappingConstitutive, &
    phase_vacancyInstance, &
    vacancyState
  use math, only : &
+   math_mul33x33, &
+   math_mul66x6, &
+   math_Mandel33to6, &
    math_Mandel6to33, &
-   math_trace33
+   math_transpose33, &
+   math_trace33, &
+   math_I3
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -304,31 +310,39 @@ subroutine vacancy_generation_dotState(nSlip, accumulatedSlip, Tstar_v, Temperat
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
    el                                                                                               !< element
+ real(pReal), intent(in) :: &
+   Fe(3,3), &
+   C (6,6)
  real(pReal), dimension(nSlip), intent(in) :: &
    accumulatedSlip
- real(pReal),  intent(in), dimension(6) :: &
-   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor (Mandel)
  real(pReal),  intent(in) :: &
    Temperature                                                                                      !< 2nd Piola Kirchhoff stress tensor (Mandel)
+ real(pReal),  intent(in) :: &
+   subdt
  real(pReal) :: &
    pressure, &
-   energyBarrier                                                                     
+   energyBarrier, &
+   stress(6), &
+   strain(6)                                                                     
  integer(pInt) :: &
    instance, phase, constituent 
 
  phase = mappingConstitutive(2,ipc,ip,el)
  constituent = mappingConstitutive(1,ipc,ip,el)
  instance = phase_vacancyInstance(phase)
- pressure = math_trace33(math_Mandel6to33(Tstar_v))
+
- energyBarrier = (vacancy_generation_formationEnergy(instance) - &
+ strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
-                  pressure)*vacancy_generation_atomicVol(instance) - &
+ stress = math_mul66x6(C,strain) 
+ pressure = math_trace33(math_Mandel6to33(stress))
+ energyBarrier = (vacancy_generation_formationEnergy(instance) - pressure)* &
+                 vacancy_generation_atomicVol(instance) - &
                  sum(accumulatedSlip)*vacancy_generation_plasticityCoeff(instance)
 
- vacancyState(phase)%dotState(1,constituent) = &
+ vacancyState(phase)%state(1,constituent) = &
-    vacancy_generation_freq(instance)* &
+   vacancyState(phase)%subState0(1,constituent) + &
-    exp(-energyBarrier/(kB*Temperature))
+   subdt*vacancy_generation_freq(instance)*exp(-energyBarrier/(kB*Temperature))
 
-end subroutine vacancy_generation_dotState
+end subroutine vacancy_generation_microstructure
 
 !--------------------------------------------------------------------------------------------------
 !> @brief returns vacancy concentration based on state layout