added history variable for brittle damage to prevent healing. renamed/removed some lattice damage tensors and symmetrize them.

code/constitutive.f90

@@ -424,6 +424,10 @@ subroutine constitutive_microstructure(temperature, Tstar_v, Fe, Fp, ipc, ip, el
  real(pReal),   intent(in), dimension(3,3) :: &
    Fe, &                                                                                            !< elastic deformation gradient
    Fp                                                                                               !< plastic deformation gradient
+ real(pReal) :: damage, Tstar_v_effective(6)
+ 
+ damage = constitutive_getNonlocalDamage(ipc,ip,el)
+ Tstar_v_effective = Tstar_v/damage*damage
 
  select case (phase_plasticity(material_phase(ipc,ip,el)))
        
@@ -440,7 +444,7 @@ subroutine constitutive_microstructure(temperature, Tstar_v, Fe, Fp, ipc, ip, el
  
  select case (phase_damage(material_phase(ipc,ip,el)))
    case (LOCAL_DAMAGE_BRITTLE_ID)
-     call damage_brittle_microstructure(Tstar_v, Fe, ipc, ip, el)
+     call damage_brittle_microstructure(Tstar_v_effective, Fe, ipc, ip, el)
 
  end select
 
@@ -496,7 +500,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip
  real(pReal) :: damage, Tstar_v_effective(6)
  
  damage = constitutive_getNonlocalDamage(ipc,ip,el)
- Tstar_v_effective = damage*damage*Tstar_v
+ Tstar_v_effective = Tstar_v/damage*damage
  select case (phase_plasticity(material_phase(ipc,ip,el)))
  
    case (PLASTICITY_NONE_ID)

code/damage_brittle.f90

@@ -25,8 +25,9 @@ module damage_brittle
  integer(pInt),                       dimension(:),           allocatable,         private :: &
    damage_brittle_Noutput                                                                   !< number of outputs per instance of this damage 
 
- real(pReal),                         dimension(:),     allocatable,         public :: &
+ real(pReal),                         dimension(:),     allocatable,         private :: &
-   damage_brittle_aTol
+   damage_brittle_aTol, &
+   damage_brittle_critStrainEnergy
 
  enum, bind(c) 
    enumerator :: undefined_ID, &
@@ -117,6 +118,7 @@ subroutine damage_brittle_init(fileUnit)
           damage_brittle_output = ''
  allocate(damage_brittle_outputID(maxval(phase_Noutput),maxNinstance),      source=undefined_ID)
  allocate(damage_brittle_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(damage_brittle_critStrainEnergy(maxNinstance),                    source=0.0_pReal) 
  allocate(damage_brittle_aTol(maxNinstance),                                source=0.0_pReal) 
 
  rewind(fileUnit)
@@ -150,6 +152,9 @@ subroutine damage_brittle_init(fileUnit)
                                                        IO_lc(IO_stringValue(line,positions,2_pInt))
           end select
 
+       case ('critical_strain_energy')
+         damage_brittle_critStrainEnergy(instance) = IO_floatValue(line,positions,2_pInt)
+
        case ('atol_damage')
          damage_brittle_aTol(instance) = IO_floatValue(line,positions,2_pInt)
 
@@ -222,7 +227,7 @@ subroutine damage_brittle_stateInit(phase,instance)
  real(pReal), dimension(damageState(phase)%sizeState) :: tempState
 
  tempState(1) = 1.0_pReal
- tempState(2) = 0.0_pReal
+ tempState(2) = 1.0_pReal
  damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:)))
  damageState(phase)%state0 = damageState(phase)%state
  damageState(phase)%partionedState0 = damageState(phase)%state
@@ -283,16 +288,10 @@ subroutine damage_brittle_microstructure(Tstar_v, Fe, ipc, ip, el)
    phase_damageInstance, &
    damageState
  use math, only: &
-   math_Mandel66to3333, &
    math_Mandel6to33, &
    math_mul33x33, &
-   math_mul3333xx33, &
    math_transpose33, &
-   math_trace33, &
    math_I3
- use lattice, only: &
-   lattice_damageDiffusion33, &
-   lattice_C66
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -311,11 +310,10 @@ subroutine damage_brittle_microstructure(Tstar_v, Fe, ipc, ip, el)
  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)%state(2,constituent) = min(1.0_pReal, &
+ damageState(phase)%state(2,constituent) = min(damageState(phase)%state(2,constituent), &
-                                               math_trace33(lattice_damageDiffusion33(1:3,1:3,phase))/ &
+                                               damage_brittle_critStrainEnergy(phase)/ &
-                                               sum(abs(stress*strain)))
+                                               sum(abs(math_Mandel6to33(Tstar_v)*strain)))
 
 end subroutine damage_brittle_microstructure
 

code/lattice.f90

@@ -741,8 +741,7 @@ module lattice
  real(pReal),                              dimension(:,:,:),   allocatable, public, protected :: &
    lattice_thermalConductivity33, &
    lattice_thermalExpansion33, &
-   lattice_damageDiffusion33, &
+   lattice_damageDiffusion33
-   lattice_surfaceEnergy33
  real(pReal),                              dimension(:),       allocatable, public, protected :: &
    lattice_damageMobility, &
    lattice_massDensity, &
@@ -1004,7 +1003,6 @@ subroutine lattice_init
  allocate(lattice_damageMobility      (    Nphases), source=0.0_pReal)
  allocate(lattice_massDensity         (    Nphases), source=0.0_pReal)
  allocate(lattice_specificHeat        (    Nphases), source=0.0_pReal)
- allocate(lattice_surfaceEnergy33      (3,3,Nphases), source=0.0_pReal)
  allocate(lattice_referenceTemperature     (Nphases), source=0.0_pReal)
 
  allocate(lattice_mu(Nphases),       source=0.0_pReal)
@@ -1110,11 +1108,11 @@ subroutine lattice_init
        aM(section) = IO_floatValue(line,positions,2_pInt)
      case ('cm', 'c_m', 'c_martensite')
        cM(section) = IO_floatValue(line,positions,2_pInt)
-     case ('k11')
+     case ('thermal_conductivity11')
        lattice_thermalConductivity33(1,1,section) = IO_floatValue(line,positions,2_pInt)
-     case ('k22')
+     case ('thermal_conductivity22')
        lattice_thermalConductivity33(2,2,section) = IO_floatValue(line,positions,2_pInt)
-     case ('k33')
+     case ('thermal_conductivity33')
        lattice_thermalConductivity33(3,3,section) = IO_floatValue(line,positions,2_pInt)
      case ('thermal_expansion11')
        lattice_thermalExpansion33(1,1,section) = IO_floatValue(line,positions,2_pInt)
@@ -1122,38 +1120,20 @@ subroutine lattice_init
        lattice_thermalExpansion33(2,2,section) = IO_floatValue(line,positions,2_pInt)
      case ('thermal_expansion33')
        lattice_thermalExpansion33(3,3,section) = IO_floatValue(line,positions,2_pInt)
-     case ('g11')
-       lattice_surfaceEnergy33(1,1,section) = IO_floatValue(line,positions,2_pInt)
-     case ('g22')
-       lattice_surfaceEnergy33(2,2,section) = IO_floatValue(line,positions,2_pInt)
-     case ('g33')
-       lattice_surfaceEnergy33(3,3,section) = IO_floatValue(line,positions,2_pInt)
-     case ('reference_temperature')
-       lattice_referenceTemperature(section) = IO_floatValue(line,positions,2_pInt)
-     case ('k_d11')
-       lattice_DamageDiffusion33(1,1,section) = IO_floatValue(line,positions,2_pInt)
-     case ('k_d12')
-       lattice_DamageDiffusion33(1,2,section) = IO_floatValue(line,positions,2_pInt)
-     case ('k_d13')
-       lattice_DamageDiffusion33(1,3,section) = IO_floatValue(line,positions,2_pInt)
-     case ('k_d21')
-       lattice_DamageDiffusion33(2,1,section) = IO_floatValue(line,positions,2_pInt)
-     case ('k_d22')
-       lattice_DamageDiffusion33(2,3,section) = IO_floatValue(line,positions,2_pInt)
-     case ('k_d23')
-       lattice_DamageDiffusion33(2,3,section) = IO_floatValue(line,positions,2_pInt)
-     case ('k_d31')
-       lattice_DamageDiffusion33(3,1,section) = IO_floatValue(line,positions,2_pInt)
-     case ('k_d32')
-       lattice_DamageDiffusion33(3,2,section) = IO_floatValue(line,positions,2_pInt)
-     case ('k_d33')
-       lattice_DamageDiffusion33(3,3,section) = IO_floatValue(line,positions,2_pInt)
-     case ('damage_mobility')
-       lattice_DamageMobility(section) = IO_floatValue(line,positions,2_pInt)
      case ('specific_heat')
        lattice_specificHeat(section) = IO_floatValue(line,positions,2_pInt)
      case ('mass_density')
        lattice_massDensity(section) = IO_floatValue(line,positions,2_pInt)
+     case ('reference_temperature')
+       lattice_referenceTemperature(section) = IO_floatValue(line,positions,2_pInt)
+     case ('damage_difusion11')
+       lattice_DamageDiffusion33(1,1,section) = IO_floatValue(line,positions,2_pInt)
+     case ('damage_difusion22')
+       lattice_DamageDiffusion33(2,2,section) = IO_floatValue(line,positions,2_pInt)
+     case ('damage_difusion33')
+       lattice_DamageDiffusion33(3,3,section) = IO_floatValue(line,positions,2_pInt)
+     case ('damage_mobility')
+       lattice_DamageMobility(section) = IO_floatValue(line,positions,2_pInt)
      end select
    endif
  enddo
@@ -1241,8 +1221,8 @@ subroutine lattice_initializeStructure(myPhase,CoverA,aA,aM,cM)
                                                                        lattice_thermalConductivity33(1:3,1:3,myPhase))
  lattice_thermalExpansion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
                                                                     lattice_thermalExpansion33(1:3,1:3,myPhase))
- lattice_surfaceEnergy33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
+ lattice_DamageDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
-                                                                 lattice_surfaceEnergy33(1:3,1:3,myPhase))
+                                                                 lattice_DamageDiffusion33(1:3,1:3,myPhase))
  
  select case(lattice_structure(myPhase))
 !--------------------------------------------------------------------------------------------------