From c4b9ddf45da8f329ccc9aebdcf34b51a6e180650 Mon Sep 17 00:00:00 2001
From: Luv	Sharma <l.sharma@mpie.de>
Date: Fri, 8 Aug 2014 15:37:32 +0000
Subject: [PATCH] added strain decomposition to uncouple damage from negative
 hydrostatic part

---
 code/constitutive.f90    | 36 +++++++++++++++++++++++++++++-------
 code/damage_gradient.f90 | 12 ++++++++++--
 2 files changed, 39 insertions(+), 9 deletions(-)

diff --git a/code/constitutive.f90 b/code/constitutive.f90
index 87220667a..2dc684855 100644
--- a/code/constitutive.f90
+++ b/code/constitutive.f90
@@ -455,7 +455,8 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el)
    math_mul3333xx33, &
    math_Mandel66to3333, &
    math_transpose33, &
-   MATH_I3
+   MATH_I3, &
+   math_trace33
  use material, only: &
    mappingConstitutive
  use lattice, only: &
@@ -479,21 +480,42 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el)
    dT_dFe                                                                                           !< dT/dFe
  
  integer(pInt) :: i, j, k, l
- real(pReal), dimension(3,3)     :: FeT
- real(pReal), dimension(3,3,3,3) :: C
+ real(pReal), dimension(3,3)     :: FeT,strain, CxxDel, CxxDel_undamaged
+ real(pReal), dimension(3,3,3,3) :: C, C_undamged
+ real(pReal) :: strain_trace
 
- C = constitutive_damageValue(ipc,ip,el)*&
+ C_undamged = math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el))
+ C  = constitutive_damageValue(ipc,ip,el)*&
      math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el))
 
  FeT = math_transpose33(Fe)
- T = math_mul3333xx33(C,0.5_pReal*(math_mul33x33(FeT,Fe)-MATH_I3) - &
+ strain = 0.5_pReal*(math_mul33x33(FeT,Fe)-MATH_I3)
+ strain_trace = math_trace33(strain)
+
+
+ if(strain_trace>=0.0_pReal) then
+ T = math_mul3333xx33(C,strain) - &
                         lattice_thermalExpansion33(1:3,1:3,mappingConstitutive(2,ipc,ip,el))* &
                         (constitutive_temperature(ipc,ip,el) - &
-                         lattice_referenceTemperature(mappingConstitutive(2,ipc,ip,el))))
+                         lattice_referenceTemperature(mappingConstitutive(2,ipc,ip,el)))
  dT_dFe = 0.0_pReal
  forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt, k=1_pInt:3_pInt, l=1_pInt:3_pInt) &
    dT_dFe(i,j,k,l) = math_mul3x3(C(i,j,l,1:3),Fe(k,1:3))                                            ! dT*_ij/dFe_kl
-
+ else  
+  strain = strain - (strain_trace/3)*math_I3                                                        ! removing (negative) volumetric part
+  
+ T = math_mul3333xx33(C,strain) + &
+           math_mul3333xx33(C_undamged,((strain_trace/3)*math_I3 )) - &
+                        lattice_thermalExpansion33(1:3,1:3,mappingConstitutive(2,ipc,ip,el))* &
+                        (constitutive_temperature(ipc,ip,el) - &
+                         lattice_referenceTemperature(mappingConstitutive(2,ipc,ip,el)))
+  CxxDel           = math_mul3333xx33(C,math_I3)                                                    ! temporary variable C_ijxy * del_xy (damaged part)
+  CxxDel_undamaged = math_mul3333xx33(C_undamged,math_I3)                                                    ! temporary variable C_ijxy * del_xy (Undamaged part)
+  forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt, k=1_pInt:3_pInt, l=1_pInt:3_pInt) &
+    dT_dFe(i,j,k,l) = math_mul3x3(C(i,j,l,1:3),Fe(k,1:3)) - &
+                              0.5*CxxDel(i,j)*Fe(k,l) + &
+                               0.5*CxxDel_undamaged(i,j)*Fe(k,l)                                    ! dT*_ij/dFe_kl = C_ijlt*Fe_kt - 0.5*Cxxdel_ij*Fe_kl
+ endif
 end subroutine constitutive_hooke_TandItsTangent
 
 
diff --git a/code/damage_gradient.f90 b/code/damage_gradient.f90
index 33789591f..b731553c7 100644
--- a/code/damage_gradient.f90
+++ b/code/damage_gradient.f90
@@ -287,14 +287,22 @@ subroutine damage_gradient_microstructure(Tstar_v, Fe, ipc, ip, el)
  integer(pInt) :: &
    phase, constituent 
  real(pReal) :: &
-   strainEnergy, strain(3,3)
+   strainEnergy, strain(3,3),strain_trace
 
  phase = mappingConstitutive(2,ipc,ip,el)
  constituent = mappingConstitutive(1,ipc,ip,el)
 
  strain = 0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
- strainEnergy = 0.5_pReal*sum(strain*math_mul3333xx33(math_Mandel66to3333(lattice_C66(1:6,1:6,phase)), &
+ strain_trace = math_trace33(strain)
+
+ if(strain_trace < 0.0_pReal) then
+  strain = strain  - (strain_trace/3.0_pReal)*math_I3
+  strainEnergy = 0.5_pReal*sum(strain*math_mul3333xx33(math_Mandel66to3333(lattice_C66(1:6,1:6,phase)), &
+                                                      strain ))
+ else
+  strainEnergy = 0.5_pReal*sum(strain*math_mul3333xx33(math_Mandel66to3333(lattice_C66(1:6,1:6,phase)), &
                                                       strain))
+ endif
 
  damageState(phase)%state(2,constituent) = abs(strainEnergy)/ &
                                            (math_trace33(lattice_surfaceEnergy33(1:3,1:3,phase))/3.0_pReal) + &