From b6ccbc0fe9029b89a1b7edb481eaaff01bcd797a Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Sun, 10 Aug 2014 10:45:07 +0000 Subject: [PATCH] reworked energy splitting for damage --- code/constitutive.f90 | 35 ++++++++++++++++++----------------- code/damage_gradient.f90 | 13 ++++++------- code/damage_local.f90 | 7 ++++--- 3 files changed, 28 insertions(+), 27 deletions(-) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 7a5a1ab38..085dadffb 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -480,30 +480,31 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el) dT_dFe !< dT/dFe integer(pInt) :: i, j, k, l - real(pReal) :: damage, negative_volStrain, negative_volStress, Csum - real(pReal), dimension(3,3) :: strain + real(pReal) :: damage, pressure real(pReal), dimension(3,3,3,3) :: C C = 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_temperature(ipc,ip,el) - & + lattice_referenceTemperature(mappingConstitutive(2,ipc,ip,el)))) + + pressure = math_trace33(T)/3.0_pReal 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 + T = damage*T + if (pressure < 0.0_pReal) T = T + (1.0_pReal - damage)*pressure*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) & - dT_dFe(i,j,k,l) = math_mul3x3(C(i,j,l,1:3),Fe(k,1:3)) ! dT*_ij/dFe_kl - + dT_dFe(i,j,k,l) = sum(C(i,j,l,1:3)*Fe(k,1:3)) ! dT*_ij/dFe_kl + + if (pressure < 0.0_pReal) then + do i=1_pInt, 3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt; do j=1_pInt,3_pInt + dT_dFe(i,i,k,l) = dT_dFe(i,i,k,l) + & + (1.0_pReal - damage)*math_trace33(C(1:3,1:3,l,j))*Fe(k,j)/3.0_pReal + enddo; enddo; enddo; enddo + endif + end subroutine constitutive_hooke_TandItsTangent diff --git a/code/damage_gradient.f90 b/code/damage_gradient.f90 index c35417d02..72f3bbd8c 100644 --- a/code/damage_gradient.f90 +++ b/code/damage_gradient.f90 @@ -286,18 +286,17 @@ subroutine damage_gradient_microstructure(Tstar_v, Fe, ipc, ip, el) integer(pInt) :: & phase, constituent real(pReal) :: & - strainEnergy, strain(3,3), negative_volStrain + strainEnergy, pressure, stress(3,3), strain(3,3) phase = mappingConstitutive(2,ipc,ip,el) constituent = mappingConstitutive(1,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) - - 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)/ & + stress = math_mul3333xx33(math_Mandel66to3333(lattice_C66(1:6,1:6,phase)),strain) + pressure = math_trace33(stress)/3.0_pReal + if (pressure < 0.0_pReal) stress = stress - pressure*math_I3 + strainEnergy = 0.5*sum(abs(stress*strain)) + damageState(phase)%state(2,constituent) = strainEnergy/ & (math_trace33(lattice_surfaceEnergy33(1:3,1:3,phase))/3.0_pReal) + & damageState(phase)%state(1,constituent) diff --git a/code/damage_local.f90 b/code/damage_local.f90 index f51283734..72f528298 100644 --- a/code/damage_local.f90 +++ b/code/damage_local.f90 @@ -282,18 +282,19 @@ subroutine damage_local_dotState(Tstar_v, Fe, Lp, ipc, ip, el) integer(pInt) :: & phase, constituent, instance real(pReal) :: & - trialDamage, strain(3,3), stress(3,3), negative_volStrain + trialDamage, strain(3,3), stress(3,3), pressure 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) stress = math_mul3333xx33(math_Mandel66to3333(lattice_C66(1:6,1:6,phase)),strain) + pressure = math_trace33(stress)/3.0_pReal + if (pressure < 0.0_pReal) stress = stress - pressure*math_I3 trialDamage = min(1.0_pReal, & (math_trace33(lattice_surfaceEnergy33(1:3,1:3,phase))/3.0_pReal)/ & - (abs(sum((strain - negative_volStrain*math_I3)*stress)) + & + (abs(sum(strain*stress)) + & damageState(phase)%state(1,constituent))) damageState(phase)%dotState(1,constituent) = &