From 87d42bf4475771dab4ba09bd036c6a2c8c24d9c3 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Fri, 24 Jul 2015 14:47:18 +0000 Subject: [PATCH] initialise Fi correctly for initial field values away from equilibrium --- code/constitutive.f90 | 59 +++++++++++++++++++++++++++ code/crystallite.f90 | 6 ++- code/kinematics_hydrogen_strain.f90 | 38 +++++++++++++++++ code/kinematics_thermal_expansion.f90 | 12 +++--- code/kinematics_vacancy_strain.f90 | 37 +++++++++++++++++ code/lattice.f90 | 8 +++- 6 files changed, 149 insertions(+), 11 deletions(-) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 201f7d199..5c9015f02 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -23,6 +23,7 @@ module constitutive constitutive_microstructure, & constitutive_LpAndItsTangent, & constitutive_LiAndItsTangent, & + constitutive_initialFi, & constitutive_TandItsTangent, & constitutive_collectDotState, & constitutive_collectDeltaState, & @@ -713,6 +714,64 @@ enddo end subroutine constitutive_LiAndItsTangent +!-------------------------------------------------------------------------------------------------- +!> @brief collects initial intermediate deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function constitutive_initialFi(ipc, ip, el) + use prec, only: & + pReal + use math, only: & + math_I3, & + math_inv33, & + math_mul33x33 + use material, only: & + phase_kinematics, & + phase_Nkinematics, & + material_phase, & + KINEMATICS_thermal_expansion_ID, & + KINEMATICS_vacancy_strain_ID, & + KINEMATICS_hydrogen_strain_ID + use kinematics_thermal_expansion, only: & + kinematics_thermal_expansion_initialStrain + use kinematics_vacancy_strain, only: & + kinematics_vacancy_strain_initialStrain + use kinematics_hydrogen_strain, only: & + kinematics_hydrogen_strain_initialStrain + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + constitutive_initialFi !< composite initial intermediate deformation gradient + real(pReal), dimension(3,3) :: & + my_Fi !< individual intermediate deformation gradients + integer(pInt) :: & + kinematics + + constitutive_initialFi = math_I3 + + do kinematics = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) !< Warning: small initial strain assumption + select case (phase_kinematics(kinematics,material_phase(ipc,ip,el))) + case (KINEMATICS_thermal_expansion_ID) + constitutive_initialFi = & + constitutive_initialFi + kinematics_thermal_expansion_initialStrain(ipc, ip, el) + + case (KINEMATICS_vacancy_strain_ID) + constitutive_initialFi = & + constitutive_initialFi + kinematics_vacancy_strain_initialStrain(ipc, ip, el) + + case (KINEMATICS_hydrogen_strain_ID) + constitutive_initialFi = & + constitutive_initialFi + kinematics_hydrogen_strain_initialStrain(ipc, ip, el) + + end select +enddo + +end function constitutive_initialFi + + !-------------------------------------------------------------------------------------------------- !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> the elastic deformation gradient depending on the selected elastic law (so far no case switch diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 019d1f96a..699c5c177 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -179,6 +179,7 @@ subroutine crystallite_init IO_EOF use material use constitutive, only: & + constitutive_initialFi, & constitutive_microstructure ! derived (shortcut) quantities of given state implicit none @@ -415,10 +416,11 @@ subroutine crystallite_init myNgrains = homogenization_Ngrains(mesh_element(3,e)) ! look up homogenization-->grainCount forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1_pInt:myNgrains) crystallite_Fp0(1:3,1:3,g,i,e) = math_EulerToR(material_EulerAngles(1:3,g,i,e)) ! plastic def gradient reflects init orientation - crystallite_Fi0(1:3,1:3,g,i,e) = math_I3 + crystallite_Fi0(1:3,1:3,g,i,e) = constitutive_initialFi(g,i,e) crystallite_F0(1:3,1:3,g,i,e) = math_I3 crystallite_localPlasticity(g,i,e) = phase_localPlasticity(material_phase(g,i,e)) - crystallite_Fe(1:3,1:3,g,i,e) = math_transpose33(crystallite_Fp0(1:3,1:3,g,i,e)) + crystallite_Fe(1:3,1:3,g,i,e) = math_inv33(math_mul33x33(crystallite_Fi0(1:3,1:3,g,i,e), & + crystallite_Fp0(1:3,1:3,g,i,e))) ! assuming that euler angles are given in internal strain free configuration crystallite_Fp(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) crystallite_Fi(1:3,1:3,g,i,e) = crystallite_Fi0(1:3,1:3,g,i,e) crystallite_requested(g,i,e) = .true. diff --git a/code/kinematics_hydrogen_strain.f90 b/code/kinematics_hydrogen_strain.f90 index 79f6677e5..47e6d7724 100644 --- a/code/kinematics_hydrogen_strain.f90 +++ b/code/kinematics_hydrogen_strain.f90 @@ -31,6 +31,7 @@ module kinematics_hydrogen_strain public :: & kinematics_hydrogen_strain_init, & + kinematics_hydrogen_strain_initialStrain, & kinematics_hydrogen_strain_LiAndItsTangent, & kinematics_hydrogen_strain_ChemPotAndItsTangent @@ -142,6 +143,43 @@ subroutine kinematics_hydrogen_strain_init(fileUnit) end subroutine kinematics_hydrogen_strain_init +!-------------------------------------------------------------------------------------------------- +!> @brief report initial hydrogen strain based on current hydrogen conc deviation from +!> equillibrium (0) +!-------------------------------------------------------------------------------------------------- +pure function kinematics_hydrogen_strain_initialStrain(ipc, ip, el) + use math, only: & + math_I3 + use material, only: & + material_phase, & + material_homog, & + hydrogenConc, & + hydrogenfluxMapping + use lattice, only: & + lattice_equilibriumHydrogenConcentration + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + kinematics_hydrogen_strain_initialStrain !< initial thermal strain (should be small strain, though) + integer(pInt) :: & + phase, & + homog, offset, instance + + phase = material_phase(ipc,ip,el) + instance = kinematics_hydrogen_strain_instance(phase) + homog = material_homog(ip,el) + offset = hydrogenfluxMapping(homog)%p(ip,el) + + kinematics_hydrogen_strain_initialStrain = & + (hydrogenConc(homog)%p(offset) - lattice_equilibriumHydrogenConcentration(phase)) * & + kinematics_hydrogen_strain_coeff(instance)* math_I3 + +end function kinematics_hydrogen_strain_initialStrain + !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- diff --git a/code/kinematics_thermal_expansion.f90 b/code/kinematics_thermal_expansion.f90 index 1e2560890..b4119a31e 100644 --- a/code/kinematics_thermal_expansion.f90 +++ b/code/kinematics_thermal_expansion.f90 @@ -32,7 +32,7 @@ module kinematics_thermal_expansion ! end enum public :: & kinematics_thermal_expansion_init, & - kinematics_thermal_expansion_initialFi, & + kinematics_thermal_expansion_initialStrain, & kinematics_thermal_expansion_LiAndItsTangent contains @@ -153,14 +153,12 @@ end subroutine kinematics_thermal_expansion_init !-------------------------------------------------------------------------------------------------- !> @brief report initial thermal strain based on current temperature deviation from reference !-------------------------------------------------------------------------------------------------- -pure function kinematics_thermal_expansion_initialFi(ipc, ip, el) +pure function kinematics_thermal_expansion_initialStrain(ipc, ip, el) use material, only: & material_phase, & material_homog, & temperature, & thermalMapping -use math, only: & - math_I3 use lattice, only: & lattice_thermalExpansion33, & lattice_referenceTemperature @@ -171,7 +169,7 @@ use math, only: & ip, & !< integration point number el !< element number real(pReal), dimension(3,3) :: & - kinematics_thermal_expansion_initialFi !< initial thermal strain (should be small strain, though) + kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) integer(pInt) :: & phase, & homog, offset @@ -180,11 +178,11 @@ use math, only: & homog = material_homog(ip,el) offset = thermalMapping(homog)%p(ip,el) - kinematics_thermal_expansion_initialFi = math_I3 + & + kinematics_thermal_expansion_initialStrain = & (temperature(homog)%p(offset) - lattice_referenceTemperature(phase)) * & lattice_thermalExpansion33(1:3,1:3,phase) -end function kinematics_thermal_expansion_initialFi +end function kinematics_thermal_expansion_initialStrain !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient diff --git a/code/kinematics_vacancy_strain.f90 b/code/kinematics_vacancy_strain.f90 index 403b1b730..f013acee4 100644 --- a/code/kinematics_vacancy_strain.f90 +++ b/code/kinematics_vacancy_strain.f90 @@ -31,6 +31,7 @@ module kinematics_vacancy_strain public :: & kinematics_vacancy_strain_init, & + kinematics_vacancy_strain_initialStrain, & kinematics_vacancy_strain_LiAndItsTangent, & kinematics_vacancy_strain_ChemPotAndItsTangent @@ -142,6 +143,42 @@ subroutine kinematics_vacancy_strain_init(fileUnit) end subroutine kinematics_vacancy_strain_init +!-------------------------------------------------------------------------------------------------- +!> @brief report initial vacancy strain based on current vacancy conc deviation from equillibrium +!-------------------------------------------------------------------------------------------------- +pure function kinematics_vacancy_strain_initialStrain(ipc, ip, el) + use math, only: & + math_I3 + use material, only: & + material_phase, & + material_homog, & + vacancyConc, & + vacancyfluxMapping + use lattice, only: & + lattice_equilibriumVacancyConcentration + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + kinematics_vacancy_strain_initialStrain !< initial thermal strain (should be small strain, though) + integer(pInt) :: & + phase, & + homog, offset, instance + + phase = material_phase(ipc,ip,el) + instance = kinematics_vacancy_strain_instance(phase) + homog = material_homog(ip,el) + offset = vacancyfluxMapping(homog)%p(ip,el) + + kinematics_vacancy_strain_initialStrain = & + (vacancyConc(homog)%p(offset) - lattice_equilibriumVacancyConcentration(phase)) * & + kinematics_vacancy_strain_coeff(instance)* math_I3 + +end function kinematics_vacancy_strain_initialStrain + !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- diff --git a/code/lattice.f90 b/code/lattice.f90 index 774ad6713..87b95754c 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -1041,7 +1041,8 @@ module lattice lattice_hydrogenSurfaceEnergy, & lattice_hydrogenVol, & lattice_referenceTemperature, & - lattice_equilibriumVacancyConcentration + lattice_equilibriumVacancyConcentration, & + lattice_equilibriumHydrogenConcentration enum, bind(c) enumerator :: LATTICE_undefined_ID, & LATTICE_iso_ID, & @@ -1344,8 +1345,9 @@ subroutine lattice_init allocate(lattice_hydrogenFormationEnergy( Nphases), source=0.0_pReal) allocate(lattice_hydrogenSurfaceEnergy ( Nphases), source=0.0_pReal) allocate(lattice_hydrogenVol ( Nphases), source=0.0_pReal) - allocate(lattice_referenceTemperature ( Nphases), source=0.0_pReal) + allocate(lattice_referenceTemperature ( Nphases), source=300.0_pReal) allocate(lattice_equilibriumVacancyConcentration(Nphases), source=0.0_pReal) + allocate(lattice_equilibriumHydrogenConcentration(Nphases),source=0.0_pReal) allocate(lattice_mu(Nphases), source=0.0_pReal) allocate(lattice_nu(Nphases), source=0.0_pReal) @@ -1551,6 +1553,8 @@ subroutine lattice_init lattice_hydrogenfluxMobility33(3,3,section) = IO_floatValue(line,positions,2_pInt) case ('vacancy_eqcv') lattice_equilibriumVacancyConcentration(section) = IO_floatValue(line,positions,2_pInt) + case ('hydrogen_eqch') + lattice_equilibriumHydrogenConcentration(section) = IO_floatValue(line,positions,2_pInt) end select endif enddo