potential driving force for conservative cahn hilliard vacancy diffusion

code/constitutive.f90

@@ -48,6 +48,7 @@ module constitutive
    constitutive_getVacancyConcentration, &
    constitutive_getVacancyDiffusion33, &
    constitutive_getVacancyMobility33, &
+   constitutive_getVacancyPotentialDrivingForce, &
    constitutive_postResults
  
  private :: &
@@ -1730,8 +1731,6 @@ end function constitutive_getVacancyDiffusion33
 function constitutive_getVacancyMobility33(ipc, ip, el)
  use prec, only: &
    pReal
- use lattice, only: &
-   lattice_VacancyDiffusion33
  use material, only: &
    material_phase, &
    LOCAL_VACANCY_generation_ID, &
@@ -1762,6 +1761,42 @@ function constitutive_getVacancyMobility33(ipc, ip, el)
 
 end function constitutive_getVacancyMobility33
 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns vacancy chemical potential driving force
+!--------------------------------------------------------------------------------------------------
+real(pReal) function constitutive_getVacancyPotentialDrivingForce(ipc, ip, el)
+ use prec, only: &
+   pReal
+ use material, only: &
+   material_phase, &
+   LOCAL_VACANCY_generation_ID, &
+   phase_vacancy
+ use vacancy_generation, only: &
+   vacancy_generation_getVacancyPotentialDrivingForce
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), dimension(3,3) :: &
+   constitutive_getVacancyMobility33
+ real(pReal), dimension(:), allocatable :: &
+   accumulatedSlip
+ integer(pInt) :: &
+   nSlip
+ 
+ select case(phase_vacancy(material_phase(ipc,ip,el)))                                                   
+   case (LOCAL_VACANCY_generation_ID)
+    call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,ipc,ip,el)
+    constitutive_getVacancyMobility33 = &
+      vacancy_generation_getVacancyPotentialDrivingForce(constitutive_getDamage(ipc, ip, el), &
+                                                         ipc,ip,el)
+    
+ end select
+
+end function constitutive_getVacancyPotentialDrivingForce
+
 !--------------------------------------------------------------------------------------------------
 !> @brief returns accumulated slip on each system defined
 !--------------------------------------------------------------------------------------------------

code/homogenization.f90

@@ -74,8 +74,9 @@ module homogenization
    field_getThermalConductivity33, &
    field_getMassDensity, &
    field_getSpecificHeat, &
-   field_getVacancyMobility, &
+   field_getVacancyMobility33, &
    field_getVacancyDiffusion33, &
+   field_getVacancyPotentialDrivingForce, &
    materialpoint_postResults, &
    field_postResults
  private :: &
@@ -1104,8 +1105,6 @@ end function field_getDamageMobility
 function field_getVacancyDiffusion33(ip,el)
  use mesh, only: &
    mesh_element
- use lattice, only: &
-   lattice_vacancyDiffusion33
  use material, only: &
    material_homog, &
    field_vacancy_type, &
@@ -1140,20 +1139,59 @@ function field_getVacancyDiffusion33(ip,el)
                                homogenization_Ngrains(mesh_element(3,el))
 
 end function field_getVacancyDiffusion33
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Returns average mobility for vacancy field at each integration point 
 !--------------------------------------------------------------------------------------------------
-real(pReal) function field_getVacancyMobility(ip,el)
+function field_getVacancyMobility33(ip,el)
  use mesh, only: &
    mesh_element
- use lattice, only: &
-   lattice_vacancyMobility
  use material, only: &
-   material_phase, &
    material_homog, &
    field_vacancy_type, &
    FIELD_VACANCY_NONLOCAL_ID, &
    homogenization_Ngrains
+ use constitutive, only: &
+   constitutive_getVacancyMobility33
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), dimension(3,3) :: &
+   field_getVacancyMobility33
+ integer(pInt) :: &
+   ipc
+   
+ 
+ field_getVacancyMobility33 = 0.0_pReal
+                                                
+ select case(field_vacancy_type(material_homog(ip,el)))                                                   
+   case (FIELD_VACANCY_NONLOCAL_ID)
+     do ipc = 1, homogenization_Ngrains(mesh_element(3,el))
+       field_getVacancyMobility33 = field_getVacancyMobility33 + constitutive_getVacancyMobility33(ipc,ip,el)
+     enddo
+      
+ end select   
+
+ field_getVacancyMobility33 = field_getVacancyMobility33/ &
+                              homogenization_Ngrains(mesh_element(3,el))
+
+end function field_getVacancyMobility33
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Returns average driving for vacancy chemical potential at each integration point 
+!--------------------------------------------------------------------------------------------------
+real(pReal) function field_getVacancyPotentialDrivingForce(ip,el)
+ use mesh, only: &
+   mesh_element
+ use material, only: &
+   material_homog, &
+   field_vacancy_type, &
+   FIELD_VACANCY_NONLOCAL_ID, &
+   homogenization_Ngrains
+ use constitutive, only: &
+   constitutive_getVacancyPotentialDrivingForce
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -1163,20 +1201,21 @@ real(pReal) function field_getVacancyMobility(ip,el)
    ipc
    
  
- field_getVacancyMobility =0.0_pReal
+ field_getVacancyPotentialDrivingForce = 0.0_pReal
                                                 
  select case(field_vacancy_type(material_homog(ip,el)))                                                   
    case (FIELD_VACANCY_NONLOCAL_ID)
      do ipc = 1, homogenization_Ngrains(mesh_element(3,el))
-       field_getVacancyMobility = field_getVacancyMobility + lattice_VacancyMobility(material_phase(ipc,ip,el))
+       field_getVacancyPotentialDrivingForce = field_getVacancyPotentialDrivingForce + &
+         constitutive_getVacancyPotentialDrivingForce(ipc,ip,el)
      enddo
       
  end select   
 
- field_getVacancyMobility = field_getVacancyMobility/ &
+ field_getVacancyPotentialDrivingForce = field_getVacancyPotentialDrivingForce/ &
                                          homogenization_Ngrains(mesh_element(3,el))
 
-end function field_getVacancyMobility
+end function field_getVacancyPotentialDrivingForce
 
 !--------------------------------------------------------------------------------------------------
 !> @brief ToDo

code/vacancy_generation.f90

@@ -55,6 +55,7 @@ module vacancy_generation
    vacancy_generation_getConcentration, &
    vacancy_generation_getVacancyDiffusion33, &
    vacancy_generation_getVacancyMobility33, &
+   vacancy_generation_getVacancyPotentialDrivingForce, &
    vacancy_generation_postResults
 
 contains
@@ -427,6 +428,8 @@ end function vacancy_generation_getVacancyDiffusion33
 !> @brief returns generation vacancy mobility tensor 
 !--------------------------------------------------------------------------------------------------
 function vacancy_generation_getVacancyMobility33(nSlip,accumulatedSlip,temperature,ipc,ip,el)
+ use math, only: &
+   math_I3
  use material, only: &
    mappingConstitutive, &
    phase_vacancyInstance
@@ -451,6 +454,7 @@ function vacancy_generation_getVacancyMobility33(nSlip,accumulatedSlip,temperatu
  instance = phase_vacancyInstance(phase)
 
  vacancy_generation_getVacancyMobility33 = &
+   math_I3* &
    vacancy_generation_surfaceEnergy(instance)* &
    vacancy_generation_diffusionCoeff0(instance)* &
    exp(-vacancy_generation_migrationEnergy(instance)/(kB*temperature))/ &
@@ -458,6 +462,34 @@ function vacancy_generation_getVacancyMobility33(nSlip,accumulatedSlip,temperatu
     
 end function vacancy_generation_getVacancyMobility33
 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns generation vacancy mobility tensor 
+!--------------------------------------------------------------------------------------------------
+real(pReal) function vacancy_generation_getVacancyPotentialDrivingForce(damage,ipc,ip,el)
+ use material, only: &
+   mappingConstitutive, &
+   phase_vacancyInstance
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), intent(in) :: &
+   damage
+ integer(pInt) :: &
+   phase, constituent, instance
+ 
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ instance = phase_vacancyInstance(phase)
+
+ vacancy_generation_getVacancyPotentialDrivingForce = &
+   1.0_pReal - damage - damage*damage*vacancy_generation_formationEnergy(instance)/ &
+                                      vacancy_generation_surfaceEnergy(instance)
+    
+end function vacancy_generation_getVacancyPotentialDrivingForce
+
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of constitutive results
 !--------------------------------------------------------------------------------------------------