From 20437b2ae09d45328d85bcd2b5728b97620ecea7 Mon Sep 17 00:00:00 2001
From: Pratheek Shanthraj
Date: Tue, 25 Nov 2014 17:23:37 +0000
Subject: [PATCH] added brittle(elastic energy release rate)/ductile(vacancy
condensation) phase field damage model coupled to vacancy concentration
simplified vacancy_generation
---
code/Makefile | 5 +-
code/constitutive.f90 | 138 +++++++---
code/damage_phaseField.f90 | 520 ++++++++++++++++++++++++++++++++++++
code/material.f90 | 7 +-
code/vacancy_generation.f90 | 235 +++++++---------
5 files changed, 735 insertions(+), 170 deletions(-)
create mode 100644 code/damage_phaseField.f90
diff --git a/code/Makefile b/code/Makefile
index f4d6d1c19..4d5816ee2 100644
--- a/code/Makefile
+++ b/code/Makefile
@@ -326,7 +326,7 @@ COMPILE =$(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATI
COMPILE_MAXOPTI =$(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90))
###################################################################################################
DAMAGE_FILES = \
- damage_none.o damage_isoBrittle.o damage_isoDuctile.o damage_gurson.o damage_anisoBrittle.o damage_anisoDuctile.o
+ damage_none.o damage_isoBrittle.o damage_isoDuctile.o damage_gurson.o damage_anisoBrittle.o damage_anisoDuctile.o damage_phaseField.o
THERMAL_FILES = \
thermal_isothermal.o thermal_adiabatic.o
@@ -496,6 +496,9 @@ damage_anisoBrittle.o: damage_anisoBrittle.f90 \
damage_anisoDuctile.o: damage_anisoDuctile.f90 \
lattice.o
+damage_phaseField.o: damage_phaseField.f90 \
+ lattice.o
+
damage_gurson.o: damage_gurson.f90 \
lattice.o
diff --git a/code/constitutive.f90 b/code/constitutive.f90
index 4fb19b0ab..e42b73923 100644
--- a/code/constitutive.f90
+++ b/code/constitutive.f90
@@ -47,6 +47,7 @@ module constitutive
constitutive_putLocalVacancyConcentration, &
constitutive_getVacancyConcentration, &
constitutive_getVacancyDiffusion33, &
+ constitutive_getVacancyMobility33, &
constitutive_postResults
private :: &
@@ -131,6 +132,7 @@ subroutine constitutive_init(temperature_init)
LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_DAMAGE_anisoDuctile_ID, &
LOCAL_DAMAGE_gurson_ID, &
+ LOCAL_DAMAGE_phaseField_ID, &
LOCAL_THERMAL_isothermal_ID, &
LOCAL_THERMAL_adiabatic_ID, &
LOCAL_VACANCY_constant_ID, &
@@ -141,6 +143,7 @@ subroutine constitutive_init(temperature_init)
LOCAL_DAMAGE_anisoBrittle_LABEL, &
LOCAL_DAMAGE_anisoDuctile_LABEL, &
LOCAL_DAMAGE_gurson_LABEL, &
+ LOCAL_DAMAGE_phaseField_label, &
LOCAL_THERMAL_isothermal_label, &
LOCAL_THERMAL_adiabatic_label, &
LOCAL_VACANCY_constant_label, &
@@ -165,6 +168,7 @@ subroutine constitutive_init(temperature_init)
use damage_anisoDuctile
use damage_anisoBrittle
use damage_gurson
+ use damage_phaseField
use thermal_isothermal
use thermal_adiabatic
use vacancy_constant
@@ -211,6 +215,7 @@ subroutine constitutive_init(temperature_init)
if (any(phase_damage == LOCAL_DAMAGE_anisoBrittle_ID)) call damage_anisoBrittle_init(FILEUNIT)
if (any(phase_damage == LOCAL_DAMAGE_anisoductile_ID)) call damage_anisoDuctile_init(FILEUNIT)
if (any(phase_damage == LOCAL_DAMAGE_gurson_ID)) call damage_gurson_init(FILEUNIT)
+ if (any(phase_damage == LOCAL_DAMAGE_phaseField_ID)) call damage_phaseField_init(FILEUNIT)
close(FILEUNIT)
!--------------------------------------------------------------------------------------------------
@@ -324,6 +329,11 @@ subroutine constitutive_init(temperature_init)
thisNoutput => damage_gurson_Noutput
thisOutput => damage_gurson_output
thisSize => damage_gurson_sizePostResult
+ case (LOCAL_DAMAGE_phaseField_ID)
+ outputName = LOCAL_DAMAGE_phaseField_label
+ thisNoutput => damage_phaseField_Noutput
+ thisOutput => damage_phaseField_output
+ thisSize => damage_phaseField_sizePostResult
case default
knownDamage = .false.
end select
@@ -509,9 +519,12 @@ function constitutive_damagedC(ipc,ip,el)
use material, only: &
material_phase, &
LOCAL_DAMAGE_isoBrittle_ID, &
+ LOCAL_DAMAGE_phaseField_ID, &
phase_damage
use damage_isoBrittle, only: &
damage_isoBrittle_getDamagedC66
+ use damage_phaseField, only: &
+ damage_phaseField_getDamagedC66
implicit none
real(pReal), dimension(6,6) :: constitutive_damagedC
@@ -524,6 +537,9 @@ function constitutive_damagedC(ipc,ip,el)
case (LOCAL_DAMAGE_isoBrittle_ID)
constitutive_damagedC = damage_isoBrittle_getDamagedC66(constitutive_homogenizedC(ipc,ip,el), &
ipc,ip,el)
+ case (LOCAL_DAMAGE_phaseField_ID)
+ constitutive_damagedC = damage_phaseField_getDamagedC66(constitutive_homogenizedC(ipc,ip,el), &
+ ipc,ip,el)
case default
constitutive_damagedC = constitutive_homogenizedC(ipc,ip,el)
@@ -547,7 +563,8 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el)
PLASTICITY_nonlocal_ID, &
LOCAL_DAMAGE_isoBrittle_ID, &
LOCAL_DAMAGE_isoDuctile_ID, &
- LOCAL_DAMAGE_gurson_ID
+ LOCAL_DAMAGE_gurson_ID, &
+ LOCAL_DAMAGE_phaseField_ID
use constitutive_titanmod, only: &
constitutive_titanmod_microstructure
@@ -559,6 +576,10 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el)
constitutive_dislokmc_microstructure
use damage_gurson, only: &
damage_gurson_microstructure
+ use damage_gurson, only: &
+ damage_gurson_microstructure
+ use damage_phaseField, only: &
+ damage_phaseField_microstructure
implicit none
integer(pInt), intent(in) :: &
@@ -587,6 +608,10 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el)
select case (phase_damage(material_phase(ipc,ip,el)))
case (LOCAL_DAMAGE_gurson_ID)
call damage_gurson_microstructure(ipc, ip, el)
+ case (LOCAL_DAMAGE_phaseField_ID)
+ call damage_phaseField_microstructure(constitutive_homogenizedC(ipc,ip,el), Fe, &
+ constitutive_getVacancyConcentration(ipc, ip, el), &
+ ipc, ip, el)
end select
@@ -1037,6 +1062,7 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
LOCAL_DAMAGE_anisoDuctile_ID, &
LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_DAMAGE_gurson_ID, &
+ LOCAL_DAMAGE_phaseField_ID, &
LOCAL_VACANCY_generation_ID
use constitutive_j2, only: &
constitutive_j2_dotState
@@ -1060,6 +1086,8 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
damage_anisoDuctile_dotState
use damage_gurson, only: &
damage_gurson_dotState
+ use damage_phaseField, only: &
+ damage_phaseField_dotState
use vacancy_generation, only: &
vacancy_generation_dotState
@@ -1121,6 +1149,8 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
call damage_anisoDuctile_dotState(nSlip, accumulatedSlip, ipc, ip, el)
case (LOCAL_DAMAGE_gurson_ID)
call damage_gurson_dotState(Tstar_v, Lp, ipc, ip, el)
+ case (LOCAL_DAMAGE_phaseField_ID)
+ call damage_phaseField_dotState(constitutive_getVacancyConcentration(ipc, ip, el), ipc, ip, el)
end select
select case (phase_vacancy(material_phase(ipc,ip,el)))
@@ -1216,6 +1246,7 @@ function constitutive_getLocalDamage(ipc, ip, el)
LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_DAMAGE_anisoDuctile_ID, &
LOCAL_DAMAGE_gurson_ID, &
+ LOCAL_DAMAGE_phaseField_ID, &
phase_damage
use damage_isoBrittle, only: &
damage_isoBrittle_getLocalDamage
@@ -1227,6 +1258,8 @@ function constitutive_getLocalDamage(ipc, ip, el)
damage_anisoDuctile_getLocalDamage
use damage_gurson, only: &
damage_gurson_getLocalDamage
+ use damage_phaseField, only: &
+ damage_phaseField_getLocalDamage
implicit none
integer(pInt), intent(in) :: &
@@ -1248,13 +1281,15 @@ function constitutive_getLocalDamage(ipc, ip, el)
case (LOCAL_DAMAGE_anisoBrittle_ID)
constitutive_getLocalDamage = damage_anisoBrittle_getLocalDamage(ipc, ip, el)
- case (LOCAL_DAMAGE_anisoDuctile_ID)
-
+ case (LOCAL_DAMAGE_anisoDuctile_ID)
constitutive_getLocalDamage = damage_anisoDuctile_getLocalDamage(ipc, ip, el)
case (LOCAL_DAMAGE_gurson_ID)
constitutive_getLocalDamage = damage_gurson_getLocalDamage(ipc, ip, el)
+ case (LOCAL_DAMAGE_phaseField_ID)
+ constitutive_getLocalDamage = damage_phaseField_getLocalDamage(ipc, ip, el)
+
end select
end function constitutive_getLocalDamage
@@ -1272,6 +1307,7 @@ subroutine constitutive_putLocalDamage(ipc, ip, el, localDamage)
LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_DAMAGE_anisoDuctile_ID, &
LOCAL_DAMAGE_gurson_ID, &
+ LOCAL_DAMAGE_phaseField_ID, &
phase_damage
use damage_isoBrittle, only: &
damage_isoBrittle_putLocalDamage
@@ -1283,6 +1319,8 @@ subroutine constitutive_putLocalDamage(ipc, ip, el, localDamage)
damage_anisoDuctile_putLocalDamage
use damage_gurson, only: &
damage_gurson_putLocalDamage
+ use damage_phaseField, only: &
+ damage_phaseField_putLocalDamage
implicit none
integer(pInt), intent(in) :: &
@@ -1308,6 +1346,9 @@ subroutine constitutive_putLocalDamage(ipc, ip, el, localDamage)
case (LOCAL_DAMAGE_gurson_ID)
call damage_gurson_putLocalDamage(ipc, ip, el, localDamage)
+ case (LOCAL_DAMAGE_phaseField_ID)
+ call damage_phaseField_putLocalDamage(ipc, ip, el, localDamage)
+
end select
end subroutine constitutive_putLocalDamage
@@ -1326,6 +1367,7 @@ function constitutive_getDamage(ipc, ip, el)
LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_DAMAGE_anisoDuctile_ID, &
LOCAL_DAMAGE_gurson_ID, &
+ LOCAL_DAMAGE_phaseField_ID, &
phase_damage
use damage_isoBrittle, only: &
damage_isoBrittle_getDamage
@@ -1337,6 +1379,8 @@ function constitutive_getDamage(ipc, ip, el)
damage_anisoDuctile_getDamage
use damage_gurson, only: &
damage_gurson_getDamage
+ use damage_phaseField, only: &
+ damage_phaseField_getDamage
implicit none
integer(pInt), intent(in) :: &
@@ -1364,6 +1408,9 @@ function constitutive_getDamage(ipc, ip, el)
case (LOCAL_DAMAGE_gurson_ID)
constitutive_getDamage = damage_gurson_getDamage(ipc, ip, el)
+ case (LOCAL_DAMAGE_phaseField_ID)
+ constitutive_getDamage = damage_phaseField_getDamage(ipc, ip, el)
+
end select
end function constitutive_getDamage
@@ -1424,14 +1471,13 @@ function constitutive_getDamageDiffusion33(ipc, ip, el)
lattice_DamageDiffusion33
use material, only: &
material_phase, &
- LOCAL_DAMAGE_none_ID, &
+ phase_damage, &
LOCAL_DAMAGE_isoBrittle_ID, &
- LOCAL_DAMAGE_isoDuctile_ID, &
- LOCAL_DAMAGE_anisoBrittle_ID, &
- LOCAL_DAMAGE_gurson_ID, &
- phase_damage
+ LOCAL_DAMAGE_phaseField_ID
use damage_isoBrittle, only: &
damage_isoBrittle_getDamageDiffusion33
+ use damage_phaseField, only: &
+ damage_phaseField_getDamageDiffusion33
implicit none
integer(pInt), intent(in) :: &
@@ -1445,6 +1491,8 @@ function constitutive_getDamageDiffusion33(ipc, ip, el)
select case(phase_damage(material_phase(ipc,ip,el)))
case (LOCAL_DAMAGE_isoBrittle_ID)
constitutive_getDamageDiffusion33 = damage_isoBrittle_getDamageDiffusion33(ipc, ip, el)
+ case (LOCAL_DAMAGE_phaseField_ID)
+ constitutive_getDamageDiffusion33 = damage_phaseField_getDamageDiffusion33(ipc, ip, el)
end select
@@ -1559,7 +1607,7 @@ function constitutive_getLocalVacancyConcentration(ipc, ip, el)
LOCAL_VACANCY_generation_ID, &
phase_vacancy
use vacancy_generation, only: &
- vacancy_generation_getConcentration
+ vacancy_generation_getLocalConcentration
use lattice, only: &
lattice_equilibriumVacancyConcentration
@@ -1576,7 +1624,7 @@ function constitutive_getLocalVacancyConcentration(ipc, ip, el)
lattice_equilibriumVacancyConcentration(material_phase(ipc,ip,el))
case (LOCAL_VACANCY_generation_ID)
- constitutive_getLocalVacancyConcentration = vacancy_generation_getConcentration(ipc, ip, el)
+ constitutive_getLocalVacancyConcentration = vacancy_generation_getLocalConcentration(ipc, ip, el)
end select
end function constitutive_getLocalVacancyConcentration
@@ -1592,7 +1640,7 @@ subroutine constitutive_putLocalVacancyConcentration(ipc, ip, el, localVacancyCo
LOCAL_VACANCY_generation_ID, &
phase_vacancy
use vacancy_generation, only: &
- vacancy_generation_putConcentration
+ vacancy_generation_putLocalConcentration
implicit none
integer(pInt), intent(in) :: &
@@ -1604,7 +1652,7 @@ subroutine constitutive_putLocalVacancyConcentration(ipc, ip, el, localVacancyCo
select case (phase_vacancy(material_phase(ipc,ip,el)))
case (LOCAL_VACANCY_generation_ID)
- call vacancy_generation_putConcentration(ipc, ip, el, localVacancyConcentration)
+ call vacancy_generation_putLocalConcentration(ipc, ip, el, localVacancyConcentration)
end select
@@ -1617,13 +1665,12 @@ function constitutive_getVacancyConcentration(ipc, ip, el)
use prec, only: &
pReal
use material, only: &
- mappingHomogenization, &
material_phase, &
- fieldVacancy, &
- field_vacancy_type, &
- FIELD_VACANCY_local_ID, &
- FIELD_VACANCY_nonlocal_ID, &
- material_homog
+ LOCAL_VACANCY_constant_ID, &
+ LOCAL_VACANCY_generation_ID, &
+ phase_vacancy
+ use vacancy_generation, only: &
+ vacancy_generation_getConcentration
use lattice, only: &
lattice_equilibriumVacancyConcentration
implicit none
@@ -1634,14 +1681,14 @@ function constitutive_getVacancyConcentration(ipc, ip, el)
el !< element number
real(pReal) :: constitutive_getVacancyConcentration
- select case(field_vacancy_type(material_homog(ip,el)))
- case (FIELD_VACANCY_local_ID)
- constitutive_getVacancyConcentration = constitutive_getLocalVacancyConcentration(ipc, ip, el)
-
- case (FIELD_VACANCY_nonlocal_ID)
- constitutive_getVacancyConcentration = fieldVacancy(material_homog(ip,el))% &
- field(1,mappingHomogenization(1,ip,el)) ! Taylor type
- end select
+ select case (phase_vacancy(material_phase(ipc,ip,el)))
+ case (LOCAL_VACANCY_constant_ID)
+ constitutive_getVacancyConcentration = &
+ lattice_equilibriumVacancyConcentration(material_phase(ipc,ip,el))
+
+ case (LOCAL_VACANCY_generation_ID)
+ constitutive_getVacancyConcentration = vacancy_generation_getConcentration(ipc, ip, el)
+ end select
end function constitutive_getVacancyConcentration
@@ -1667,22 +1714,53 @@ function constitutive_getVacancyDiffusion33(ipc, ip, el)
el !< element number
real(pReal), dimension(3,3) :: &
constitutive_getVacancyDiffusion33
+
+ select case(phase_vacancy(material_phase(ipc,ip,el)))
+ case (LOCAL_VACANCY_generation_ID)
+ constitutive_getVacancyDiffusion33 = &
+ vacancy_generation_getVacancyDiffusion33(ipc,ip,el)
+
+ end select
+
+end function constitutive_getVacancyDiffusion33
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns vacancy diffusion tensor
+!--------------------------------------------------------------------------------------------------
+function constitutive_getVacancyMobility33(ipc, ip, el)
+ use prec, only: &
+ pReal
+ use lattice, only: &
+ lattice_VacancyDiffusion33
+ use material, only: &
+ material_phase, &
+ LOCAL_VACANCY_generation_ID, &
+ phase_vacancy
+ use vacancy_generation, only: &
+ vacancy_generation_getVacancyMobility33
+
+ 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
- constitutive_getVacancyDiffusion33 = lattice_VacancyDiffusion33(1:3,1:3,material_phase(ipc,ip,el))
select case(phase_vacancy(material_phase(ipc,ip,el)))
case (LOCAL_VACANCY_generation_ID)
call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,ipc,ip,el)
- constitutive_getVacancyDiffusion33 = &
- vacancy_generation_getVacancyDiffusion33(nSlip,accumulatedSlip,constitutive_getTemperature(ipc,ip,el), &
+ constitutive_getVacancyMobility33 = &
+ vacancy_generation_getVacancyMobility33(nSlip,accumulatedSlip,constitutive_getTemperature(ipc,ip,el), &
ipc,ip,el)
end select
-end function constitutive_getVacancyDiffusion33
+end function constitutive_getVacancyMobility33
!--------------------------------------------------------------------------------------------------
!> @brief returns accumulated slip on each system defined
diff --git a/code/damage_phaseField.f90 b/code/damage_phaseField.f90
new file mode 100644
index 000000000..298acfdae
--- /dev/null
+++ b/code/damage_phaseField.f90
@@ -0,0 +1,520 @@
+!--------------------------------------------------------------------------------------------------
+! $Id: damage_phaseField.f90 3736 2014-11-21 13:12:54Z MPIE\l.sharma $
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine incoprorating isotropic brittle damage
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module damage_phaseField
+ use prec, only: &
+ pReal, &
+ pInt
+
+ implicit none
+ private
+ integer(pInt), dimension(:), allocatable, public, protected :: &
+ damage_phaseField_sizePostResults !< cumulative size of post results
+
+ integer(pInt), dimension(:,:), allocatable, target, public :: &
+ damage_phaseField_sizePostResult !< size of each post result output
+
+ character(len=64), dimension(:,:), allocatable, target, public :: &
+ damage_phaseField_output !< name of each post result output
+
+ integer(pInt), dimension(:), allocatable, target, public :: &
+ damage_phaseField_Noutput !< number of outputs per instance of this damage
+
+ real(pReal), dimension(:), allocatable, private :: &
+ damage_phaseField_aTol, &
+ damage_phaseField_surfaceEnergy, &
+ damage_phaseField_vacancyFormationEnergy
+
+ enum, bind(c)
+ enumerator :: undefined_ID, &
+ local_damage_ID
+ end enum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 ToDo
+
+ integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
+ damage_phaseField_outputID !< ID of each post result output
+
+
+ public :: &
+ damage_phaseField_init, &
+ damage_phaseField_stateInit, &
+ damage_phaseField_aTolState, &
+ damage_phaseField_dotState, &
+ damage_phaseField_microstructure, &
+ damage_phaseField_getDamage, &
+ damage_phaseField_putLocalDamage, &
+ damage_phaseField_getLocalDamage, &
+ damage_phaseField_getDamageDiffusion33, &
+ damage_phaseField_getDamagedC66, &
+ damage_phaseField_postResults
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_init(fileUnit)
+ use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+ debug_level,&
+ debug_constitutive,&
+ debug_levelBasic
+ use mesh, only: &
+ mesh_maxNips, &
+ mesh_NcpElems
+ use IO, only: &
+ IO_read, &
+ IO_lc, &
+ IO_getTag, &
+ IO_isBlank, &
+ IO_stringPos, &
+ IO_stringValue, &
+ IO_floatValue, &
+ IO_intValue, &
+ IO_warning, &
+ IO_error, &
+ IO_timeStamp, &
+ IO_EOF
+ use material, only: &
+ homogenization_maxNgrains, &
+ phase_damage, &
+ phase_damageInstance, &
+ phase_Noutput, &
+ LOCAL_damage_phaseField_label, &
+ LOCAL_damage_phaseField_ID, &
+ material_phase, &
+ damageState, &
+ MATERIAL_partPhase
+ use numerics,only: &
+ worldrank, &
+ numerics_integrator
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
+ integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions
+ integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,o
+ integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: NofMyPhase
+ character(len=65536) :: &
+ tag = '', &
+ line = ''
+
+ mainProcess: if (worldrank == 0) then
+ write(6,'(/,a)') ' <<<+- damage_'//LOCAL_damage_phaseField_label//' init -+>>>'
+ write(6,'(a)') ' $Id: damage_phaseField.f90 3736 2014-11-21 13:12:54Z MPIE\l.sharma $'
+ write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+
+ maxNinstance = int(count(phase_damage == LOCAL_damage_phaseField_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+ write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+
+ allocate(damage_phaseField_sizePostResults(maxNinstance), source=0_pInt)
+ allocate(damage_phaseField_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(damage_phaseField_output(maxval(phase_Noutput),maxNinstance))
+ damage_phaseField_output = ''
+ allocate(damage_phaseField_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
+ 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_aTol(maxNinstance), source=0.0_pReal)
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to
+ line = IO_read(fileUnit)
+ enddo
+
+ parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
+ line = IO_read(fileUnit)
+ if (IO_isBlank(line)) cycle ! skip empty lines
+ if (IO_getTag(line,'<','>') /= '') then ! stop at next part
+ line = IO_read(fileUnit, .true.) ! reset IO_read
+ exit
+ endif
+ if (IO_getTag(line,'[',']') /= '') then ! next phase section
+ phase = phase + 1_pInt ! advance phase section counter
+ cycle ! skip to next line
+ endif
+ if (phase > 0_pInt ) then; if (phase_damage(phase) == LOCAL_damage_phaseField_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+ instance = phase_damageInstance(phase) ! which instance of my damage is present phase
+ positions = IO_stringPos(line,MAXNCHUNKS)
+ tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
+ select case(tag)
+ case ('(output)')
+ select case(IO_lc(IO_stringValue(line,positions,2_pInt)))
+ case ('local_damage')
+ damage_phaseField_Noutput(instance) = damage_phaseField_Noutput(instance) + 1_pInt
+ damage_phaseField_outputID(damage_phaseField_Noutput(instance),instance) = local_damage_ID
+ damage_phaseField_output(damage_phaseField_Noutput(instance),instance) = &
+ IO_lc(IO_stringValue(line,positions,2_pInt))
+ end select
+
+ case ('surfaceenergy')
+ damage_phaseField_surfaceEnergy(instance) = IO_floatValue(line,positions,2_pInt)
+
+ case ('vacancyformationenergy')
+ damage_phaseField_vacancyFormationEnergy(instance) = IO_floatValue(line,positions,2_pInt)
+
+ case ('atol_damage')
+ damage_phaseField_aTol(instance) = IO_floatValue(line,positions,2_pInt)
+
+ end select
+ endif; endif
+ enddo parsingFile
+
+
+ sanityChecks: do phase = 1_pInt, size(phase_damage)
+ myPhase: if (phase_damage(phase) == LOCAL_damage_phaseField_ID) then
+ NofMyPhase=count(material_phase==phase)
+ instance = phase_damageInstance(phase)
+! sanity checks
+ if (damage_phaseField_aTol(instance) < 0.0_pReal) &
+ damage_phaseField_aTol(instance) = 1.0e-3_pReal ! default absolute tolerance 1e-3
+ if (damage_phaseField_surfaceEnergy(instance) <= 0.0_pReal) &
+ call IO_error(211_pInt,el=instance,ext_msg='surfaceEnergy ('//LOCAL_damage_phaseField_LABEL//')')
+ endif myPhase
+ enddo sanityChecks
+
+ initializeInstances: do phase = 1_pInt, size(phase_damage)
+ if (phase_damage(phase) == LOCAL_damage_phaseField_ID) then
+ NofMyPhase=count(material_phase==phase)
+ instance = phase_damageInstance(phase)
+!--------------------------------------------------------------------------------------------------
+! Determine size of postResults array
+ outputsLoop: do o = 1_pInt,damage_phaseField_Noutput(instance)
+ select case(damage_phaseField_outputID(o,instance))
+ case(local_damage_ID)
+ mySize = 1_pInt
+ end select
+
+ if (mySize > 0_pInt) then ! any meaningful output found
+ damage_phaseField_sizePostResult(o,instance) = mySize
+ damage_phaseField_sizePostResults(instance) = damage_phaseField_sizePostResults(instance) + mySize
+ endif
+ enddo outputsLoop
+! Determine size of state array
+ sizeDotState = 1_pInt
+ sizeState = 2_pInt
+
+ damageState(phase)%sizeState = sizeState
+ damageState(phase)%sizeDotState = sizeDotState
+ damageState(phase)%sizePostResults = damage_phaseField_sizePostResults(instance)
+ allocate(damageState(phase)%aTolState (sizeState), source=0.0_pReal)
+ allocate(damageState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
+ allocate(damageState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
+ allocate(damageState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
+ allocate(damageState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal)
+ allocate(damageState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
+
+ allocate(damageState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
+ allocate(damageState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal)
+ allocate(damageState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
+ if (any(numerics_integrator == 1_pInt)) then
+ allocate(damageState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
+ allocate(damageState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
+ endif
+ if (any(numerics_integrator == 4_pInt)) &
+ allocate(damageState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
+ if (any(numerics_integrator == 5_pInt)) &
+ allocate(damageState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
+
+ call damage_phaseField_stateInit(phase)
+ call damage_phaseField_aTolState(phase,instance)
+ endif
+
+ enddo initializeInstances
+end subroutine damage_phaseField_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief sets the relevant NEW state values for a given instance of this damage
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_stateInit(phase)
+ use material, only: &
+ damageState
+
+ implicit none
+ integer(pInt), intent(in) :: phase !< number specifying the phase of the damage
+
+ real(pReal), dimension(damageState(phase)%sizeState) :: tempState
+
+ tempState = 1.0_pReal
+ damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:)))
+ damageState(phase)%state0 = damageState(phase)%state
+ damageState(phase)%partionedState0 = damageState(phase)%state
+end subroutine damage_phaseField_stateInit
+
+!--------------------------------------------------------------------------------------------------
+!> @brief sets the relevant state values for a given instance of this damage
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_aTolState(phase,instance)
+ use material, only: &
+ damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ phase, &
+ instance ! number specifying the current instance of the damage
+ real(pReal), dimension(damageState(phase)%sizeState) :: tempTol
+
+ tempTol = damage_phaseField_aTol(instance)
+ damageState(phase)%aTolState = tempTol
+end subroutine damage_phaseField_aTolState
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_dotState(Cv, ipc, ip, el)
+ use material, only: &
+ mappingConstitutive, &
+ damageState
+ use lattice, only: &
+ lattice_DamageMobility
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ ipc, & !< component-ID of integration point
+ ip, & !< integration point
+ el !< element
+ real(pReal), intent(in) :: &
+ Cv
+ integer(pInt) :: &
+ phase, constituent
+
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+
+ damageState(phase)%dotState(1,constituent) = &
+ (damageState(phase)%state(2,constituent)*(1.0_pReal - Cv) - &
+ damageState(phase)%state(1,constituent))/ &
+ lattice_DamageMobility(phase)
+
+end subroutine damage_phaseField_dotState
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_microstructure(C, Fe, Cv, ipc, ip, el)
+ use material, only: &
+ mappingConstitutive, &
+ phase_damageInstance, &
+ damageState
+ use math, only : &
+ math_mul33x33, &
+ math_mul66x6, &
+ math_Mandel33to6, &
+ math_transpose33, &
+ math_I3
+ use lattice, only: &
+ lattice_DamageMobility
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ ipc, & !< component-ID of integration point
+ ip, & !< integration point
+ el !< element
+ real(pReal), intent(in), dimension(3,3) :: &
+ Fe
+ real(pReal), intent(in), dimension(6,6) :: &
+ C
+ real(pReal), intent(in) :: &
+ Cv
+ integer(pInt) :: &
+ phase, constituent, instance
+ real(pReal) :: &
+ strain(6), &
+ stress(6)
+
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ instance = phase_damageInstance(phase)
+
+ strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
+ stress = math_mul66x6(C,strain)
+
+ damageState(phase)%state(2,constituent) = &
+ min(damageState(phase)%state0(2,constituent), &
+ damage_phaseField_surfaceEnergy(instance)/ &
+ (2.0_pReal*(sum(abs(stress*strain)) + Cv*damage_phaseField_vacancyFormationEnergy(instance))))
+
+end subroutine damage_phaseField_microstructure
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns damage
+!--------------------------------------------------------------------------------------------------
+function damage_phaseField_getDamage(ipc, ip, el)
+ use material, only: &
+ material_homog, &
+ mappingHomogenization, &
+ fieldDamage, &
+ field_damage_type, &
+ FIELD_DAMAGE_LOCAL_ID, &
+ FIELD_DAMAGE_NONLOCAL_ID
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ ipc, & !< grain number
+ ip, & !< integration point number
+ el !< element number
+ real(pReal) :: damage_phaseField_getDamage
+
+ select case(field_damage_type(material_homog(ip,el)))
+ case (FIELD_DAMAGE_LOCAL_ID)
+ damage_phaseField_getDamage = damage_phaseField_getLocalDamage(ipc, ip, el)
+
+ case (FIELD_DAMAGE_NONLOCAL_ID)
+ damage_phaseField_getDamage = fieldDamage(material_homog(ip,el))% &
+ field(1,mappingHomogenization(1,ip,el)) ! Taylor type
+
+ end select
+
+end function damage_phaseField_getDamage
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns temperature based on local damage model state layout
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_putLocalDamage(ipc, ip, el, localDamage)
+ use material, only: &
+ mappingConstitutive, &
+ damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ ipc, & !< grain number
+ ip, & !< integration point number
+ el !< element number
+ real(pReal), intent(in) :: localDamage
+
+ damageState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el)) = &
+ localDamage
+
+end subroutine damage_phaseField_putLocalDamage
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns local damage
+!--------------------------------------------------------------------------------------------------
+function damage_phaseField_getLocalDamage(ipc, ip, el)
+ use material, only: &
+ mappingConstitutive, &
+ damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ ipc, & !< grain number
+ ip, & !< integration point number
+ el !< element number
+ real(pReal) :: damage_phaseField_getLocalDamage
+
+ damage_phaseField_getLocalDamage = &
+ damageState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el))
+
+end function damage_phaseField_getLocalDamage
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns brittle damage diffusion tensor
+!--------------------------------------------------------------------------------------------------
+function damage_phaseField_getDamageDiffusion33(ipc, ip, el)
+ use lattice, only: &
+ lattice_DamageDiffusion33
+ use material, only: &
+ mappingConstitutive, &
+ damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ ipc, & !< grain number
+ ip, & !< integration point number
+ el !< element number
+ real(pReal), dimension(3,3) :: &
+ damage_phaseField_getDamageDiffusion33
+ integer(pInt) :: &
+ phase, constituent
+
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ damage_phaseField_getDamageDiffusion33 = &
+ damageState(phase)%state(2,constituent)* &
+ lattice_DamageDiffusion33(1:3,1:3,phase)
+
+end function damage_phaseField_getDamageDiffusion33
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns brittle damaged stiffness tensor
+!--------------------------------------------------------------------------------------------------
+function damage_phaseField_getDamagedC66(C, ipc, ip, el)
+ use material, only: &
+ mappingConstitutive, &
+ damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ ipc, & !< grain number
+ ip, & !< integration point number
+ el !< element number
+ real(pReal), intent(in), dimension(6,6) :: &
+ C
+ real(pReal), dimension(6,6) :: &
+ damage_phaseField_getDamagedC66
+ integer(pInt) :: &
+ phase, constituent
+ real(pReal) :: &
+ damage
+
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ damage = damage_phaseField_getDamage(ipc, ip, el)
+
+ damage_phaseField_getDamagedC66 = &
+ damage*damage*C
+
+end function damage_phaseField_getDamagedC66
+
+!--------------------------------------------------------------------------------------------------
+!> @brief return array of constitutive results
+!--------------------------------------------------------------------------------------------------
+function damage_phaseField_postResults(ipc,ip,el)
+ use material, only: &
+ mappingConstitutive, &
+ phase_damageInstance,&
+ damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ ipc, & !< component-ID of integration point
+ ip, & !< integration point
+ el !< element
+ real(pReal), dimension(damage_phaseField_sizePostResults(phase_damageInstance(mappingConstitutive(2,ipc,ip,el)))) :: &
+ damage_phaseField_postResults
+
+ integer(pInt) :: &
+ instance, phase, constituent, o, c
+
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ instance = phase_damageInstance(phase)
+
+ c = 0_pInt
+ damage_phaseField_postResults = 0.0_pReal
+
+ do o = 1_pInt,damage_phaseField_Noutput(instance)
+ select case(damage_phaseField_outputID(o,instance))
+ case (local_damage_ID)
+ damage_phaseField_postResults(c+1_pInt) = damageState(phase)%state(2,constituent)
+ c = c + 1
+
+ end select
+ enddo
+end function damage_phaseField_postResults
+
+end module damage_phaseField
diff --git a/code/material.f90 b/code/material.f90
index 65ca5d96d..005956d58 100644
--- a/code/material.f90
+++ b/code/material.f90
@@ -33,6 +33,7 @@ module material
LOCAL_DAMAGE_anisoBrittle_LABEL= 'anisobrittle', &
LOCAL_DAMAGE_anisoDuctile_LABEL= 'anisoductile', &
LOCAL_DAMAGE_gurson_LABEL = 'gurson', &
+ LOCAL_DAMAGE_phaseField_LABEL = 'phasefield', &
LOCAL_THERMAL_isothermal_label = 'isothermal', &
LOCAL_THERMAL_adiabatic_label = 'adiabatic', &
LOCAL_VACANCY_constant_label = 'constant', &
@@ -69,7 +70,8 @@ module material
LOCAL_DAMAGE_isoDuctile_ID, &
LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_DAMAGE_anisoDuctile_ID, &
- LOCAL_DAMAGE_gurson_ID
+ LOCAL_DAMAGE_gurson_ID, &
+ LOCAL_DAMAGE_phaseField_ID
end enum
enum, bind(c)
enumerator :: LOCAL_THERMAL_isothermal_ID, &
@@ -244,6 +246,7 @@ module material
LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_DAMAGE_anisoDuctile_ID, &
LOCAL_DAMAGE_gurson_ID, &
+ LOCAL_DAMAGE_phaseField_ID, &
LOCAL_THERMAL_isothermal_ID, &
LOCAL_THERMAL_adiabatic_ID, &
LOCAL_VACANCY_constant_ID, &
@@ -859,6 +862,8 @@ subroutine material_parsePhase(fileUnit,myPart)
phase_damage(section) = LOCAL_DAMAGE_anisoDuctile_ID
case (LOCAL_DAMAGE_gurson_label)
phase_damage(section) = LOCAL_DAMAGE_gurson_ID
+ case (LOCAL_DAMAGE_phaseField_label)
+ phase_damage(section) = LOCAL_DAMAGE_phaseField_ID
case default
call IO_error(200_pInt,ext_msg=trim(IO_stringValue(line,positions,2_pInt)))
end select
diff --git a/code/vacancy_generation.f90 b/code/vacancy_generation.f90
index 33ee42f13..e132050d3 100644
--- a/code/vacancy_generation.f90
+++ b/code/vacancy_generation.f90
@@ -28,21 +28,11 @@ module vacancy_generation
vacancy_generation_aTol, &
vacancy_generation_freq, &
vacancy_generation_formationEnergy, &
- vacancy_generation_diffusionEnergy, &
- vacancy_generation_diffusionCoeff0, & !< the temperature-independent pre-exponential of diffusion coefficient D_0
- vacancy_generation_stressCoeff, &
- vacancy_generation_jogHeight, & !< the height of jogs in Burgers vectors
- vacancy_generation_jogSeparation, & !< the jog seperation
- vacancy_generation_nLatticeSites, & !< the number of lattice sites per unit volume
- vacancy_generation_burgersVec, & !< the Burgers vector
- vacancy_generation_dislocationCoeff, &
- vacancy_generation_equilibConcentration !< the equilibrium concentration of vacancy
-
- real(pReal), dimension(:), allocatable, public :: &
- pore_nucleation_surfaceEnergy, & !< surface energy of metal which controls the necleation of pores
- pore_nucleation_atomVolume, & !< the volume of atom
- pore_nucleation_shellThickness, & !< the thickness of spherical shell surrounding the pore
- pore_nucleation_concentrationCoeff0 !< the pre-exponential of equilibrium concentration of critical pore
+ vacancy_generation_migrationEnergy, &
+ vacancy_generation_diffusionCoeff0, & !< the temperature-independent diffusion coefficient D_0
+ vacancy_generation_atomicVol, &
+ vacancy_generation_surfaceEnergy, &
+ vacancy_generation_plasticityCoeff
real(pReal), parameter, private :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
@@ -60,9 +50,11 @@ module vacancy_generation
vacancy_generation_stateInit, &
vacancy_generation_aTolState, &
vacancy_generation_dotState, &
+ vacancy_generation_getLocalConcentration, &
+ vacancy_generation_putLocalConcentration, &
vacancy_generation_getConcentration, &
- vacancy_generation_putConcentration, &
vacancy_generation_getVacancyDiffusion33, &
+ vacancy_generation_getVacancyMobility33, &
vacancy_generation_postResults
contains
@@ -141,21 +133,10 @@ 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_diffusionEnergy(maxNinstance), source=0.0_pReal)
- allocate(vacancy_generation_stressCoeff(maxNinstance), source=0.0_pReal)
- allocate(vacancy_generation_jogHeight(maxNinstance), source=0.0_pReal)
- allocate(vacancy_generation_jogSeparation(maxNinstance), source=0.0_pReal)
- allocate(vacancy_generation_nLatticeSites(maxNinstance), source=0.0_pReal)
- allocate(vacancy_generation_burgersVec(maxNinstance), source=0.0_pReal)
- allocate(vacancy_generation_diffusionCoeff0(maxNinstance), source=0.0_pReal)
- allocate(vacancy_generation_equilibConcentration(maxNinstance), source=0.0_pReal)
-
- allocate(vacancy_generation_dislocationCoeff(maxNinstance), source=0.0_pReal)
-
- allocate(pore_nucleation_surfaceEnergy(maxNinstance), source=0.0_pReal)
- allocate(pore_nucleation_atomVolume(maxNinstance), source=0.0_pReal)
- allocate(pore_nucleation_shellThickness(maxNinstance), source=0.0_pReal)
- allocate(pore_nucleation_concentrationCoeff0(maxNinstance), source=0.0_pReal)
+ allocate(vacancy_generation_migrationEnergy(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)
rewind(fileUnit)
phase = 0_pInt
@@ -190,50 +171,29 @@ subroutine vacancy_generation_init(fileUnit)
IO_lc(IO_stringValue(line,positions,2_pInt))
end select
- case ('atol_vacancygeneration')
+ case ('atolvacancygeneration')
vacancy_generation_aTol(instance) = IO_floatValue(line,positions,2_pInt)
- case ('vacancy_frequency')
+ case ('debyefrequency')
vacancy_generation_freq(instance) = IO_floatValue(line,positions,2_pInt)
- case ('vacancy_formationenergy')
+ case ('vacancyformationenergy')
vacancy_generation_formationEnergy(instance) = IO_floatValue(line,positions,2_pInt)
- case ('vacancy_equilibconcentration')
- vacancy_generation_equilibConcentration(instance) = IO_floatValue(line,positions,2_pInt)
+ case ('vacancymigrationenergy')
+ vacancy_generation_migrationEnergy(instance) = IO_floatValue(line,positions,2_pInt)
- case ('vacancy_diffusionenergy')
- vacancy_generation_diffusionEnergy(instance) = IO_floatValue(line,positions,2_pInt)
-
- case ('vacancy_diffusioncoeff0')
+ case ('vacancydiffusioncoeff0')
vacancy_generation_diffusionCoeff0(instance) = IO_floatValue(line,positions,2_pInt)
- case ('vacancy_stresscoeff')
- vacancy_generation_stressCoeff(instance) = IO_floatValue(line,positions,2_pInt)
+ case ('atomicvolume')
+ vacancy_generation_atomicVol(instance) = IO_floatValue(line,positions,2_pInt)
- case ('vacancy_jogheight')
- vacancy_generation_jogHeight(instance) = IO_floatValue(line,positions,2_pInt)
+ case ('surfaceenergy')
+ vacancy_generation_surfaceEnergy(instance) = IO_floatValue(line,positions,2_pInt)
- case ('vacancy_jogseparation')
- vacancy_generation_jogSeparation(instance) = IO_floatValue(line,positions,2_pInt)
-
- case ('vacancy_nlatticesites')
- vacancy_generation_nLatticeSites(instance) = IO_floatValue(line,positions,2_pInt)
-
- case ('vacancy_burgersvec')
- vacancy_generation_burgersVec(instance) = IO_floatValue(line,positions,2_pInt)
-
- case ('pore_surfacefnergy')
- pore_nucleation_surfaceEnergy(instance) = IO_floatValue(line,positions,2_pInt)
-
- case ('pore_atomvolume')
- pore_nucleation_atomVolume(instance) = IO_floatValue(line,positions,2_pInt)
-
- case ('pore_shellthickness')
- pore_nucleation_shellThickness(instance) = IO_floatValue(line,positions,2_pInt)
-
- case ('pore_concentrationcoeff0')
- pore_nucleation_concentrationCoeff0(instance) = IO_floatValue(line,positions,2_pInt)
+ case ('vacancyplasticitycoeff')
+ vacancy_generation_plasticityCoeff(instance) = IO_floatValue(line,positions,2_pInt)
end select
endif; endif
@@ -244,14 +204,6 @@ subroutine vacancy_generation_init(fileUnit)
NofMyPhase=count(material_phase==phase)
instance = phase_vacancyInstance(phase)
-!--------------------------------------------------------------------------------------------------
-! Calculate the coefficient for dislocation motion induced vacancy generation
- vacancy_generation_dislocationCoeff(instance) = vacancy_generation_jogHeight(instance)/ &
- vacancy_generation_jogSeparation(instance)/ &
- vacancy_generation_nLatticeSites(instance)/ &
- vacancy_generation_burgersVec(instance)/ &
- vacancy_generation_burgersVec(instance)
-
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,vacancy_generation_Noutput(instance)
@@ -337,17 +289,13 @@ end subroutine vacancy_generation_aTolState
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine vacancy_generation_dotState(nSlip, accumulatedSlip, Tstar_v, Temperature, ipc, ip, el)
- use lattice, only: &
- lattice_massDensity, &
- lattice_specificHeat
use material, only: &
mappingConstitutive, &
phase_vacancyInstance, &
vacancyState
use math, only: &
math_Mandel6to33, &
- math_trace33, &
- pi
+ math_trace33
implicit none
integer(pInt), intent(in) :: &
@@ -362,70 +310,29 @@ subroutine vacancy_generation_dotState(nSlip, accumulatedSlip, Tstar_v, Temperat
real(pReal), intent(in) :: &
Temperature !< 2nd Piola Kirchhoff stress tensor (Mandel)
real(pReal) :: &
- pressure !< 2nd Piola Kirchhoff stress tensor (Mandel)
+ pressure, &
+ energyBarrier
integer(pInt) :: &
instance, phase, constituent
- real(pReal) :: &
- vacancyConcentration, & !< current vacancy concentration
- vacancyDiffusion, & !< the diffusion coefficient D_v
- poleZeldovichCoeff, & !< Zeldovich factor of pore nucleation
- vacancyAbsorpRateCoeff, & !< vacancy absorption rate
- chemicalPotential, & !< the chemical potential due to vacancy concentration
- criticalRadius, & !< the critical pore radius
- Gibbs4Pore, & !< the Gibbs free energy for generating a critical pore
- equilibPoreConcentration, & !< the equilibrium pore concentration
- nucleationRatePore, & !< the nucleation rate of pore
- ratioCvCve !< the ratio of Cv with respect to Cve
- real(pReal) :: &
- threshold4ratioCvCve = 2.0_pReal !< the threshold value for Cv/Cve
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) - &
+ pressure*vacancy_generation_atomicVol(instance) - &
+ sum(accumulatedSlip)*vacancy_generation_plasticityCoeff(instance)
-!--------------------------------------------------------------------------------------------------
- vacancyConcentration = vacancy_generation_getConcentration(ipc, ip, el)
- ratioCvCve = vacancyConcentration/vacancy_generation_equilibConcentration(instance)
-
- if(ratioCvCve < threshold4ratioCvCve) then
- nucleationRatePore = 0.0_pReal
- else
-! Calculate nucleation rate of pore
- vacancyDiffusion = vacancy_generation_diffusionCoeff0(instance)* &
- exp( -vacancy_generation_diffusionEnergy(instance)/(kB*temperature) )
- chemicalPotential = kB*Temperature * log(vacancyConcentration/ &
- vacancy_generation_equilibConcentration(instance))
- criticalRadius = 2_pReal/chemicalPotential* &
- pore_nucleation_surfaceEnergy(instance) * pore_nucleation_atomVolume(instance)
- Gibbs4Pore = 4_pReal/3_pReal * pi * pore_nucleation_surfaceEnergy(instance)* &
- criticalRadius * criticalRadius
- equilibPoreConcentration = pore_nucleation_concentrationCoeff0(instance)* &
- exp( -Gibbs4Pore/(kB*temperature) )
-
- vacancyAbsorpRateCoeff = 2_pReal/pore_nucleation_shellThickness(instance) * &
- vacancyDiffusion * vacancyConcentration
- poleZeldovichCoeff = pore_nucleation_atomVolume(instance)* &
- sqrt( pore_nucleation_surfaceEnergy(instance)/(kB*temperature) )
- nucleationRatePore = poleZeldovichCoeff * vacancyAbsorpRateCoeff* equilibPoreConcentration
- endif
-
-!--------------------------------------------------------------------------------------------------
-! the net generating rate vacancy
vacancyState(phase)%dotState(1,constituent) = &
vacancy_generation_freq(instance)* &
- exp(-(vacancy_generation_formationEnergy(instance) - vacancy_generation_stressCoeff(instance)*pressure)/ &
- (kB*Temperature)) + &
- sum(accumulatedSlip) * vacancy_generation_dislocationCoeff(instance)- & !< Induced by dislocation motion
- nucleationRatePore * (4_pReal/3_pReal * pi * criticalRadius**3_pReal)/ & !< Reduced by the formation of pore
- pore_nucleation_atomVolume(instance)
+ exp(-energyBarrier/(kB*Temperature))
end subroutine vacancy_generation_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns vacancy concentration based on state layout
!--------------------------------------------------------------------------------------------------
-function vacancy_generation_getConcentration(ipc, ip, el)
+function vacancy_generation_getLocalConcentration(ipc, ip, el)
use material, only: &
mappingConstitutive, &
vacancyState
@@ -435,17 +342,17 @@ function vacancy_generation_getConcentration(ipc, ip, el)
ipc, & !< grain number
ip, & !< integration point number
el !< element number
- real(pReal) :: vacancy_generation_getConcentration
+ real(pReal) :: vacancy_generation_getLocalConcentration
- vacancy_generation_getConcentration = &
+ vacancy_generation_getLocalConcentration = &
vacancyState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el))
-end function vacancy_generation_getConcentration
+end function vacancy_generation_getLocalConcentration
!--------------------------------------------------------------------------------------------------
!> @brief returns temperature based on local damage model state layout
!--------------------------------------------------------------------------------------------------
-subroutine vacancy_generation_putConcentration(ipc, ip, el, localVacancyConcentration)
+subroutine vacancy_generation_putLocalConcentration(ipc, ip, el, localVacancyConcentration)
use material, only: &
mappingConstitutive, &
vacancyState
@@ -461,18 +368,68 @@ subroutine vacancy_generation_putConcentration(ipc, ip, el, localVacancyConcentr
vacancyState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el))= &
localVacancyConcentration
-end subroutine vacancy_generation_putConcentration
+end subroutine vacancy_generation_putLocalConcentration
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns vacancy concentration based on state layout
+!--------------------------------------------------------------------------------------------------
+function vacancy_generation_getConcentration(ipc, ip, el)
+ use material, only: &
+ mappingHomogenization, &
+ material_phase, &
+ fieldVacancy, &
+ field_vacancy_type, &
+ FIELD_VACANCY_local_ID, &
+ FIELD_VACANCY_nonlocal_ID, &
+ material_homog
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ ipc, & !< grain number
+ ip, & !< integration point number
+ el !< element number
+ real(pReal) :: vacancy_generation_getConcentration
+
+ select case(field_vacancy_type(material_homog(ip,el)))
+ case (FIELD_VACANCY_local_ID)
+ vacancy_generation_getConcentration = vacancy_generation_getLocalConcentration(ipc, ip, el)
+
+ case (FIELD_VACANCY_nonlocal_ID)
+ vacancy_generation_getConcentration = fieldVacancy(material_homog(ip,el))% &
+ field(1,mappingHomogenization(1,ip,el)) ! Taylor type
+ end select
+
+end function vacancy_generation_getConcentration
!--------------------------------------------------------------------------------------------------
!> @brief returns generation vacancy diffusion tensor
!--------------------------------------------------------------------------------------------------
-function vacancy_generation_getVacancyDiffusion33(nSlip,accumulatedSlip,temperature,ipc,ip,el)
+function vacancy_generation_getVacancyDiffusion33(ipc,ip,el)
use lattice, only: &
lattice_VacancyDiffusion33
+ use material, only: &
+ mappingConstitutive
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ ipc, & !< grain number
+ ip, & !< integration point number
+ el !< element number
+ real(pReal), dimension(3,3) :: &
+ vacancy_generation_getVacancyDiffusion33
+
+ vacancy_generation_getVacancyDiffusion33 = &
+ lattice_VacancyDiffusion33(1:3,1:3,mappingConstitutive(2,ipc,ip,el))
+
+end function vacancy_generation_getVacancyDiffusion33
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns generation vacancy mobility tensor
+!--------------------------------------------------------------------------------------------------
+function vacancy_generation_getVacancyMobility33(nSlip,accumulatedSlip,temperature,ipc,ip,el)
use material, only: &
mappingConstitutive, &
- phase_vacancyInstance, &
- vacancyState
+ phase_vacancyInstance
implicit none
integer(pInt), intent(in) :: &
@@ -481,7 +438,7 @@ function vacancy_generation_getVacancyDiffusion33(nSlip,accumulatedSlip,temperat
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
- vacancy_generation_getVacancyDiffusion33
+ vacancy_generation_getVacancyMobility33
real(pReal), dimension(nSlip) :: &
accumulatedSlip
real(pReal) :: &
@@ -493,11 +450,13 @@ function vacancy_generation_getVacancyDiffusion33(nSlip,accumulatedSlip,temperat
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_vacancyInstance(phase)
- vacancy_generation_getVacancyDiffusion33 = &
- lattice_VacancyDiffusion33(1:3,1:3,phase)* &
- exp(-vacancy_generation_diffusionEnergy(instance)/(kB*temperature))
+ vacancy_generation_getVacancyMobility33 = &
+ vacancy_generation_surfaceEnergy(instance)* &
+ vacancy_generation_diffusionCoeff0(instance)* &
+ exp(-vacancy_generation_migrationEnergy(instance)/(kB*temperature))/ &
+ (kB*temperature)
-end function vacancy_generation_getVacancyDiffusion33
+end function vacancy_generation_getVacancyMobility33
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results