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

reworked energy splitting for damage

This commit is contained in:
Pratheek Shanthraj 2014-08-10 10:45:07 +00:00
parent ec71d77038
commit b6ccbc0fe9
3 changed files with 28 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
code/constitutive.f90
View File

@ -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

13
code/damage_gradient.f90
View File

@ -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)

7
code/damage_local.f90
View File

@ -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) = &