updated vacancy-damage model

2015-01-16 17:04:01 +00:00 · 2015-01-16 17:04:01 +00:00 · 6411d36633
parent 93e50366bb
commit 6411d36633
6 changed files with 132 additions and 94 deletions
--- a/code/constitutive.f90
+++ b/code/constitutive.f90
@ -49,7 +49,7 @@ module constitutive
   constitutive_getVacancyConcentration, &
   constitutive_getVacancyDiffusion33, &
   constitutive_getVacancyMobility33, &
-   constitutive_getVacancyPotentialDrivingForce, &
+   constitutive_getVacancyEnergy, &
   constitutive_postResults
 
 private :: &
@ -666,8 +666,9 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el)

 select case (phase_vacancy(material_phase(ipc,ip,el)))
   case (LOCAL_VACANCY_generation_ID)
-     call vacancy_generation_microstructure(constitutive_homogenizedC(ipc,ip,el), Fe, &
+     call vacancy_generation_microstructure(Tstar_v, &
                                            constitutive_getTemperature(ipc,ip,el), &
+                                            constitutive_getDamage(ipc, ip, el), &
                                            subdt,ipc,ip,el)

 end select
@ -1828,13 +1829,10 @@ function constitutive_getVacancyMobility33(ipc, ip, el)
 real(pReal), dimension(3,3) :: &
   constitutive_getVacancyMobility33
 
- integer(pInt) :: &
-   nSlip
- 
 select case(phase_vacancy(material_phase(ipc,ip,el)))                                                   
   case (LOCAL_VACANCY_generation_ID)
    constitutive_getVacancyMobility33 = &
-      vacancy_generation_getVacancyMobility33(nSlip,constitutive_getTemperature(ipc,ip,el), &
+      vacancy_generation_getVacancyMobility33(constitutive_getTemperature(ipc,ip,el), &
                                              ipc,ip,el)
    
 end select
@ -1844,7 +1842,7 @@ end function constitutive_getVacancyMobility33
 !--------------------------------------------------------------------------------------------------
 !> @brief returns vacancy chemical potential driving force
 !--------------------------------------------------------------------------------------------------
-real(pReal) function constitutive_getVacancyPotentialDrivingForce(ipc, ip, el)
+real(pReal) function constitutive_getVacancyEnergy(ipc, ip, el)
 use prec, only: &
   pReal
 use material, only: &
@ -1852,7 +1850,7 @@ real(pReal) function constitutive_getVacancyPotentialDrivingForce(ipc, ip, el)
   LOCAL_VACANCY_generation_ID, &
   phase_vacancy
 use vacancy_generation, only: &
-   vacancy_generation_getVacancyPotentialDrivingForce
+   vacancy_generation_getVacancyEnergy

 implicit none
 integer(pInt), intent(in) :: &
@ -1862,13 +1860,12 @@ real(pReal) function constitutive_getVacancyPotentialDrivingForce(ipc, ip, el)
 
 select case(phase_vacancy(material_phase(ipc,ip,el)))                                                   
   case (LOCAL_VACANCY_generation_ID)
-    constitutive_getVacancyPotentialDrivingForce = &
-      vacancy_generation_getVacancyPotentialDrivingForce(constitutive_getDamage(ipc, ip, el), &
-                                                         ipc,ip,el)
+    constitutive_getVacancyEnergy = &
+      vacancy_generation_getVacancyEnergy(ipc,ip,el)
    
 end select

-end function constitutive_getVacancyPotentialDrivingForce
+end function constitutive_getVacancyEnergy


 !--------------------------------------------------------------------------------------------------
--- a/code/damage_phaseField.f90
+++ b/code/damage_phaseField.f90
@ -28,7 +28,9 @@ module damage_phaseField
 real(pReal),                         dimension(:),           allocatable,         private :: &
   damage_phaseField_aTol, &
   damage_phaseField_surfaceEnergy, &
-   damage_phaseField_vacancyFormationEnergy
+   damage_phaseField_vacancyFormationEnergy, &
+   damage_phaseField_atomicVol, &
+   damage_phaseField_specificVacancyFormationEnergy

 enum, bind(c) 
   enumerator :: undefined_ID, &
@ -123,6 +125,8 @@ subroutine damage_phaseField_init(fileUnit)
 allocate(damage_phaseField_Noutput(maxNinstance),                             source=0_pInt) 
 allocate(damage_phaseField_surfaceEnergy(maxNinstance),                       source=0.0_pReal) 
 allocate(damage_phaseField_vacancyFormationEnergy(maxNinstance),              source=0.0_pReal) 
+ allocate(damage_phaseField_atomicVol(maxNinstance),                           source=0.0_pReal) 
+ allocate(damage_phaseField_specificVacancyFormationEnergy(maxNinstance),      source=0.0_pReal) 
 allocate(damage_phaseField_aTol(maxNinstance),                                source=0.0_pReal) 

 rewind(fileUnit)
@ -162,6 +166,9 @@ subroutine damage_phaseField_init(fileUnit)
       case ('vacancyformationenergy')
         damage_phaseField_vacancyFormationEnergy(instance) = IO_floatValue(line,positions,2_pInt)

+       case ('atomicvolume')
+         damage_phaseField_atomicVol(instance) = IO_floatValue(line,positions,2_pInt)
+         
       case ('atol_damage')
         damage_phaseField_aTol(instance) = IO_floatValue(line,positions,2_pInt)

@ -186,6 +193,12 @@ subroutine damage_phaseField_init(fileUnit)
   if (phase_damage(phase) == LOCAL_damage_phaseField_ID) then
     NofMyPhase=count(material_phase==phase)
     instance = phase_damageInstance(phase)
+
+!--------------------------------------------------------------------------------------------------
+!  pre-calculating derived material parameters
+     damage_phaseField_specificVacancyFormationEnergy(instance) = &
+       damage_phaseField_vacancyFormationEnergy(instance)/damage_phaseField_atomicVol(instance)
+
 !--------------------------------------------------------------------------------------------------
 !  Determine size of postResults array
     outputsLoop: do o = 1_pInt,damage_phaseField_Noutput(instance)
@ -199,6 +212,7 @@ subroutine damage_phaseField_init(fileUnit)
          damage_phaseField_sizePostResults(instance)  = damage_phaseField_sizePostResults(instance) + mySize
       endif
     enddo outputsLoop
+
 ! Determine size of state array
     sizeDotState              =   0_pInt
     sizeState                 =   2_pInt
@ -303,7 +317,8 @@ subroutine damage_phaseField_microstructure(C, Fe, Cv, subdt, ipc, ip, el)
   phase, constituent, instance
 real(pReal) :: &
   strain(6), &
-   stress(6)
+   stress(6), &
+   drivingForce

 phase = mappingConstitutive(2,ipc,ip,el)
 constituent = mappingConstitutive(1,ipc,ip,el)
@ -312,16 +327,16 @@ subroutine damage_phaseField_microstructure(C, Fe, Cv, subdt, ipc, ip, el)
 strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
 stress = math_mul66x6(C,strain) 
 
+ drivingForce = (1.0_pReal - Cv)*(1.0_pReal - Cv) + &
+                (Cv*damage_phaseField_specificVacancyFormationEnergy(instance) + &
+                 sum(abs(stress*strain)))/damage_phaseField_surfaceEnergy(instance)           
 damageState(phase)%state(2,constituent) = &
-   max(residualStiffness, &
-       min(damageState(phase)%state0(2,constituent), &
-           (1.0_pReal - Cv)*damage_phaseField_surfaceEnergy(instance)/ &
-           (2.0_pReal*(sum(abs(stress*strain)) + Cv*damage_phaseField_vacancyFormationEnergy(instance)))))
+   (1.0_pReal - Cv)*(1.0_pReal - Cv)/drivingForce
       
 damageState(phase)%state(1,constituent) = &
   damageState(phase)%state(2,constituent) + &
   (damageState(phase)%subState0(1,constituent) - damageState(phase)%state(2,constituent))* &
-   exp(-subdt/(damageState(phase)%state(2,constituent)*lattice_DamageMobility(phase)))
+   exp(-2.0_pReal*subdt*drivingForce/lattice_DamageMobility(phase))

 end subroutine damage_phaseField_microstructure
 
--- a/code/homogenization.f90
+++ b/code/homogenization.f90
@ -64,6 +64,7 @@ module homogenization
   homogenization_init, &
   materialpoint_stressAndItsTangent, &
   field_getLocalDamage, &
+   field_getFieldDamage, &
   field_putFieldDamage, &
   field_getLocalTemperature, &
   field_putFieldTemperature, &
@ -77,7 +78,7 @@ module homogenization
   field_getSpecificHeat, &
   field_getVacancyMobility33, &
   field_getVacancyDiffusion33, &
-   field_getVacancyPotentialDrivingForce, &
+   field_getVacancyEnergy, &
   materialpoint_postResults, &
   field_postResults
 private :: &
@ -1161,6 +1162,8 @@ function field_getVacancyMobility33(ip,el)
   field_vacancy_type, &
   FIELD_VACANCY_NONLOCAL_ID, &
   homogenization_Ngrains
+ use crystallite, only: &
+   crystallite_push33ToRef
 use constitutive, only: &
   constitutive_getVacancyMobility33

@ -1179,7 +1182,9 @@ function field_getVacancyMobility33(ip,el)
 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)
+       field_getVacancyMobility33 = field_getVacancyMobility33 + &
+        crystallite_push33ToRef(ipc,ip,el, &
+                                constitutive_getVacancyMobility33(ipc,ip,el))
     enddo
      
 end select   
@ -1192,7 +1197,7 @@ end function field_getVacancyMobility33
 !--------------------------------------------------------------------------------------------------
 !> @brief Returns average driving for vacancy chemical potential at each integration point 
 !--------------------------------------------------------------------------------------------------
-real(pReal) function field_getVacancyPotentialDrivingForce(ip,el)
+real(pReal) function field_getVacancyEnergy(ip,el)
 use mesh, only: &
   mesh_element
 use material, only: &
@ -1201,7 +1206,7 @@ real(pReal) function field_getVacancyPotentialDrivingForce(ip,el)
   FIELD_VACANCY_NONLOCAL_ID, &
   homogenization_Ngrains
 use constitutive, only: &
-   constitutive_getVacancyPotentialDrivingForce
+   constitutive_getVacancyEnergy

 implicit none
 integer(pInt), intent(in) :: &
@ -1211,21 +1216,21 @@ real(pReal) function field_getVacancyPotentialDrivingForce(ip,el)
   ipc
   
 
- field_getVacancyPotentialDrivingForce = 0.0_pReal
+ field_getVacancyEnergy = 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_getVacancyPotentialDrivingForce = field_getVacancyPotentialDrivingForce + &
-         constitutive_getVacancyPotentialDrivingForce(ipc,ip,el)
+       field_getVacancyEnergy = field_getVacancyEnergy + &
+         constitutive_getVacancyEnergy(ipc,ip,el)
     enddo
      
 end select   

- field_getVacancyPotentialDrivingForce = field_getVacancyPotentialDrivingForce/ &
+ field_getVacancyEnergy = field_getVacancyEnergy/ &
                          homogenization_Ngrains(mesh_element(3,el))

-end function field_getVacancyPotentialDrivingForce
+end function field_getVacancyEnergy

 !--------------------------------------------------------------------------------------------------
 !> @brief ToDo
@ -1257,6 +1262,36 @@ real(pReal) function field_getLocalDamage(ip,el)

 end function field_getLocalDamage

+!--------------------------------------------------------------------------------------------------
+!> @brief ToDo
+!--------------------------------------------------------------------------------------------------
+real(pReal) function field_getFieldDamage(ip,el)
+ use mesh, only: &
+   mesh_element
+ use material, only: &
+   homogenization_Ngrains
+ use constitutive, only: &
+   constitutive_getDamage
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ integer(pInt) :: &
+   ipc
+
+!--------------------------------------------------------------------------------------------------
+! computing the damage value needed to be passed to field solver
+ field_getFieldDamage = 0.0_pReal
+                                                
+ do ipc = 1, homogenization_Ngrains(mesh_element(3,el))
+   field_getFieldDamage = field_getFieldDamage + constitutive_getDamage(ipc,ip,el)
+ enddo
+
+ field_getFieldDamage = field_getFieldDamage/homogenization_Ngrains(mesh_element(3,el))
+
+end function field_getFieldDamage
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Sets the regularised damage value in field state
 !--------------------------------------------------------------------------------------------------
--- a/code/numerics.f90
+++ b/code/numerics.f90
@ -150,14 +150,14 @@ module numerics
                             &-thermal_mg_levels_ksp_type chebyshev &
                             &-thermal_mg_levels_ksp_chebyshev_estimate_eigenvalues 0,0.1,0,1.1 &
                             &-thermal_mg_levels_pc_type sor &
-                             &-vacancyDiffusion_snes_type newtonls &
-                             &-vacancyDiffusion_snes_linesearch_type cp &
-                             &-vacancyDiffusion_ksp_type fgmres &
-                             &-vacancyDiffusion_snes_atol 1e-6 &
-                             &-vacancyDiffusion_pc_type ml &
-                             &-vacancyDiffusion_mg_levels_ksp_type chebyshev &
-                             &-vacancyDiffusion_mg_levels_ksp_chebyshev_estimate_eigenvalues 0,0.1,0,1.1 &
-                             &-vacancyDiffusion_mg_levels_pc_type sor '
+                             &-vacancy_snes_type newtonls &
+                             &-vacancy_snes_linesearch_type cp &
+                             &-vacancy_ksp_type fgmres &
+                             &-vacancy_snes_atol 1e-9 &
+                             &-vacancy_pc_type ml &
+                             &-vacancy_mg_levels_ksp_type chebyshev &
+                             &-vacancy_mg_levels_ksp_chebyshev_estimate_eigenvalues 0,0.1,0,1.1 &
+                             &-vacancy_mg_levels_pc_type sor '
 integer(pInt), protected, public :: &
   itmaxFEM                   =  25_pInt, &                                                         !< maximum number of iterations
   itminFEM                   =  1_pInt, &                                                          !< minimum number of iterations
--- a/code/vacancy_generation.f90
+++ b/code/vacancy_generation.f90
@ -24,18 +24,20 @@ module vacancy_generation
 integer(pInt),                       dimension(:),           allocatable, target, public :: &
   vacancy_generation_Noutput                                                                   !< number of outputs per instance of this damage 

- real(pReal),                         dimension(:),           allocatable,         public :: &
+ real(pReal),                         dimension(:),           allocatable,        private :: &
   vacancy_generation_aTol, &
   vacancy_generation_freq, &
   vacancy_generation_formationEnergy, &
+   vacancy_generation_specificFormationEnergy, &
   vacancy_generation_migrationEnergy, &
   vacancy_generation_diffusionCoeff0, &                                                        !< the temperature-independent diffusion coefficient D_0
   vacancy_generation_atomicVol, &
   vacancy_generation_surfaceEnergy, &
-   vacancy_generation_plasticityCoeff
+   vacancy_generation_plasticityCoeff, &
+   vacancy_generation_kBCoeff

 real(pReal),                                                 parameter,           private :: &
-   kB = 1.38e-23_pReal                                                                          !< Boltzmann constant in J/Kelvin
+   kB = 1.3806488e-23_pReal                                                                          !< Boltzmann constant in J/Kelvin

 enum, bind(c) 
   enumerator :: undefined_ID, &
@ -55,7 +57,7 @@ module vacancy_generation
   vacancy_generation_getConcentration, &
   vacancy_generation_getVacancyDiffusion33, &
   vacancy_generation_getVacancyMobility33, &
-   vacancy_generation_getVacancyPotentialDrivingForce, &
+   vacancy_generation_getVacancyEnergy, &
   vacancy_generation_postResults

 contains
@ -134,10 +136,13 @@ subroutine vacancy_generation_init(fileUnit)
 allocate(vacancy_generation_aTol(maxNinstance),                                source=0.0_pReal) 
 allocate(vacancy_generation_freq(maxNinstance),                                source=0.0_pReal) 
 allocate(vacancy_generation_formationEnergy(maxNinstance),                     source=0.0_pReal) 
+ allocate(vacancy_generation_specificFormationEnergy(maxNinstance),             source=0.0_pReal) 
 allocate(vacancy_generation_migrationEnergy(maxNinstance),                     source=0.0_pReal) 
+ allocate(vacancy_generation_diffusionCoeff0(maxNinstance),                     source=0.0_pReal) 
 allocate(vacancy_generation_atomicVol(maxNinstance),                           source=0.0_pReal) 
 allocate(vacancy_generation_surfaceEnergy(maxNinstance),                       source=0.0_pReal)
 allocate(vacancy_generation_plasticityCoeff(maxNinstance),                     source=0.0_pReal)
+ allocate(vacancy_generation_kBCoeff(maxNinstance),                             source=0.0_pReal)

 rewind(fileUnit)
 phase = 0_pInt
@ -175,7 +180,7 @@ subroutine vacancy_generation_init(fileUnit)
       case ('atolvacancygeneration')
         vacancy_generation_aTol(instance) = IO_floatValue(line,positions,2_pInt)

-       case ('debyefrequency')
+       case ('vacancyformationfreq')
         vacancy_generation_freq(instance) = IO_floatValue(line,positions,2_pInt)

       case ('vacancyformationenergy')
@ -205,6 +210,12 @@ subroutine vacancy_generation_init(fileUnit)
     NofMyPhase=count(material_phase==phase)
     instance = phase_vacancyInstance(phase)

+!--------------------------------------------------------------------------------------------------
+!  pre-calculating derived material parameters
+     vacancy_generation_kBCoeff(instance) = kB/vacancy_generation_atomicVol(instance)
+     vacancy_generation_specificFormationEnergy(instance) = &
+       vacancy_generation_formationEnergy(instance)/vacancy_generation_atomicVol(instance)
+     
 !--------------------------------------------------------------------------------------------------
 !  Determine size of postResults array
     outputsLoop: do o = 1_pInt,vacancy_generation_Noutput(instance)
@ -263,10 +274,9 @@ subroutine vacancy_generation_stateInit(phase)
 integer(pInt), intent(in) :: phase                                                    !< number specifying the phase of the vacancy
 real(pReal), dimension(vacancyState(phase)%sizeState) :: tempState

- tempState = lattice_equilibriumVacancyConcentration(phase)
- vacancyState(phase)%state = spread(tempState,2,size(vacancyState(phase)%state(1,:)))
- vacancyState(phase)%state0 = vacancyState(phase)%state
- vacancyState(phase)%partionedState0 = vacancyState(phase)%state
+ tempState(1) = lattice_equilibriumVacancyConcentration(phase)
+ vacancyState(phase)%state0 = spread(tempState,2,size(vacancyState(phase)%state(1,:)))
+
 end subroutine vacancy_generation_stateInit

 !--------------------------------------------------------------------------------------------------
@ -289,7 +299,7 @@ end subroutine vacancy_generation_aTolState
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-subroutine vacancy_generation_microstructure(C, Fe, Temperature, subdt, &
+subroutine vacancy_generation_microstructure(Tstar_v, temperature, damage, subdt, &
                                             ipc, ip, el)
 use material, only: &
   mappingConstitutive, &
@ -297,13 +307,8 @@ subroutine vacancy_generation_microstructure(C, Fe, Temperature, subdt, &
   plasticState, &
   vacancyState
 use math, only : &
-   math_mul33x33, &
-   math_mul66x6, &
-   math_Mandel33to6, &
   math_Mandel6to33, &
-   math_transpose33, &
-   math_trace33, &
-   math_I3
+   math_trace33

 implicit none
 integer(pInt), intent(in) :: &
@ -311,17 +316,13 @@ subroutine vacancy_generation_microstructure(C, Fe, Temperature, subdt, &
   ip, &                                                                                            !< integration point
   el                                                                                               !< element
 real(pReal), intent(in) :: &
-   Fe(3,3), &
-   C (6,6)
- real(pReal),  intent(in) :: &
-   Temperature                                                                                      !< 2nd Piola Kirchhoff stress tensor (Mandel)
- real(pReal),  intent(in) :: &
+   Tstar_v(6), &
+   temperature, &                                                                                   !< 2nd Piola Kirchhoff stress tensor (Mandel)
+   damage, &
   subdt
 real(pReal) :: &
   pressure, &
-   energyBarrier, &
-   stress(6), &
-   strain(6)                                                                     
+   stressBarrier
 integer(pInt) :: &
   instance, phase, constituent 

@ -329,12 +330,12 @@ subroutine vacancy_generation_microstructure(C, Fe, Temperature, subdt, &
 constituent = mappingConstitutive(1,ipc,ip,el)
 instance = phase_vacancyInstance(phase)

- strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
- stress = math_mul66x6(C,strain) 
- pressure = math_trace33(math_Mandel6to33(stress))
- energyBarrier = (vacancy_generation_formationEnergy(instance) - pressure)* &
-                 vacancy_generation_atomicVol(instance) - &
-                 sum(plasticState(phase)%accumulatedSlip(:,constituent))*vacancy_generation_plasticityCoeff(instance)
+ pressure = math_trace33(math_Mandel6to33(Tstar_v))/3.0_pReal
+ stressBarrier = max(0.0_pReal, &
+                     vacancy_generation_specificFormationEnergy(instance) - &
+                     pressure - &
+                     vacancy_generation_plasticityCoeff(instance)* &
+                     sum(plasticState(phase)%accumulatedSlip(:,constituent))) 

 vacancyState(phase)%state(1,constituent) = &
   vacancyState(phase)%subState0(1,constituent) + &
@ -439,16 +440,16 @@ end function vacancy_generation_getVacancyDiffusion33
 !--------------------------------------------------------------------------------------------------
 !> @brief returns generation vacancy mobility tensor 
 !--------------------------------------------------------------------------------------------------
-function vacancy_generation_getVacancyMobility33(nSlip,temperature,ipc,ip,el)
+function vacancy_generation_getVacancyMobility33(temperature,ipc,ip,el)
 use math, only: &
   math_I3
 use material, only: &
   mappingConstitutive, &
-   phase_vacancyInstance
+   phase_vacancyInstance, &
+   plasticState

 implicit none
 integer(pInt), intent(in) :: &
-   nSlip, &
   ipc, &                                                                                           !< grain number
   ip, &                                                                                            !< integration point number
   el                                                                                               !< element number
@ -464,18 +465,14 @@ function vacancy_generation_getVacancyMobility33(nSlip,temperature,ipc,ip,el)
 instance = phase_vacancyInstance(phase)

 vacancy_generation_getVacancyMobility33 = &
-   math_I3* &
-   vacancy_generation_surfaceEnergy(instance)* &
-   vacancy_generation_diffusionCoeff0(instance)* &
-   exp(-vacancy_generation_migrationEnergy(instance)/(kB*temperature))/ &
-   (kB*temperature)
+   math_I3*(1.0_pReal + sum(plasticState(phase)%accumulatedSlip(:,constituent)))
       
 end function vacancy_generation_getVacancyMobility33

 !--------------------------------------------------------------------------------------------------
 !> @brief returns generation vacancy mobility tensor 
 !--------------------------------------------------------------------------------------------------
-real(pReal) function vacancy_generation_getVacancyPotentialDrivingForce(damage,ipc,ip,el)
+real(pReal) function vacancy_generation_getVacancyEnergy(ipc,ip,el)
 use material, only: &
   mappingConstitutive, &
   phase_vacancyInstance
@ -485,20 +482,14 @@ real(pReal) function vacancy_generation_getVacancyPotentialDrivingForce(damage,i
   ipc, &                                                                                           !< grain number
   ip, &                                                                                            !< integration point number
   el                                                                                               !< element number
- real(pReal), intent(in) :: &
-   damage
 integer(pInt) :: &
-   phase, constituent, instance
+   instance
   
- phase = mappingConstitutive(2,ipc,ip,el)
- constituent = mappingConstitutive(1,ipc,ip,el)
- instance = phase_vacancyInstance(phase)
+ instance = phase_vacancyInstance(mappingConstitutive(2,ipc,ip,el))   
+ vacancy_generation_getVacancyEnergy = &
+   vacancy_generation_specificFormationEnergy(instance)/vacancy_generation_surfaceEnergy(instance)
    
- vacancy_generation_getVacancyPotentialDrivingForce = &
-   1.0_pReal - damage - damage*damage*vacancy_generation_formationEnergy(instance)/ &
-                                      vacancy_generation_surfaceEnergy(instance)
-    
-end function vacancy_generation_getVacancyPotentialDrivingForce
+end function vacancy_generation_getVacancyEnergy

 !--------------------------------------------------------------------------------------------------
 !> @brief return array of constitutive results