From 3d94ae4074d291bd4e9252fa80fb70c55335e968 Mon Sep 17 00:00:00 2001
From: Pratheek	Shanthraj <p.shanthraj@mpie.de>
Date: Wed, 5 Aug 2015 14:12:54 +0000
Subject: [PATCH] finite strain correction for respective expansions

---
 code/kinematics_hydrogen_strain.f90   | 14 +++++++++++---
 code/kinematics_thermal_expansion.f90 | 13 +++++++++++--
 code/kinematics_vacancy_strain.f90    | 16 +++++++++++++---
 3 files changed, 35 insertions(+), 8 deletions(-)

diff --git a/code/kinematics_hydrogen_strain.f90 b/code/kinematics_hydrogen_strain.f90
index 47e6d7724..6ee6b6124 100644
--- a/code/kinematics_hydrogen_strain.f90
+++ b/code/kinematics_hydrogen_strain.f90
@@ -187,10 +187,13 @@ subroutine kinematics_hydrogen_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, i
  use material, only: &
    material_phase, &
    material_homog, &
+   hydrogenConc, &
    hydrogenConcRate, &
    hydrogenfluxMapping
  use math, only: &
    math_I3
+ use lattice, only: &
+   lattice_equilibriumHydrogenConcentration
  
  implicit none
  integer(pInt), intent(in) :: &
@@ -205,15 +208,20 @@ subroutine kinematics_hydrogen_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, i
    phase, &
    instance, &
    homog, offset
+ real(pReal) :: &
+   Ch, ChEq, ChDot  
    
  phase = material_phase(ipc,ip,el)
  instance = kinematics_hydrogen_strain_instance(phase)
  homog = material_homog(ip,el)
  offset = hydrogenfluxMapping(homog)%p(ip,el)
+ Ch = hydrogenConc(homog)%p(offset)
+ ChDot = hydrogenConcRate(homog)%p(offset)
+ ChEq = lattice_equilibriumHydrogenConcentration(phase)
  
- Li = hydrogenConcRate(homog)%p(offset)* &
-      kinematics_hydrogen_strain_coeff(instance)* &
-      math_I3
+ Li = ChDot*math_I3* &
+      kinematics_hydrogen_strain_coeff(instance)/ &
+      (1.0_pReal + kinematics_hydrogen_strain_coeff(instance)*(Ch - ChEq))
  dLi_dTstar3333 = 0.0_pReal
   
 end subroutine kinematics_hydrogen_strain_LiAndItsTangent
diff --git a/code/kinematics_thermal_expansion.f90 b/code/kinematics_thermal_expansion.f90
index b4119a31e..49adf17fd 100644
--- a/code/kinematics_thermal_expansion.f90
+++ b/code/kinematics_thermal_expansion.f90
@@ -191,10 +191,12 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc,
  use material, only: &
    material_phase, &
    material_homog, &
+   temperature, &
    temperatureRate, &
    thermalMapping
  use lattice, only: &
-   lattice_thermalExpansion33
+   lattice_thermalExpansion33, &
+   lattice_referenceTemperature
  
  implicit none
  integer(pInt), intent(in) :: &
@@ -208,12 +210,19 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc,
  integer(pInt) :: &
    phase, &
    homog, offset
+ real(pReal) :: &
+   T, TRef, TDot  
    
  phase = material_phase(ipc,ip,el)
  homog = material_homog(ip,el)
  offset = thermalMapping(homog)%p(ip,el)
+ T = temperature(homog)%p(offset)
+ TDot = temperatureRate(homog)%p(offset)
+ TRef = lattice_referenceTemperature(phase)
  
- Li = temperatureRate(homog)%p(offset) * lattice_thermalExpansion33(1:3,1:3,phase)
+ Li = TDot* &
+      lattice_thermalExpansion33(1:3,1:3,phase)/ &
+      (1.0_pReal + lattice_thermalExpansion33(1:3,1:3,phase)*(T - TRef))
  dLi_dTstar3333 = 0.0_pReal 
   
 end subroutine kinematics_thermal_expansion_LiAndItsTangent
diff --git a/code/kinematics_vacancy_strain.f90 b/code/kinematics_vacancy_strain.f90
index f013acee4..e0de4532d 100644
--- a/code/kinematics_vacancy_strain.f90
+++ b/code/kinematics_vacancy_strain.f90
@@ -186,10 +186,13 @@ subroutine kinematics_vacancy_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip
  use material, only: &
    material_phase, &
    material_homog, &
+   vacancyConc, &
    vacancyConcRate, &
    vacancyfluxMapping
  use math, only: &
    math_I3
+ use lattice, only: &
+   lattice_equilibriumVacancyConcentration
  
  implicit none
  integer(pInt), intent(in) :: &
@@ -204,15 +207,22 @@ subroutine kinematics_vacancy_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip
    phase, &
    instance, &
    homog, offset
+ real(pReal) :: &
+   Cv, CvEq, CvDot  
    
  phase = material_phase(ipc,ip,el)
  instance = kinematics_vacancy_strain_instance(phase)
  homog = material_homog(ip,el)
  offset = vacancyfluxMapping(homog)%p(ip,el)
  
- Li = vacancyConcRate(homog)%p(offset)* &
-      kinematics_vacancy_strain_coeff(instance)* &
-      math_I3
+ Cv = vacancyConc(homog)%p(offset)
+ CvDot = vacancyConcRate(homog)%p(offset)
+ CvEq = lattice_equilibriumvacancyConcentration(phase)
+ 
+ Li = CvDot*math_I3* &
+      kinematics_vacancy_strain_coeff(instance)/ &
+      (1.0_pReal + kinematics_vacancy_strain_coeff(instance)*(Cv - CvEq))
+      
  dLi_dTstar3333 = 0.0_pReal
   
 end subroutine kinematics_vacancy_strain_LiAndItsTangent