From 9779127521b0cd045e9359e5dcac4d15bb105ee2 Mon Sep 17 00:00:00 2001
From: Pratheek Shanthraj
Date: Mon, 20 Oct 2014 15:43:28 +0000
Subject: [PATCH] added intermediate configuration to crystal plasticity
kinematics intended handle intermediate deformation stages consistently in a
finite strain framework
currently implemented for thermal strains, but will expand on this to add damage strains, phase transformation strains etc.
---
code/constitutive.f90 | 37 +++++++++++++++++++++++++++++++------
code/crystallite.f90 | 29 +++++++++++++++++++++--------
2 files changed, 52 insertions(+), 14 deletions(-)
diff --git a/code/constitutive.f90 b/code/constitutive.f90
index 81e738f69..07abb188d 100644
--- a/code/constitutive.f90
+++ b/code/constitutive.f90
@@ -36,6 +36,7 @@ module constitutive
constitutive_getAdiabaticTemperature, &
constitutive_putAdiabaticTemperature, &
constitutive_getTemperature, &
+ constitutive_getThermalStrain, &
constitutive_getLocalVacancyConcentration, &
constitutive_putLocalVacancyConcentration, &
constitutive_getVacancyConcentration, &
@@ -712,10 +713,7 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el)
damage = constitutive_getDamage(ipc,ip,el)
C = damage*damage*math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el))
- T = math_mul3333xx33(C,0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3) - &
- lattice_thermalExpansion33(1:3,1:3,mappingConstitutive(2,ipc,ip,el))* &
- (constitutive_getTemperature(ipc,ip,el) - &
- lattice_referenceTemperature(mappingConstitutive(2,ipc,ip,el))))
+ T = math_mul3333xx33(C,0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3))
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) &
@@ -1169,8 +1167,6 @@ function constitutive_getTemperature(ipc, ip, el)
FIELD_THERMAL_local_ID, &
FIELD_THERMAL_nonlocal_ID, &
material_homog
- use lattice, only: &
- lattice_referenceTemperature
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
@@ -1191,6 +1187,35 @@ function constitutive_getTemperature(ipc, ip, el)
end function constitutive_getTemperature
+!--------------------------------------------------------------------------------------------------
+!> @brief returns thermal deformation gradient
+!--------------------------------------------------------------------------------------------------
+function constitutive_getThermalStrain(ipc, ip, el)
+ use prec, only: &
+ pReal
+ use math, only: &
+ math_I3
+ use lattice, only: &
+ lattice_referenceTemperature, &
+ lattice_thermalExpansion33
+ use material, only: &
+ mappingConstitutive
+
+ implicit none
+ integer(pInt), intent(in) :: &
+ ipc, & !< grain number
+ ip, & !< integration point number
+ el !< element number
+ real(pReal), dimension(3,3) :: &
+ constitutive_getThermalStrain
+
+ constitutive_getThermalStrain = math_I3 + &
+ (constitutive_getTemperature(ipc, ip, el) - &
+ lattice_referenceTemperature(mappingConstitutive(2,ipc,ip,el)))* &
+ lattice_thermalExpansion33(1:3,1:3,mappingConstitutive(2,ipc,ip,el))
+
+end function constitutive_getThermalStrain
+
!--------------------------------------------------------------------------------------------------
!> @brief returns local vacancy concentration
!--------------------------------------------------------------------------------------------------
diff --git a/code/crystallite.f90 b/code/crystallite.f90
index a5bce9854..eb3a11ca1 100644
--- a/code/crystallite.f90
+++ b/code/crystallite.f90
@@ -3296,7 +3296,8 @@ logical function crystallite_integrateStress(&
debug_cumLpTicks, &
debug_StressLoopDistribution
use constitutive, only: constitutive_LpAndItsTangent, &
- constitutive_TandItsTangent
+ constitutive_TandItsTangent, &
+ constitutive_getThermalStrain
use math, only: math_mul33x33, &
math_mul33xx33, &
math_mul66x6, &
@@ -3339,7 +3340,10 @@ logical function crystallite_integrateStress(&
Tstar,& ! 2nd Piola-Kirchhoff Stress
A,&
B, &
- Fe ! elastic deformation gradient
+ Fe, & ! elastic deformation gradient
+ Fi, & ! gradient of intermediate deformation stages
+ Ci, & ! stretch of intermediate deformation stages
+ invFi
real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation
real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK
integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK
@@ -3354,7 +3358,8 @@ logical function crystallite_integrateStress(&
steplength0, &
steplength, &
dt, & ! time increment
- aTol
+ aTol, &
+ detFi
logical error ! flag indicating an error
integer(pInt) NiterationStress, & ! number of stress integrations
ierr, & ! error indicator for LAPACK
@@ -3418,6 +3423,11 @@ logical function crystallite_integrateStress(&
return
endif
A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
+ Fi = math_I3 ! intermediate configuration, assume decomposition as F = Fe Fi Fp
+ Fi = math_mul33x33(Fi,constitutive_getThermalStrain(g,i,e))
+ Ci = math_mul33x33(math_transpose33(Fi),Fi) ! non-plastic stretch tensor (neglecting elastic contributions)
+ invFi = math_inv33(Fi)
+ detFi = math_det33(Fi)
!* start LpLoop with normal step length
@@ -3443,8 +3453,10 @@ logical function crystallite_integrateStress(&
!* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
B = math_I3 - dt*Lpguess
- Fe = math_mul33x33(A,B) ! current elastic deformation tensor
- call constitutive_TandItsTangent(Tstar, dT_dFe3333, Fe, g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative
+ Fe = math_mul33x33(math_mul33x33(A,B),invFi) ! current elastic deformation tensor
+ call constitutive_TandItsTangent(Tstar, dT_dFe3333, Fe, g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration
+ Tstar = math_mul33x33(Ci, math_mul33x33(invFi, &
+ math_mul33x33(Tstar,math_transpose33(invFi))))/detFi ! push Tstar forward to plastic (lattice) configuration
Tstar_v = math_Mandel33to6(Tstar)
@@ -3574,13 +3586,14 @@ logical function crystallite_integrateStress(&
#endif
return
endif
- Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe
+ Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi) ! calc resulting Fe
!* calculate 1st Piola-Kirchhoff stress
- crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_new, math_mul33x33(math_Mandel6to33(Tstar_v), &
- math_transpose33(invFp_new)))
+ crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(Fe_new,Fi), &
+ math_mul33x33(math_Mandel6to33(Tstar_v), &
+ math_transpose33(invFp_new)))*detFi
!* store local values in global variables