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