223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

cleaned up strain energy splitting. coupled damage to plasticity for J2 and phenopowerlaw

This commit is contained in:
Pratheek Shanthraj 2014-08-08 22:29:38 +00:00
parent 9c40f187a0
commit 3700e132b0
4 changed files with 59 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
code/constitutive.f90
View File

@ -397,9 +397,9 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip
Lp = 0.0_pReal
dLp_dTstar = math_identity2nd(9)
case (PLASTICITY_J2_ID)
call constitutive_j2_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,ipc,ip,el)
call constitutive_j2_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,constitutive_damageValue(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID)
call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,ipc,ip,el)
call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,constitutive_damageValue(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID)
call constitutive_nonlocal_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,temperature, ip,el)
case (PLASTICITY_DISLOTWIN_ID)
@ -455,8 +455,8 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el)
math_mul3333xx33, &
math_Mandel66to3333, &
math_transpose33, &
MATH_I3, &
math_trace33
math_trace33, &
math_I3
use material, only: &
mappingConstitutive
use lattice, only: &
@ -480,42 +480,30 @@ 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,strain, CxxDel, CxxDel_undamaged
real(pReal), dimension(3,3,3,3) :: C, C_undamged
real(pReal) :: strain_trace
real(pReal) :: damage, negative_volStrain, negative_volStress, Csum
real(pReal), dimension(3,3) :: strain
real(pReal), dimension(3,3,3,3) :: C
C_undamged = math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el))
C = constitutive_damageValue(ipc,ip,el)*&
math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el))
C = math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el))
damage = constitutive_damageValue(ipc,ip,el)
strain = 0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
negative_volStrain = min(0.0_pReal,math_trace33(strain)/3.0_pReal)
negative_volStress = math_trace33(math_mul3333xx33(C,negative_volStrain*math_I3))/3.0_pReal
T = damage* &
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)))) + &
(1.0_pReal - damage)*negative_volStress*math_I3
Csum = sum(math_I3*math_mul3333xx33(C,math_I3))/9.0_pReal
C = damage*C
forall (i = 1_pInt:3_pInt) &
C(1:3,1:3,i,i) = C(1:3,1:3,i,i) + (1.0_pReal - damage)*Csum*math_I3
FeT = math_transpose33(Fe)
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)))
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
@ -549,6 +537,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature,
PLASTICITY_DISLOKMC_ID, &
PLASTICITY_TITANMOD_ID, &
PLASTICITY_NONLOCAL_ID
use constitutive_damage, only: &
constitutive_damageValue
use constitutive_j2, only: &
constitutive_j2_dotState
use constitutive_phenopowerlaw, only: &
@ -587,9 +577,9 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature,
select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_J2_ID)
call constitutive_j2_dotState (Tstar_v,ipc,ip,el)
call constitutive_j2_dotState (Tstar_v,constitutive_damageValue(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID)
call constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
call constitutive_phenopowerlaw_dotState(Tstar_v,constitutive_damageValue(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID)
call constitutive_dislotwin_dotState (Tstar_v,Temperature,ipc,ip,el)
case (PLASTICITY_DISLOKMC_ID)
@ -697,6 +687,8 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
PLASTICITY_DISLOKMC_ID, &
PLASTICITY_TITANMOD_ID, &
PLASTICITY_NONLOCAL_ID
use constitutive_damage, only: &
constitutive_damageValue
use constitutive_j2, only: &
#ifdef HDF
constitutive_j2_postResults2,&
@ -734,7 +726,7 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
case (PLASTICITY_J2_ID)
constitutive_postResults= constitutive_j2_postResults (Tstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID)
constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,constitutive_damageValue(ipc,ip,el),ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID)
constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
case (PLASTICITY_DISLOKMC_ID)

12
code/constitutive_j2.f90
View File

@ -343,7 +343,7 @@ end subroutine constitutive_j2_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,damage,ipc,ip,el)
use math, only: &
math_mul6x6, &
math_Mandel6to33, &
@ -368,6 +368,8 @@ subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
damage
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
@ -395,7 +397,7 @@ subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
dLp_dTstar99 = 0.0_pReal
else
gamma_dot = constitutive_j2_gdot0(instance) &
* (sqrt(1.5_pReal) * norm_Tstar_dev / (constitutive_j2_fTaylor(instance) * &
* (sqrt(1.5_pReal) * norm_Tstar_dev / (damage * constitutive_j2_fTaylor(instance) * &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el)))) &
**constitutive_j2_n(instance)
@ -419,7 +421,7 @@ end subroutine constitutive_j2_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine constitutive_j2_dotState(Tstar_v,ipc,ip,el)
subroutine constitutive_j2_dotState(Tstar_v,damage,ipc,ip,el)
use math, only: &
math_mul6x6
use mesh, only: &
@ -441,6 +443,8 @@ subroutine constitutive_j2_dotState(Tstar_v,ipc,ip,el)
el !< element
real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
damage
real(pReal) :: &
gamma_dot, & !< strainrate
hardening, & !< hardening coefficient
@ -465,7 +469,7 @@ subroutine constitutive_j2_dotState(Tstar_v,ipc,ip,el)
! strain rate
gamma_dot = constitutive_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!-----------------------------------------------------------------------------------
(constitutive_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**constitutive_j2_n(instance)
(damage*constitutive_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**constitutive_j2_n(instance)
!--------------------------------------------------------------------------------------------------
! hardening coefficient

30
code/constitutive_phenopowerlaw.f90
View File

@ -662,7 +662,7 @@ end subroutine constitutive_phenopowerlaw_aTolState
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,damage,ipc,ip,el)
use prec, only: &
p_vec
use math, only: &
@ -696,6 +696,8 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ip
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
damage
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
@ -753,11 +755,11 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ip
lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph)
enddo
gdot_slip_pos(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* &
((abs(tau_slip_pos(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))*&
((abs(tau_slip_pos(j))/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance))*&
sign(1.0_pReal,tau_slip_pos(j))
gdot_slip_neg(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* &
((abs(tau_slip_neg(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))*&
((abs(tau_slip_neg(j))/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance))*&
sign(1.0_pReal,tau_slip_neg(j))
Lp = Lp + (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
@ -795,7 +797,7 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ip
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph))
gdot_twin(j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau_twin(j))/plasticState(ph)%state(nSlip+j,of))**&
(abs(tau_twin(j))/(damage*plasticState(ph)%state(nSlip+j,of)))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
Lp = Lp + gdot_twin(j)*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
@ -819,7 +821,7 @@ end subroutine constitutive_phenopowerlaw_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
subroutine constitutive_phenopowerlaw_dotState(Tstar_v,damage,ipc,ip,el)
use lattice, only: &
lattice_Sslip_v, &
lattice_Stwin_v, &
@ -842,6 +844,8 @@ subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
damage
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
@ -916,8 +920,8 @@ subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
enddo
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))&
((abs(tau_slip_pos(j))/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg(j))/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos(j))
enddo
enddo slipFamiliesLoop1
@ -938,7 +942,7 @@ subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph))
gdot_twin(j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau_twin(j))/plasticState(ph)%state(nslip+j,of))**&
(abs(tau_twin(j))/(damage*plasticState(ph)%state(nslip+j,of)))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
enddo
enddo twinFamiliesLoop1
@ -988,7 +992,7 @@ end subroutine constitutive_phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function constitutive_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
function constitutive_phenopowerlaw_postResults(Tstar_v,damage,ipc,ip,el)
use prec, only: &
p_vec
use mesh, only: &
@ -1016,6 +1020,8 @@ function constitutive_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
damage
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
@ -1073,8 +1079,8 @@ function constitutive_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
enddo
constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos)/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg)/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))&
((abs(tau_slip_pos)/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg)/(damage*plasticState(ph)%state(j,of)))**constitutive_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos)
enddo
@ -1116,7 +1122,7 @@ function constitutive_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph))
constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau)/plasticState(ph)%state(j+nSlip,of))**&
(abs(tau)/(damage*plasticState(ph)%state(j+nSlip,of)))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau))
enddo
enddo twinFamiliesLoop1

4
code/damage_local.f90
View File

@ -255,6 +255,7 @@ end subroutine damage_local_aTolState
subroutine damage_local_dotState(Tstar_v, Fe, Lp, ipc, ip, el)
use material, only: &
mappingConstitutive, &
phase_damageInstance, &
damageState
use math, only: &
math_Mandel66to3333, &
@ -279,12 +280,13 @@ subroutine damage_local_dotState(Tstar_v, Fe, Lp, ipc, ip, el)
Lp, &
Fe
integer(pInt) :: &
phase, constituent
phase, constituent, instance
real(pReal) :: &
trialDamage, strain(3,3), stress(3,3), negative_volStrain
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase) ! which instance of my damage is present phase
strain = 0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
negative_volStrain = min(0.0_pReal,math_trace33(strain)/3.0_pReal)