From dd16851ab78b52ac924016633f101e384a6ab738 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Fri, 10 Oct 2014 20:55:09 +0000 Subject: [PATCH] implemented possibly diffusive vacancy physics. to be coupled with micro void nucleation and ductile damage --- code/DAMASK_abaqus_exp.f | 2 + code/DAMASK_abaqus_std.f | 2 + code/DAMASK_marc.f90 | 2 + code/Makefile | 16 +- code/constitutive.f90 | 197 ++++++++++++++++-- code/constitutive_nonlocal.f90 | 1 - code/homogenization.f90 | 175 ++++++++++++++-- code/lattice.f90 | 32 ++- code/material.f90 | 86 +++++++- code/thermal_adiabatic.f90 | 5 +- code/vacancy_constant.f90 | 103 ++++++++++ code/vacancy_generation.f90 | 365 +++++++++++++++++++++++++++++++++ 12 files changed, 931 insertions(+), 55 deletions(-) create mode 100644 code/vacancy_constant.f90 create mode 100644 code/vacancy_generation.f90 diff --git a/code/DAMASK_abaqus_exp.f b/code/DAMASK_abaqus_exp.f index 37c9e5c17..5887ad2a4 100644 --- a/code/DAMASK_abaqus_exp.f +++ b/code/DAMASK_abaqus_exp.f @@ -88,6 +88,8 @@ end module DAMASK_interface #include "damage_gurson.f90" #include "thermal_isothermal.f90" #include "thermal_adiabatic.f90" +#include "vacancy_constant.f90" +#include "vacancy_generation.f90" #include "constitutive_none.f90" #include "constitutive_j2.f90" #include "constitutive_phenopowerlaw.f90" diff --git a/code/DAMASK_abaqus_std.f b/code/DAMASK_abaqus_std.f index a87e95f2f..e11144f0c 100644 --- a/code/DAMASK_abaqus_std.f +++ b/code/DAMASK_abaqus_std.f @@ -88,6 +88,8 @@ end module DAMASK_interface #include "damage_gurson.f90" #include "thermal_isothermal.f90" #include "thermal_adiabatic.f90" +#include "vacancy_constant.f90" +#include "vacancy_generation.f90" #include "constitutive_none.f90" #include "constitutive_j2.f90" #include "constitutive_phenopowerlaw.f90" diff --git a/code/DAMASK_marc.f90 b/code/DAMASK_marc.f90 index 5924547b8..549c6ffcd 100644 --- a/code/DAMASK_marc.f90 +++ b/code/DAMASK_marc.f90 @@ -117,6 +117,8 @@ end module DAMASK_interface #include "damage_gurson.f90" #include "thermal_isothermal.f90" #include "thermal_adiabatic.f90" +#include "vacancy_constant.f90" +#include "vacancy_generation.f90" #include "constitutive_none.f90" #include "constitutive_j2.f90" #include "constitutive_phenopowerlaw.f90" diff --git a/code/Makefile b/code/Makefile index 5ee266fd9..4104c9554 100644 --- a/code/Makefile +++ b/code/Makefile @@ -331,6 +331,9 @@ DAMAGE_FILES = \ THERMAL_FILES = \ thermal_isothermal.o thermal_adiabatic.o +VACANCY_FILES = \ + vacancy_constant.o vacancy_generation.o + CONSTITUTIVE_FILES = \ constitutive_dislotwin.o constitutive_dislokmc.o constitutive_j2.o constitutive_phenopowerlaw.o \ constitutive_titanmod.o constitutive_nonlocal.o constitutive_none.o constitutive.o @@ -350,7 +353,7 @@ DAMASK_spectral.exe: INTERFACENAME := DAMASK_spectral_interface.f90 SPECTRAL_FILES = prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o \ FEsolving.o mesh.o material.o lattice.o \ - $(DAMAGE_FILES) $(THERMAL_FILES) $(CONSTITUTIVE_FILES) \ + $(DAMAGE_FILES) $(THERMAL_FILES) $(VACANCY_FILES) $(CONSTITUTIVE_FILES) \ crystallite.o $(HOMOGENIZATION_FILES) CPFEM.o \ DAMASK_spectral_utilities.o DAMASK_spectral_solverBasic.o \ DAMASK_spectral_solverAL.o DAMASK_spectral_solverBasicPETSc.o DAMASK_spectral_solverPolarisation.o @@ -393,7 +396,7 @@ DAMASK_FEM.exe: INCLUDE_DIRS += -I./ FEM_FILES = prec.o DAMASK_interface.o FEZoo.o IO.o libs.o numerics.o debug.o math.o \ FEsolving.o mesh.o material.o lattice.o \ - $(DAMAGE_FILES) $(THERMAL_FILES) $(CONSTITUTIVE_FILES) \ + $(DAMAGE_FILES) $(THERMAL_FILES) $(VACANCY_FILES) $(CONSTITUTIVE_FILES) \ crystallite.o $(HOMOGENIZATION_FILES) CPFEM.o \ FEM_utilities.o FEM_mech.o FEM_thermal.o FEM_damage.o @@ -451,7 +454,8 @@ constitutive.o: constitutive.f90 \ constitutive_j2.o \ constitutive_none.o \ $(DAMAGE_FILES) \ - $(THERMAL_FILES) + $(THERMAL_FILES) \ + $(VACANCY_FILES) constitutive_nonlocal.o: constitutive_nonlocal.f90 \ lattice.o @@ -492,6 +496,12 @@ thermal_isothermal.o: thermal_isothermal.f90 \ thermal_adiabatic.o: thermal_adiabatic.f90 \ lattice.o +vacancy_constant.o: vacancy_constant.f90 \ + lattice.o + +vacancy_generation.o: vacancy_generation.f90 \ + lattice.o + lattice.o: lattice.f90 \ material.o diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 16021b68d..76fb5f795 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -17,7 +17,9 @@ module constitutive constitutive_damage_maxSizePostResults, & constitutive_damage_maxSizeDotState, & constitutive_thermal_maxSizePostResults, & - constitutive_thermal_maxSizeDotState + constitutive_thermal_maxSizeDotState, & + constitutive_vacancy_maxSizePostResults, & + constitutive_vacancy_maxSizeDotState public :: & constitutive_init, & @@ -33,6 +35,9 @@ module constitutive constitutive_getAdiabaticTemperature, & constitutive_putAdiabaticTemperature, & constitutive_getTemperature, & + constitutive_getLocalVacancyConcentration, & + constitutive_putLocalVacancyConcentration, & + constitutive_getVacancyConcentration, & constitutive_postResults private :: & @@ -89,6 +94,8 @@ subroutine constitutive_init phase_damageInstance, & phase_thermal, & phase_thermalInstance, & + phase_vacancy, & + phase_vacancyInstance, & phase_Noutput, & homogenization_Ngrains, & homogenization_maxNgrains, & @@ -114,15 +121,20 @@ subroutine constitutive_init LOCAL_DAMAGE_gurson_ID, & LOCAL_THERMAL_isothermal_ID, & LOCAL_THERMAL_adiabatic_ID, & + LOCAL_VACANCY_constant_ID, & + LOCAL_VACANCY_generation_ID, & LOCAL_DAMAGE_NONE_label, & LOCAL_DAMAGE_brittle_label, & LOCAL_DAMAGE_ductile_label, & LOCAL_DAMAGE_gurson_label, & LOCAL_THERMAL_isothermal_label, & LOCAL_THERMAL_adiabatic_label, & + LOCAL_VACANCY_constant_label, & + LOCAL_VACANCY_generation_label, & plasticState, & damageState, & thermalState, & + vacancyState, & mappingConstitutive @@ -139,6 +151,9 @@ subroutine constitutive_init use damage_gurson use thermal_isothermal use thermal_adiabatic + use vacancy_constant + use vacancy_generation + implicit none integer(pInt), parameter :: FILEUNIT = 200_pInt integer(pInt) :: & @@ -150,7 +165,7 @@ subroutine constitutive_init integer(pInt), dimension(:) , pointer :: thisNoutput character(len=64), dimension(:,:), pointer :: thisOutput character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready - logical :: knownPlasticity, knownDamage, knownThermal, nonlocalConstitutionPresent + logical :: knownPlasticity, knownDamage, knownThermal, knownVacancy, nonlocalConstitutionPresent nonlocalConstitutionPresent = .false. !-------------------------------------------------------------------------------------------------- @@ -173,18 +188,26 @@ subroutine constitutive_init ! parse damage from config file if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file - if (any(phase_damage == LOCAL_DAMAGE_none_ID)) call damage_none_init(FILEUNIT) - if (any(phase_damage == LOCAL_DAMAGE_brittle_ID)) call damage_brittle_init(FILEUNIT) - if (any(phase_damage == LOCAL_DAMAGE_ductile_ID)) call damage_ductile_init(FILEUNIT) - if (any(phase_damage == LOCAL_DAMAGE_gurson_ID)) call damage_gurson_init(FILEUNIT) + if (any(phase_damage == LOCAL_DAMAGE_none_ID)) call damage_none_init(FILEUNIT) + if (any(phase_damage == LOCAL_DAMAGE_brittle_ID)) call damage_brittle_init(FILEUNIT) + if (any(phase_damage == LOCAL_DAMAGE_ductile_ID)) call damage_ductile_init(FILEUNIT) + if (any(phase_damage == LOCAL_DAMAGE_gurson_ID)) call damage_gurson_init(FILEUNIT) close(FILEUNIT) !-------------------------------------------------------------------------------------------------- ! parse thermal from config file if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file - if (any(phase_thermal == LOCAL_THERMAL_isothermal_ID)) call thermal_isothermal_init(FILEUNIT) - if (any(phase_thermal == LOCAL_THERMAL_adiabatic_ID)) call thermal_adiabatic_init(FILEUNIT) + if (any(phase_thermal == LOCAL_THERMAL_isothermal_ID)) call thermal_isothermal_init(FILEUNIT) + if (any(phase_thermal == LOCAL_THERMAL_adiabatic_ID)) call thermal_adiabatic_init(FILEUNIT) + close(FILEUNIT) + +!-------------------------------------------------------------------------------------------------- +! parse vacancy model from config file + if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... + call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file + if (any(phase_vacancy == LOCAL_VACANCY_constant_ID)) call vacancy_constant_init(FILEUNIT) + if (any(phase_vacancy == LOCAL_VACANCY_generation_ID)) call vacancy_generation_init(FILEUNIT) close(FILEUNIT) mainProcess: if (worldrank == 0) then @@ -286,7 +309,7 @@ subroutine constitutive_init instance = phase_thermalInstance(phase) ! which instance is present phase knownThermal = .true. select case(phase_thermal(phase)) ! split per constititution - case (LOCAL_THERMAL_ISOTHERMAL_ID) + case (LOCAL_THERMAL_isothermal_ID) outputName = LOCAL_THERMAL_ISOTHERMAL_label thisNoutput => null() thisOutput => null() @@ -307,6 +330,30 @@ subroutine constitutive_init enddo endif endif + instance = phase_vacancyInstance(phase) ! which instance is present phase + knownVacancy = .true. + select case(phase_vacancy(phase)) ! split per constititution + case (LOCAL_VACANCY_constant_ID) + outputName = LOCAL_VACANCY_constant_label + thisNoutput => null() + thisOutput => null() + thisSize => null() + case (LOCAL_VACANCY_generation_ID) + outputName = LOCAL_VACANCY_generation_label + thisNoutput => vacancy_generation_Noutput + thisOutput => vacancy_generation_output + thisSize => vacancy_generation_sizePostResult + case default + knownVacancy = .false. + end select + if (knownVacancy) then + write(FILEUNIT,'(a)') '(vacancy)'//char(9)//trim(outputName) + if (phase_vacancy(phase) /= LOCAL_VACANCY_generation_ID) then + do e = 1_pInt,thisNoutput(instance) + write(FILEUNIT,'(a,i4)') trim(thisOutput(e,instance))//char(9),thisSize(e,instance) + enddo + endif + endif #endif enddo close(FILEUNIT) @@ -317,6 +364,8 @@ subroutine constitutive_init constitutive_damage_maxSizeDotState = 0_pInt constitutive_thermal_maxSizePostResults = 0_pInt constitutive_thermal_maxSizeDotState = 0_pInt + constitutive_vacancy_maxSizePostResults = 0_pInt + constitutive_vacancy_maxSizeDotState = 0_pInt PhaseLoop2:do phase = 1_pInt,material_Nphase plasticState(phase)%partionedState0 = plasticState(phase)%State0 @@ -331,6 +380,10 @@ subroutine constitutive_init thermalState(phase)%State = thermalState(phase)%State0 constitutive_thermal_maxSizeDotState = max(constitutive_thermal_maxSizeDotState, thermalState(phase)%sizeDotState) constitutive_thermal_maxSizePostResults = max(constitutive_thermal_maxSizePostResults, thermalState(phase)%sizePostResults) + vacancyState(phase)%partionedState0 = vacancyState(phase)%State0 + vacancyState(phase)%State = vacancyState(phase)%State0 + constitutive_vacancy_maxSizeDotState = max(constitutive_vacancy_maxSizeDotState, vacancyState(phase)%sizeDotState) + constitutive_vacancy_maxSizePostResults = max(constitutive_vacancy_maxSizePostResults, vacancyState(phase)%sizePostResults) enddo PhaseLoop2 #ifdef HDF @@ -673,6 +726,7 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su phase_plasticity, & phase_damage, & phase_thermal, & + phase_vacancy, & material_phase, & homogenization_maxNgrains, & PLASTICITY_none_ID, & @@ -685,7 +739,8 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su LOCAL_DAMAGE_brittle_ID, & LOCAL_DAMAGE_ductile_ID, & LOCAL_DAMAGE_gurson_ID, & - LOCAL_THERMAL_adiabatic_ID + LOCAL_THERMAL_adiabatic_ID, & + LOCAL_VACANCY_generation_ID use constitutive_j2, only: & constitutive_j2_dotState use constitutive_phenopowerlaw, only: & @@ -706,6 +761,8 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su damage_gurson_dotState use thermal_adiabatic, only: & thermal_adiabatic_dotState + use vacancy_generation, only: & + vacancy_generation_dotState implicit none integer(pInt), intent(in) :: & @@ -761,6 +818,11 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su call thermal_adiabatic_dotState(Tstar_v, Lp, ipc, ip, el) end select + select case (phase_vacancy(material_phase(ipc,ip,el))) + case (LOCAL_VACANCY_generation_ID) + call vacancy_generation_dotState(Tstar_v, Lp, ipc, ip, el) + end select + if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) then call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) !$OMP CRITICAL (debugTimingDotState) @@ -1048,6 +1110,104 @@ function constitutive_getTemperature(ipc, ip, el) end function constitutive_getTemperature +!-------------------------------------------------------------------------------------------------- +!> @brief returns local vacancy concentration +!-------------------------------------------------------------------------------------------------- +function constitutive_getLocalVacancyConcentration(ipc, ip, el) + use prec, only: & + pReal + use material, only: & + material_phase, & + LOCAL_VACANCY_constant_ID, & + LOCAL_VACANCY_generation_ID, & + phase_vacancy + use vacancy_generation, only: & + vacancy_generation_getConcentration + use lattice, only: & + lattice_equilibriumVacancyConcentration + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal) :: constitutive_getLocalVacancyConcentration + + select case (phase_vacancy(material_phase(ipc,ip,el))) + case (LOCAL_VACANCY_constant_ID) + constitutive_getLocalVacancyConcentration = & + lattice_equilibriumVacancyConcentration(material_phase(ipc,ip,el)) + + case (LOCAL_VACANCY_generation_ID) + constitutive_getLocalVacancyConcentration = vacancy_generation_getConcentration(ipc, ip, el) + end select + +end function constitutive_getLocalVacancyConcentration + +!-------------------------------------------------------------------------------------------------- +!> @brief Puts local vacancy concentration +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_putLocalVacancyConcentration(ipc, ip, el, localVacancyConcentration) + use prec, only: & + pReal + use material, only: & + material_phase, & + LOCAL_VACANCY_generation_ID, & + phase_vacancy + use vacancy_generation, only: & + vacancy_generation_putConcentration + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + localVacancyConcentration + + select case (phase_vacancy(material_phase(ipc,ip,el))) + case (LOCAL_VACANCY_generation_ID) + call vacancy_generation_putConcentration(ipc, ip, el, localVacancyConcentration) + + end select + +end subroutine constitutive_putLocalVacancyConcentration + +!-------------------------------------------------------------------------------------------------- +!> @brief returns nonlocal vacancy concentration +!-------------------------------------------------------------------------------------------------- +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 + use lattice, only: & + lattice_equilibriumVacancyConcentration + implicit none + + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + 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 + +end function constitutive_getVacancyConcentration + !-------------------------------------------------------------------------------------------------- !> @brief returns accumulated slip on each system defined !-------------------------------------------------------------------------------------------------- @@ -1128,9 +1288,11 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) plasticState, & damageState, & thermalState, & + vacancyState, & phase_plasticity, & phase_damage, & phase_thermal, & + phase_vacancy, & material_phase, & homogenization_maxNgrains, & PLASTICITY_NONE_ID, & @@ -1143,7 +1305,8 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) LOCAL_DAMAGE_BRITTLE_ID, & LOCAL_DAMAGE_DUCTILE_ID, & LOCAL_DAMAGE_gurson_ID, & - LOCAL_THERMAL_ADIABATIC_ID + LOCAL_THERMAL_ADIABATIC_ID, & + LOCAL_VACANCY_generation_ID use constitutive_j2, only: & #ifdef HDF constitutive_j2_postResults2,& @@ -1168,6 +1331,8 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) damage_gurson_postResults use thermal_adiabatic, only: & thermal_adiabatic_postResults + use vacancy_generation, only: & + vacancy_generation_postResults #endif implicit none @@ -1178,7 +1343,8 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) #ifdef multiphysicsOut real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults + & damageState( material_phase(ipc,ip,el))%sizePostResults + & - thermalState(material_phase(ipc,ip,el))%sizePostResults) :: & + thermalState(material_phase(ipc,ip,el))%sizePostResults + & + vacancyState(material_phase(ipc,ip,el))%sizePostResults) :: & constitutive_postResults #else real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults) :: & @@ -1235,6 +1401,13 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) case (LOCAL_THERMAL_ADIABATIC_ID) constitutive_postResults(startPos:endPos) = thermal_adiabatic_postResults(ipc, ip, el) end select + + startPos = endPos + 1_pInt + endPos = endPos + vacancyState(material_phase(ipc,ip,el))%sizePostResults + select case (phase_vacancy(material_phase(ipc,ip,el))) + case (LOCAL_VACANCY_generation_ID) + constitutive_postResults(startPos:endPos) = vacancy_generation_postResults(ipc, ip, el) + end select #endif end function constitutive_postResults diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index f4d1010d2..bb7383d3c 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -3571,7 +3571,6 @@ subroutine constitutive_nonlocal_getAccumulatedSlip(nSlip,accumulatedSlip,ipc, i offset, & phase, & instance, & - offset_accshear_slip, & s offset = mappingConstitutive(1,ipc,ip,el) diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 169bfd589..368704414 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -45,7 +45,8 @@ module homogenization enum, bind(c) enumerator :: undefined_ID, & temperature_ID, & - damage_ID + damage_ID, & + vacancy_concentration_ID end enum integer(pInt), dimension(:), allocatable, private, protected :: & field_sizePostResults @@ -66,11 +67,15 @@ module homogenization field_putFieldDamage, & field_getLocalTemperature, & field_putFieldTemperature, & + field_getLocalVacancyConcentration, & + field_putFieldVacancyConcentration, & field_getDamageMobility, & field_getDamageDiffusion33, & field_getThermalConductivity33, & field_getMassDensity, & field_getSpecificHeat, & + field_getVacancyMobility, & + field_getVacancyDiffusion33, & materialpoint_postResults, & field_postResults private :: & @@ -112,7 +117,8 @@ subroutine homogenization_init() use constitutive, only: & constitutive_maxSizePostResults, & constitutive_damage_maxSizePostResults, & - constitutive_thermal_maxSizePostResults + constitutive_thermal_maxSizePostResults, & + constitutive_vacancy_maxSizePostResults use crystallite, only: & crystallite_maxSizePostResults use material @@ -203,6 +209,12 @@ subroutine homogenization_init() field_sizePostResult(field_Noutput(section),section) = 1_pInt field_sizePostResults(section) = field_sizePostResults(section) + 1_pInt field_output(field_Noutput(section),section) = IO_lc(IO_stringValue(line,positions,2_pInt)) + case('vacancy_concentration') + field_Noutput(section) = field_Noutput(section) + 1_pInt + field_outputID(field_Noutput(section),section) = vacancy_concentration_ID + field_sizePostResult(field_Noutput(section),section) = 1_pInt + field_sizePostResults(section) = field_sizePostResults(section) + 1_pInt + field_output(field_Noutput(section),section) = IO_lc(IO_stringValue(line,positions,2_pInt)) end select @@ -302,6 +314,7 @@ subroutine homogenization_init() #ifdef multiphysicsOut + constitutive_damage_maxSizePostResults & + constitutive_thermal_maxSizePostResults & + + constitutive_vacancy_maxSizePostResults & #endif + 1 + constitutive_maxSizePostResults) ! constitutive size & constitutive results allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems)) @@ -982,7 +995,6 @@ function field_getThermalConductivity33(ip,el) material_phase, & material_homog, & field_thermal_type, & - FIELD_THERMAL_local_ID, & FIELD_THERMAL_nonlocal_ID, & homogenization_Ngrains use crystallite, only: & @@ -1000,10 +1012,6 @@ function field_getThermalConductivity33(ip,el) field_getThermalConductivity33 =0.0_pReal select case(field_thermal_type(material_homog(ip,el))) - - case (FIELD_THERMAL_local_ID) - field_getThermalConductivity33 = 0.0_pReal - case (FIELD_THERMAL_nonlocal_ID) do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) field_getThermalConductivity33 = field_getThermalConductivity33 + & @@ -1027,7 +1035,6 @@ function field_getDamageDiffusion33(ip,el) material_phase, & material_homog, & field_damage_type, & - FIELD_DAMAGE_LOCAL_ID, & FIELD_DAMAGE_NONLOCAL_ID, & homogenization_Ngrains use crystallite, only: & @@ -1044,10 +1051,6 @@ function field_getDamageDiffusion33(ip,el) field_getDamageDiffusion33 =0.0_pReal select case(field_damage_type(material_homog(ip,el))) - - case (FIELD_DAMAGE_LOCAL_ID) - field_getDamageDiffusion33 = 0.0_pReal - case (FIELD_DAMAGE_NONLOCAL_ID) do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) field_getDamageDiffusion33 = field_getDamageDiffusion33 + & @@ -1059,6 +1062,7 @@ function field_getDamageDiffusion33(ip,el) field_getDamageDiffusion33 = field_getDamageDiffusion33 /homogenization_Ngrains(mesh_element(3,el)) end function field_getDamageDiffusion33 + !-------------------------------------------------------------------------------------------------- !> @brief Returns average mobility for damage field at each integration point !-------------------------------------------------------------------------------------------------- @@ -1071,7 +1075,6 @@ real(pReal) function field_getDamageMobility(ip,el) material_phase, & material_homog, & field_damage_type, & - FIELD_DAMAGE_LOCAL_ID, & FIELD_DAMAGE_NONLOCAL_ID, & homogenization_Ngrains @@ -1081,15 +1084,10 @@ real(pReal) function field_getDamageMobility(ip,el) el !< element number integer(pInt) :: & ipc - field_getDamageMobility =0.0_pReal select case(field_damage_type(material_homog(ip,el))) - - case (FIELD_DAMAGE_LOCAL_ID) - field_getDamageMobility = 0.0_pReal - case (FIELD_DAMAGE_NONLOCAL_ID) do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) field_getDamageMobility = field_getDamageMobility + lattice_DamageMobility(material_phase(ipc,ip,el)) @@ -1100,6 +1098,85 @@ real(pReal) function field_getDamageMobility(ip,el) field_getDamageMobility = field_getDamageMobility /homogenization_Ngrains(mesh_element(3,el)) end function field_getDamageMobility + +!-------------------------------------------------------------------------------------------------- +!> @brief Returns average diffusion tensor for vacancy field at each integration point +!-------------------------------------------------------------------------------------------------- +function field_getVacancyDiffusion33(ip,el) + use mesh, only: & + mesh_element + use lattice, only: & + lattice_vacancyDiffusion33 + use material, only: & + material_phase, & + material_homog, & + field_vacancy_type, & + FIELD_VACANCY_NONLOCAL_ID, & + homogenization_Ngrains + use crystallite, only: & + crystallite_push33ToRef + + implicit none + real(pReal), dimension(3,3) :: field_getVacancyDiffusion33 + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + integer(pInt) :: & + ipc + + field_getVacancyDiffusion33 = 0.0_pReal + + select case(field_vacancy_type(material_homog(ip,el))) + case (FIELD_VACANCY_NONLOCAL_ID) + do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) + field_getVacancyDiffusion33 = field_getVacancyDiffusion33 + & + crystallite_push33ToRef(ipc,ip,el,lattice_vacancyDiffusion33(:,:,material_phase(ipc,ip,el))) + enddo + + end select + + field_getVacancyDiffusion33 = field_getVacancyDiffusion33/ & + homogenization_Ngrains(mesh_element(3,el)) + +end function field_getVacancyDiffusion33 +!-------------------------------------------------------------------------------------------------- +!> @brief Returns average mobility for vacancy field at each integration point +!-------------------------------------------------------------------------------------------------- +real(pReal) function field_getVacancyMobility(ip,el) + use mesh, only: & + mesh_element + use lattice, only: & + lattice_vacancyMobility + use material, only: & + material_phase, & + material_homog, & + field_vacancy_type, & + FIELD_VACANCY_NONLOCAL_ID, & + homogenization_Ngrains + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + integer(pInt) :: & + ipc + + + field_getVacancyMobility =0.0_pReal + + select case(field_vacancy_type(material_homog(ip,el))) + case (FIELD_VACANCY_NONLOCAL_ID) + do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) + field_getVacancyMobility = field_getVacancyMobility + lattice_VacancyMobility(material_phase(ipc,ip,el)) + enddo + + end select + + field_getVacancyMobility = field_getVacancyMobility/ & + homogenization_Ngrains(mesh_element(3,el)) + +end function field_getVacancyMobility + !-------------------------------------------------------------------------------------------------- !> @brief ToDo !-------------------------------------------------------------------------------------------------- @@ -1212,6 +1289,62 @@ subroutine field_putFieldTemperature(ip,el,fieldThermalValue) end subroutine field_putFieldTemperature +!-------------------------------------------------------------------------------------------------- +!> @brief ToDo +!-------------------------------------------------------------------------------------------------- +real(pReal) function field_getLocalVacancyConcentration(ip,el) + use mesh, only: & + mesh_element + use material, only: & + homogenization_Ngrains + use constitutive, only: & + constitutive_getLocalVacancyConcentration + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + integer(pInt) :: & + ipc + + + field_getLocalVacancyConcentration = 0.0_pReal + do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) + field_getLocalVacancyConcentration = field_getLocalVacancyConcentration + & + constitutive_getLocalVacancyConcentration(ipc,ip,el) ! array/function/subroutine which is faster + enddo + field_getLocalVacancyConcentration = field_getLocalVacancyConcentration/ & + homogenization_Ngrains(mesh_element(3,el)) + +end function field_getLocalVacancyConcentration + +!-------------------------------------------------------------------------------------------------- +!> @brief Sets the diffused vacancy concentration in field state +!-------------------------------------------------------------------------------------------------- +subroutine field_putFieldVacancyConcentration(ip,el,fieldVacancyConcentration) + use material, only: & + material_homog, & + fieldVacancy, & + mappingHomogenization, & + field_vacancy_type, & + FIELD_VACANCY_nonlocal_ID + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el + real(pReal), intent(in) :: & + fieldVacancyConcentration + + select case(field_vacancy_type(material_homog(ip,el))) + case (FIELD_VACANCY_nonlocal_ID) + fieldVacancy(material_homog(ip,el))% & + field(1,mappingHomogenization(1,ip,el)) = fieldVacancyConcentration + + end select + +end subroutine field_putFieldVacancyConcentration + !-------------------------------------------------------------------------------------------------- !> @brief return array of homogenization results for post file inclusion. call only, !> if homogenization_sizePostResults(i,e) > 0 !! @@ -1266,7 +1399,8 @@ function field_postResults(ip,el) use material, only: & mappingHomogenization, & fieldThermal, & - fieldDamage + fieldDamage, & + fieldVacancy implicit none integer(pInt), intent(in) :: & @@ -1289,6 +1423,9 @@ function field_postResults(ip,el) case (damage_ID) field_postResults(c+1_pInt) = fieldDamage(homog)%field(1,pos) c = c + 1_pInt + case (vacancy_concentration_ID) + field_postResults(c+1_pInt) = fieldVacancy(homog)%field(1,pos) + c = c + 1_pInt end select enddo diff --git a/code/lattice.f90 b/code/lattice.f90 index 46e035652..00e639ea0 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -747,12 +747,15 @@ module lattice real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_thermalConductivity33, & lattice_thermalExpansion33, & - lattice_damageDiffusion33 + lattice_damageDiffusion33, & + lattice_vacancyDiffusion33 real(pReal), dimension(:), allocatable, public, protected :: & lattice_damageMobility, & + lattice_vacancyMobility, & lattice_massDensity, & lattice_specificHeat, & - lattice_referenceTemperature + lattice_referenceTemperature, & + lattice_equilibriumVacancyConcentration enum, bind(c) enumerator :: LATTICE_undefined_ID, & LATTICE_iso_ID, & @@ -1008,11 +1011,14 @@ subroutine lattice_init allocate(lattice_C3333(3,3,3,3,Nphases), source=0.0_pReal) allocate(lattice_thermalConductivity33(3,3,Nphases), source=0.0_pReal) allocate(lattice_thermalExpansion33 (3,3,Nphases), source=0.0_pReal) - allocate(lattice_damageDiffusion33 (3,3,Nphases), source=0.0_pReal) - allocate(lattice_damageMobility ( Nphases), source=0.0_pReal) - allocate(lattice_massDensity ( Nphases), source=0.0_pReal) - allocate(lattice_specificHeat ( Nphases), source=0.0_pReal) - allocate(lattice_referenceTemperature (Nphases), source=0.0_pReal) + allocate(lattice_damageDiffusion33 (3,3,Nphases), source=0.0_pReal) + allocate(lattice_vacancyDiffusion33 (3,3,Nphases), source=0.0_pReal) + allocate(lattice_damageMobility ( Nphases), source=0.0_pReal) + allocate(lattice_vacancyMobility ( Nphases), source=0.0_pReal) + allocate(lattice_massDensity ( Nphases), source=0.0_pReal) + allocate(lattice_specificHeat ( Nphases), source=0.0_pReal) + allocate(lattice_referenceTemperature ( Nphases), source=0.0_pReal) + allocate(lattice_equilibriumVacancyConcentration(Nphases), source=0.0_pReal) allocate(lattice_mu(Nphases), source=0.0_pReal) allocate(lattice_nu(Nphases), source=0.0_pReal) @@ -1143,6 +1149,16 @@ subroutine lattice_init lattice_DamageDiffusion33(3,3,section) = IO_floatValue(line,positions,2_pInt) case ('damage_mobility') lattice_DamageMobility(section) = IO_floatValue(line,positions,2_pInt) + case ('vacancy_diffusion11') + lattice_VacancyDiffusion33(1,1,section) = IO_floatValue(line,positions,2_pInt) + case ('vacancy_diffusion22') + lattice_VacancyDiffusion33(2,2,section) = IO_floatValue(line,positions,2_pInt) + case ('vacancy_diffusion33') + lattice_VacancyDiffusion33(3,3,section) = IO_floatValue(line,positions,2_pInt) + case ('vacancy_mobility') + lattice_VacancyMobility(section) = IO_floatValue(line,positions,2_pInt) + case ('equilibrium_vacancy_concentration') + lattice_equilibriumVacancyConcentration(section) = IO_floatValue(line,positions,2_pInt) end select endif enddo @@ -1232,6 +1248,8 @@ subroutine lattice_initializeStructure(myPhase,CoverA,aA,aM,cM) lattice_thermalExpansion33(1:3,1:3,myPhase)) lattice_DamageDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_DamageDiffusion33(1:3,1:3,myPhase)) + lattice_VacancyDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_VacancyDiffusion33(1:3,1:3,myPhase)) select case(lattice_structure(myPhase)) !-------------------------------------------------------------------------------------------------- diff --git a/code/material.f90 b/code/material.f90 index ba8d288db..4788a4723 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -33,10 +33,14 @@ module material LOCAL_DAMAGE_gurson_label = 'gurson', & LOCAL_THERMAL_isothermal_label = 'isothermal', & LOCAL_THERMAL_adiabatic_label = 'adiabatic', & + LOCAL_VACANCY_constant_label = 'constant', & + LOCAL_VACANCY_generation_label = 'generation', & FIELD_DAMAGE_local_label = 'local', & FIELD_DAMAGE_nonlocal_label = 'nonlocal', & FIELD_THERMAL_local_label = 'local', & FIELD_THERMAL_nonlocal_label = 'nonlocal', & + FIELD_VACANCY_local_label = 'local', & + FIELD_VACANCY_nonlocal_label = 'nonlocal', & HOMOGENIZATION_none_label = 'none', & HOMOGENIZATION_isostrain_label = 'isostrain', & HOMOGENIZATION_rgc_label = 'rgc' @@ -57,18 +61,21 @@ module material PLASTICITY_titanmod_ID, & PLASTICITY_nonlocal_ID end enum - enum, bind(c) enumerator :: LOCAL_DAMAGE_none_ID, & LOCAL_DAMAGE_brittle_ID, & LOCAL_DAMAGE_ductile_ID, & LOCAL_DAMAGE_gurson_ID end enum - enum, bind(c) enumerator :: LOCAL_THERMAL_isothermal_ID, & LOCAL_THERMAL_adiabatic_ID end enum + enum, bind(c) + enumerator :: LOCAL_VACANCY_constant_ID, & + LOCAL_VACANCY_generation_ID + end enum + enum, bind(c) enumerator :: FIELD_DAMAGE_local_ID ,& FIELD_DAMAGE_nonlocal_ID @@ -78,6 +85,10 @@ module material enumerator :: FIELD_THERMAL_local_ID, & FIELD_THERMAL_nonlocal_ID end enum + enum, bind(c) + enumerator :: FIELD_VACANCY_local_ID, & + FIELD_VACANCY_nonlocal_ID + end enum enum, bind(c) enumerator :: HOMOGENIZATION_undefined_ID, & HOMOGENIZATION_none_ID, & @@ -94,18 +105,22 @@ module material MATERIAL_partCrystallite = 'crystallite', & !< keyword for crystallite part MATERIAL_partPhase = 'phase' !< keyword for phase part - integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable, public, protected :: & + integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable, public, protected :: & phase_elasticity !< elasticity of each phase - integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable, public, protected :: & + integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable, public, protected :: & phase_plasticity !< plasticity of each phase - integer(kind(LOCAL_DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & + integer(kind(LOCAL_DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & phase_damage !< local damage of each phase integer(kind(LOCAL_THERMAL_isothermal_ID)), dimension(:), allocatable, public, protected :: & phase_thermal !< local thermal of each phase - integer(kind(FIELD_DAMAGE_local_ID)), dimension(:), allocatable, public, protected :: & - field_damage_type !< field damage of each phase - integer(kind(FIELD_THERMAL_local_ID)), dimension(:), allocatable, public, protected :: & - field_thermal_type !< field thermal of each phase + integer(kind(LOCAL_VACANCY_constant_ID)), dimension(:), allocatable, public, protected :: & + phase_vacancy !< local vacancy model of each phase + integer(kind(FIELD_DAMAGE_local_ID)), dimension(:), allocatable, public, protected :: & + field_damage_type !< field damage of each phase + integer(kind(FIELD_THERMAL_local_ID)), dimension(:), allocatable, public, protected :: & + field_thermal_type !< field thermal of each phase + integer(kind(FIELD_VACANCY_local_ID)), dimension(:), allocatable, public, protected :: & + field_vacancy_type !< field vacancy of each phase integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: & homogenization_type !< type of each homogenization @@ -130,6 +145,7 @@ module material phase_plasticityInstance, & !< instance of particular plasticity of each phase phase_damageInstance, & !< instance of particular damage of each phase phase_thermalInstance, & !< instance of particular thermal of each phase + phase_vacancyInstance, & !< instance of particular vacancy model of each phase crystallite_Noutput, & !< number of '(output)' items per crystallite setting homogenization_typeInstance, & !< instance of particular type of each homogenization microstructure_crystallite !< crystallite setting ID of each microstructure @@ -142,12 +158,15 @@ module material plasticState, & damageState, & thermalState,& + vacancyState,& homogState type(tFieldData), allocatable, dimension(:), public :: & fieldDamage type(tFieldData), allocatable, dimension(:), public :: & fieldThermal + type(tFieldData), allocatable, dimension(:), public :: & + fieldVacancy integer(pInt), dimension(:,:,:), allocatable, public, protected :: & @@ -221,10 +240,14 @@ module material LOCAL_DAMAGE_gurson_ID, & LOCAL_THERMAL_isothermal_ID, & LOCAL_THERMAL_adiabatic_ID, & + LOCAL_VACANCY_constant_ID, & + LOCAL_VACANCY_generation_ID, & FIELD_DAMAGE_local_ID, & FIELD_DAMAGE_nonlocal_ID, & FIELD_THERMAL_local_ID, & FIELD_THERMAL_nonlocal_ID, & + FIELD_VACANCY_local_ID, & + FIELD_VACANCY_nonlocal_ID, & HOMOGENIZATION_none_ID, & HOMOGENIZATION_isostrain_ID, & #ifdef HDF @@ -309,9 +332,11 @@ subroutine material_init allocate(plasticState(material_Nphase)) allocate(damageState (material_Nphase)) allocate(thermalState(material_Nphase)) + allocate(vacancyState(material_Nphase)) allocate(homogState (material_Nhomogenization)) allocate(fieldDamage (material_Nhomogenization)) allocate(fieldThermal(material_Nhomogenization)) + allocate(fieldVacancy(material_Nhomogenization)) do m = 1_pInt,material_Nmicrostructure if(microstructure_crystallite(m) < 1_pInt .or. & microstructure_crystallite(m) > material_Ncrystallite) & @@ -406,6 +431,24 @@ subroutine material_init end select enddo + do homog = 1,material_Nhomogenization + NofMyField=count(material_homog==homog) + select case(field_vacancy_type(homog)) + case (FIELD_VACANCY_local_ID) + fieldVacancy(homog)%sizeField = 0_pInt + fieldVacancy(homog)%sizePostResults = 0_pInt + allocate(fieldVacancy(homog)%field(fieldVacancy(homog)%sizeField,NofMyField), & + source = 0.0_pReal) + + case (FIELD_VACANCY_nonlocal_ID) + fieldVacancy(homog)%sizeField = 1_pInt + fieldVacancy(homog)%sizePostResults = 1_pInt + allocate(fieldVacancy(homog)%field(fieldVacancy(homog)%sizeField,NofMyField), & + source = 0.0_pReal) + + end select + enddo + end subroutine material_init @@ -450,6 +493,7 @@ subroutine material_parseHomogenization(fileUnit,myPart) allocate(homogenization_type(Nsections), source=HOMOGENIZATION_undefined_ID) allocate(FIELD_DAMAGE_type(Nsections), source=FIELD_DAMAGE_local_ID) allocate(FIELD_THERMAL_type(Nsections), source=FIELD_THERMAL_local_ID) + allocate(FIELD_VACANCY_type(Nsections), source=FIELD_VACANCY_local_ID) allocate(homogenization_typeInstance(Nsections), source=0_pInt) allocate(homogenization_Ngrains(Nsections), source=0_pInt) allocate(homogenization_Noutput(Nsections), source=0_pInt) @@ -515,6 +559,17 @@ subroutine material_parseHomogenization(fileUnit,myPart) case default call IO_error(500_pInt,ext_msg=trim(IO_stringValue(line,positions,2_pInt))) end select + + case ('field_vacancy') + select case (IO_lc(IO_stringValue(line,positions,2_pInt))) + case(FIELD_VACANCY_local_label) + FIELD_VACANCY_type(section) = FIELD_VACANCY_local_ID + case(FIELD_VACANCY_nonlocal_label) + FIELD_VACANCY_type(section) = FIELD_VACANCY_nonlocal_ID + case default + call IO_error(500_pInt,ext_msg=trim(IO_stringValue(line,positions,2_pInt))) + end select + case ('nconstituents','ngrains') homogenization_Ngrains(section) = IO_intValue(line,positions,2_pInt) end select @@ -726,6 +781,8 @@ subroutine material_parsePhase(fileUnit,myPart) allocate(phase_damageInstance(Nsections), source=0_pInt) allocate(phase_thermal(Nsections) , source=LOCAL_THERMAL_isothermal_ID) allocate(phase_thermalInstance(Nsections), source=0_pInt) + allocate(phase_vacancy(Nsections) , source=LOCAL_VACANCY_constant_ID) + allocate(phase_vacancyInstance(Nsections), source=0_pInt) allocate(phase_Noutput(Nsections), source=0_pInt) allocate(phase_localPlasticity(Nsections), source=.false.) @@ -807,7 +864,16 @@ subroutine material_parsePhase(fileUnit,myPart) case default call IO_error(200_pInt,ext_msg=trim(IO_stringValue(line,positions,2_pInt))) end select - phase_thermalInstance(section) = count(phase_thermal(1:section) == phase_thermal(section)) ! count instances + case ('vacancy') + select case (IO_lc(IO_stringValue(line,positions,2_pInt))) + case (LOCAL_VACANCY_CONSTANT_label) + phase_vacancy(section) = LOCAL_VACANCY_constant_ID + case (LOCAL_VACANCY_GENERATION_label) + phase_vacancy(section) = LOCAL_VACANCY_generation_ID + case default + call IO_error(200_pInt,ext_msg=trim(IO_stringValue(line,positions,2_pInt))) + end select + phase_vacancyInstance(section) = count(phase_vacancy(1:section) == phase_vacancy(section)) ! count instances end select endif enddo diff --git a/code/thermal_adiabatic.f90 b/code/thermal_adiabatic.f90 index 325ff2dd7..734627956 100644 --- a/code/thermal_adiabatic.f90 +++ b/code/thermal_adiabatic.f90 @@ -1,9 +1,8 @@ !-------------------------------------------------------------------------------------------------- ! $Id: thermal_adiabatic.f90 3210 2014-06-17 15:24:44Z MPIE\m.diehl $ !-------------------------------------------------------------------------------------------------- -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine incoprorating dislocation and twinning physics +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine incoprorating local heat generation due to plastic dissipation !> @details to be done !-------------------------------------------------------------------------------------------------- module thermal_adiabatic diff --git a/code/vacancy_constant.f90 b/code/vacancy_constant.f90 new file mode 100644 index 000000000..806b0934e --- /dev/null +++ b/code/vacancy_constant.f90 @@ -0,0 +1,103 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: vacancy_constant.f90 3148 2014-05-27 14:46:03Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for constant vacancy concentration +!-------------------------------------------------------------------------------------------------- +module vacancy_constant + use prec, only: & + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + vacancy_constant_sizePostResults + + integer(pInt), dimension(:,:), allocatable, target, public :: & + vacancy_constant_sizePostResult !< size of each post result output + + public :: & + vacancy_constant_init + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine vacancy_constant_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 IO, only: & + IO_timeStamp + use numerics, only: & + worldrank, & + numerics_integrator + use material, only: & + phase_vacancy, & + phase_Noutput, & + LOCAL_VACANCY_CONSTANT_label, & + LOCAL_VACANCY_CONSTANT_ID, & + material_phase, & + vacancyState, & + MATERIAL_partPhase + + implicit none + + integer(pInt), intent(in) :: fileUnit + integer(pInt) :: & + maxNinstance, & + phase, & + NofMyPhase, & + sizeState, & + sizeDotState + + mainProcess: if (worldrank == 0) then + write(6,'(/,a)') ' <<<+- vacancy_'//LOCAL_VACANCY_CONSTANT_label//' init -+>>>' + write(6,'(a)') ' $Id: vacancy_constant.f90 3148 2014-05-27 14:46:03Z MPIE\m.diehl $' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + + maxNinstance = int(count(phase_vacancy == LOCAL_VACANCY_CONSTANT_ID),pInt) + if (maxNinstance == 0_pInt) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + + initializeInstances: do phase = 1_pInt, size(phase_vacancy) + NofMyPhase=count(material_phase==phase) + + if (phase_vacancy(phase) == LOCAL_VACANCY_CONSTANT_ID) then + sizeState = 0_pInt + vacancyState(phase)%sizeState = sizeState + sizeDotState = sizeState + vacancyState(phase)%sizeDotState = sizeDotState + vacancyState(phase)%sizePostResults = 0_pInt + allocate(vacancyState(phase)%state0 (sizeState,NofMyPhase)) + allocate(vacancyState(phase)%partionedState0(sizeState,NofMyPhase)) + allocate(vacancyState(phase)%subState0 (sizeState,NofMyPhase)) + allocate(vacancyState(phase)%state (sizeState,NofMyPhase)) + allocate(vacancyState(phase)%state_backup (sizeState,NofMyPhase)) + allocate(vacancyState(phase)%aTolState (NofMyPhase)) + allocate(vacancyState(phase)%dotState (sizeDotState,NofMyPhase)) + allocate(vacancyState(phase)%dotState_backup(sizeDotState,NofMyPhase)) + if (any(numerics_integrator == 1_pInt)) then + allocate(vacancyState(phase)%previousDotState (sizeDotState,NofMyPhase)) + allocate(vacancyState(phase)%previousDotState2 (sizeDotState,NofMyPhase)) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(vacancyState(phase)%RK4dotState (sizeDotState,NofMyPhase)) + if (any(numerics_integrator == 5_pInt)) & + allocate(vacancyState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase)) + endif + enddo initializeInstances + allocate(vacancy_constant_sizePostResults(maxNinstance), source=0_pInt) + +end subroutine vacancy_constant_init + +end module vacancy_constant diff --git a/code/vacancy_generation.f90 b/code/vacancy_generation.f90 new file mode 100644 index 000000000..3cbf2b739 --- /dev/null +++ b/code/vacancy_generation.f90 @@ -0,0 +1,365 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: vacancy_generation.f90 3210 2014-06-17 15:24:44Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for plastically generated vacancy concentrations +!> @details to be done +!-------------------------------------------------------------------------------------------------- +module vacancy_generation + use prec, only: & + pReal, & + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + vacancy_generation_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + vacancy_generation_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + vacancy_generation_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, target, public :: & + vacancy_generation_Noutput !< number of outputs per instance of this damage + + real(pReal), dimension(:), allocatable, public :: & + vacancy_generation_aTol + + enum, bind(c) + enumerator :: undefined_ID, & + vacancy_concentration_ID + end enum + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + vacancy_generation_outputID !< ID of each post result output + + + public :: & + vacancy_generation_init, & + vacancy_generation_stateInit, & + vacancy_generation_aTolState, & + vacancy_generation_dotState, & + vacancy_generation_getConcentration, & + vacancy_generation_putConcentration, & + vacancy_generation_postResults + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine vacancy_generation_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_vacancy, & + phase_vacancyInstance, & + phase_Noutput, & + LOCAL_VACANCY_GENERATION_label, & + LOCAL_VACANCY_generation_ID, & + material_phase, & + vacancyState, & + 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)') ' <<<+- vacancy_'//LOCAL_VACANCY_GENERATION_label//' init -+>>>' + write(6,'(a)') ' $Id: vacancy_generation.f90 3210 2014-06-17 15:24:44Z MPIE\m.diehl $' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + + maxNinstance = int(count(phase_vacancy == LOCAL_VACANCY_generation_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(vacancy_generation_sizePostResults(maxNinstance), source=0_pInt) + allocate(vacancy_generation_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) + allocate(vacancy_generation_output(maxval(phase_Noutput),maxNinstance)) + vacancy_generation_output = '' + allocate(vacancy_generation_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) + allocate(vacancy_generation_Noutput(maxNinstance), source=0_pInt) + allocate(vacancy_generation_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_vacancy(phase) == LOCAL_VACANCY_generation_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran + + instance = phase_vacancyInstance(phase) ! which instance of my vacancy 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 ('vacancy_concentration') + vacancy_generation_Noutput(instance) = vacancy_generation_Noutput(instance) + 1_pInt + vacancy_generation_outputID(vacancy_generation_Noutput(instance),instance) = vacancy_concentration_ID + vacancy_generation_output(vacancy_generation_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + end select + + case ('atol_vacancyGeneration') + vacancy_generation_aTol(instance) = IO_floatValue(line,positions,2_pInt) + + end select + endif; endif + enddo parsingFile + + initializeInstances: do phase = 1_pInt, size(phase_vacancy) + if (phase_vacancy(phase) == LOCAL_VACANCY_generation_ID) then + NofMyPhase=count(material_phase==phase) + instance = phase_vacancyInstance(phase) + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputsLoop: do o = 1_pInt,vacancy_generation_Noutput(instance) + select case(vacancy_generation_outputID(o,instance)) + case(vacancy_concentration_ID) + mySize = 1_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + vacancy_generation_sizePostResult(o,instance) = mySize + vacancy_generation_sizePostResults(instance) = vacancy_generation_sizePostResults(instance) + mySize + endif + enddo outputsLoop +! Determine size of state array + sizeDotState = 1_pInt + sizeState = 1_pInt + vacancyState(phase)%sizeState = sizeState + vacancyState(phase)%sizeDotState = sizeDotState + vacancyState(phase)%sizePostResults = vacancy_generation_sizePostResults(instance) + allocate(vacancyState(phase)%aTolState (sizeState), source=0.0_pReal) + allocate(vacancyState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(vacancyState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(vacancyState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(vacancyState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal) + allocate(vacancyState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) + + allocate(vacancyState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(vacancyState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(vacancyState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 1_pInt)) then + allocate(vacancyState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(vacancyState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(vacancyState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 5_pInt)) & + allocate(vacancyState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) + + call vacancy_generation_stateInit(phase,instance) + call vacancy_generation_aTolState(phase,instance) + endif + + enddo initializeInstances +end subroutine vacancy_generation_init + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant NEW state values for a given instance of this vacancy model +!-------------------------------------------------------------------------------------------------- +subroutine vacancy_generation_stateInit(phase,instance) + use material, only: & + vacancyState + use lattice, only: & + lattice_equilibriumVacancyConcentration + + implicit none + integer(pInt), intent(in) :: instance !< number specifying the instance of the vacancy + integer(pInt), intent(in) :: phase !< number specifying the phase of the vacancy + + real(pReal), dimension(vacancyState(phase)%sizeState) :: tempState + + tempState(1) = lattice_equilibriumVacancyConcentration(phase) + vacancyState(phase)%state = spread(tempState,2,size(vacancyState(phase)%state(1,:))) + vacancyState(phase)%state0 = vacancyState(phase)%state + vacancyState(phase)%partionedState0 = vacancyState(phase)%state +end subroutine vacancy_generation_stateInit + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant state values for a given instance of this vacancy model +!-------------------------------------------------------------------------------------------------- +subroutine vacancy_generation_aTolState(phase,instance) + use material, only: & + vacancyState + + implicit none + integer(pInt), intent(in) :: & + phase, & + instance ! number specifying the current instance of the vacancy + real(pReal), dimension(vacancyState(phase)%sizeState) :: tempTol + + tempTol = vacancy_generation_aTol + vacancyState(phase)%aTolState = tempTol +end subroutine vacancy_generation_aTolState + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine vacancy_generation_dotState(Tstar_v, Lp, ipc, ip, el) + use lattice, only: & + lattice_massDensity, & + lattice_specificHeat + use material, only: & + mappingConstitutive, & + phase_vacancyInstance, & + vacancyState + use math, only: & + math_Mandel6to33 + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + real(pReal), intent(in), dimension(3,3) :: & + Lp + integer(pInt) :: & + instance, phase, constituent + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_vacancyInstance(phase) + + vacancyState(phase)%dotState(1,constituent) = & + 0.95_pReal & + * sum(abs(math_Mandel6to33(Tstar_v)*Lp)) & + / (lattice_massDensity(phase)*lattice_specificHeat(phase)) + +end subroutine vacancy_generation_dotState + +!-------------------------------------------------------------------------------------------------- +!> @brief returns vacancy concentration based on state layout +!-------------------------------------------------------------------------------------------------- +function vacancy_generation_getConcentration(ipc, ip, el) + use material, only: & + mappingConstitutive, & + vacancyState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal) :: vacancy_generation_getConcentration + + vacancy_generation_getConcentration = & + vacancyState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el)) + +end function vacancy_generation_getConcentration + +!-------------------------------------------------------------------------------------------------- +!> @brief returns temperature based on local damage model state layout +!-------------------------------------------------------------------------------------------------- +subroutine vacancy_generation_putConcentration(ipc, ip, el, localVacancyConcentration) + use material, only: & + mappingConstitutive, & + vacancyState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + localVacancyConcentration + + vacancyState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el))= & + localVacancyConcentration + +end subroutine vacancy_generation_putConcentration + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +function vacancy_generation_postResults(ipc,ip,el) + use material, only: & + mappingConstitutive, & + phase_vacancyInstance, & + vacancyState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), dimension(vacancy_generation_sizePostResults(phase_vacancyInstance(mappingConstitutive(2,ipc,ip,el)))) :: & + vacancy_generation_postResults + + integer(pInt) :: & + instance, phase, constituent, o, c + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_vacancyInstance(phase) + + c = 0_pInt + vacancy_generation_postResults = 0.0_pReal + + do o = 1_pInt,vacancy_generation_Noutput(instance) + select case(vacancy_generation_outputID(o,instance)) + + case (vacancy_concentration_ID) + vacancy_generation_postResults(c+1_pInt) = vacancyState(phase)%state(1,constituent) + c = c + 1 + end select + enddo +end function vacancy_generation_postResults + +end module vacancy_generation