diff --git a/code/constitutive_damage.f90 b/code/constitutive_damage.f90 index 46f6d7e5d..d55c65bf1 100644 --- a/code/constitutive_damage.f90 +++ b/code/constitutive_damage.f90 @@ -197,7 +197,7 @@ subroutine constitutive_damage_collectDotState(Tstar_v, Fe, Lp, ipc, ip, el) call damage_local_dotState(Tstar_v, Fe, Lp, ipc, ip, el) case (DAMAGE_gradient_ID) - call damage_gradient_dotState(Tstar_v, Lp, ipc, ip, el) + call damage_gradient_dotState(Tstar_v, Fe, Lp, ipc, ip, el) end select diff --git a/code/damage_gradient.f90 b/code/damage_gradient.f90 index b731553c7..c35417d02 100644 --- a/code/damage_gradient.f90 +++ b/code/damage_gradient.f90 @@ -262,7 +262,6 @@ end subroutine damage_gradient_aTolState subroutine damage_gradient_microstructure(Tstar_v, Fe, ipc, ip, el) use material, only: & mappingConstitutive, & - phase_damageInstance, & damageState use math, only: & math_Mandel66to3333, & @@ -287,22 +286,16 @@ subroutine damage_gradient_microstructure(Tstar_v, Fe, ipc, ip, el) integer(pInt) :: & phase, constituent real(pReal) :: & - strainEnergy, strain(3,3),strain_trace + strainEnergy, strain(3,3), negative_volStrain phase = mappingConstitutive(2,ipc,ip,el) constituent = mappingConstitutive(1,ipc,ip,el) strain = 0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3) - strain_trace = math_trace33(strain) + negative_volStrain = min(0.0_pReal,math_trace33(strain)/3.0_pReal) - 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 + strainEnergy = 0.5_pReal*sum(strain*math_mul3333xx33(math_Mandel66to3333(lattice_C66(1:6,1:6,phase)), & + strain - negative_volStrain*math_I3)) damageState(phase)%state(2,constituent) = abs(strainEnergy)/ & (math_trace33(lattice_surfaceEnergy33(1:3,1:3,phase))/3.0_pReal) + & @@ -313,14 +306,19 @@ end subroutine damage_gradient_microstructure !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine damage_gradient_dotState(Tstar_v, Lp, ipc, ip, el) +subroutine damage_gradient_dotState(Tstar_v, Fe, Lp, ipc, ip, el) use material, only: & mappingConstitutive, & damageState use math, only: & - math_Mandel6to33, & - math_trace33 + math_Mandel66to3333, & + math_mul33x33, & + math_mul3333xx33, & + math_transpose33, & + math_trace33, & + math_I3 use lattice, only: & + lattice_C66, & lattice_surfaceEnergy33 implicit none @@ -331,15 +329,20 @@ subroutine damage_gradient_dotState(Tstar_v, Lp, ipc, ip, el) real(pReal), intent(in), dimension(6) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) real(pReal), intent(in), dimension(3,3) :: & - Lp + Fe, Lp integer(pInt) :: & - instance, phase, constituent + phase, constituent + real(pReal), dimension(3,3) :: & + stress, strain phase = mappingConstitutive(2,ipc,ip,el) constituent = mappingConstitutive(1,ipc,ip,el) + strain = 0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3) + stress = math_mul3333xx33(math_Mandel66to3333(lattice_C66(1:6,1:6,phase)), strain) + damageState(phase)%dotState(1,constituent) = & - sum(abs(math_Mandel6to33(Tstar_v)*Lp))/ & + sum(abs(stress*Lp))/ & (math_trace33(lattice_surfaceEnergy33(1:3,1:3,phase))/3.0_pReal) end subroutine damage_gradient_dotState diff --git a/code/damage_local.f90 b/code/damage_local.f90 index d0ab00f74..b57ebd3f8 100644 --- a/code/damage_local.f90 +++ b/code/damage_local.f90 @@ -255,7 +255,6 @@ 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, & @@ -280,22 +279,23 @@ subroutine damage_local_dotState(Tstar_v, Fe, Lp, ipc, ip, el) Lp, & Fe integer(pInt) :: & - phase, constituent, instance + phase, constituent real(pReal) :: & - trialDamage, strain(3,3) + 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) strain = 0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3) + negative_volStrain = min(0.0_pReal,math_trace33(strain)/3.0_pReal) + stress = math_mul3333xx33(math_Mandel66to3333(lattice_C66(1:6,1:6,phase)),strain) trialDamage = min(1.0_pReal, & (math_trace33(lattice_surfaceEnergy33(1:3,1:3,phase))/3.0_pReal)/ & - (abs(sum(strain*math_mul3333xx33(math_Mandel66to3333(lattice_C66(1:6,1:6,phase)),strain))) + & + (abs(sum((strain - negative_volStrain*math_I3)*stress)) + & damageState(phase)%state(1,constituent))) damageState(phase)%dotState(1,constituent) = & - sum(abs(math_Mandel6to33(Tstar_v)*Lp)) + sum(abs(stress*Lp)) damageState(phase)%dotState(2,constituent) = & damage_local_crack_mobility(instance)* & (trialDamage - damageState(phase)%state(2,constituent))