From 6411d366330c57ca52835c62402517df3e40beb1 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Fri, 16 Jan 2015 17:04:01 +0000 Subject: [PATCH] updated vacancy-damage model --- code/constitutive.f90 | 23 ++++----- code/crystallite.f90 | 4 +- code/damage_phaseField.f90 | 33 +++++++++---- code/homogenization.f90 | 55 +++++++++++++++++---- code/numerics.f90 | 16 +++---- code/vacancy_generation.f90 | 95 +++++++++++++++++-------------------- 6 files changed, 132 insertions(+), 94 deletions(-) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 7d05e07ae..1747b7080 100644 --- 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 @@ -1827,15 +1828,12 @@ function constitutive_getVacancyMobility33(ipc, ip, el) el !< element number 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), & - 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 !-------------------------------------------------------------------------------------------------- diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 1ca40bb52..cb9d2bddf 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -634,8 +634,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) damageState( mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) thermalState(mappingConstitutive(2,g,i,e))%subState0( :,mappingConstitutive(1,g,i,e)) = & thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) - vacancyState( mappingConstitutive(2,g,i,e))%subState0( :,mappingConstitutive(1,g,i,e)) = & - vacancyState( mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) + vacancyState(mappingConstitutive(2,g,i,e))%subState0( :,mappingConstitutive(1,g,i,e)) = & + vacancyState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_partionedFp0(1:3,1:3,g,i,e) ! ...plastic def grad crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_partionedLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness diff --git a/code/damage_phaseField.f90 b/code/damage_phaseField.f90 index a7cdb4753..a7ea882e9 100644 --- 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,17 +327,17 @@ 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 !-------------------------------------------------------------------------------------------------- diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 8a5a2bf1b..fa2bdea15 100644 --- 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/ & - homogenization_Ngrains(mesh_element(3,el)) + 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 !-------------------------------------------------------------------------------------------------- diff --git a/code/numerics.f90 b/code/numerics.f90 index 58e5fc600..7a8e691eb 100644 --- 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 diff --git a/code/vacancy_generation.f90 b/code/vacancy_generation.f90 index 661d535f4..f65e693df 100644 --- 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 - - 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) + instance + + instance = phase_vacancyInstance(mappingConstitutive(2,ipc,ip,el)) + vacancy_generation_getVacancyEnergy = & + vacancy_generation_specificFormationEnergy(instance)/vacancy_generation_surfaceEnergy(instance) -end function vacancy_generation_getVacancyPotentialDrivingForce +end function vacancy_generation_getVacancyEnergy !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results