initialise Fi correctly for initial field values away from equilibrium

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

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.

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  
 !--------------------------------------------------------------------------------------------------

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  

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  
 !--------------------------------------------------------------------------------------------------

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