thermodynamically consistent treatment of nonlocality over all damage models. switched from power law to exponential with viscous drag type rate formulation for damage evolution

code/damage_anisoBrittle.f90

@@ -348,7 +348,7 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el)
  real(pReal) :: &
    traction_d, traction_t, traction_n, traction_crit, &
    udotd, udott, udotn, &
-   nonlocalFactor, localDamage
+   drivingForce, nonlocalFactor, localDamage
 
  phase = mappingConstitutive(2,ipc,ip,el)
  constituent = mappingConstitutive(1,ipc,ip,el)
@@ -368,8 +368,8 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el)
      traction_t    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
      traction_n    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
      traction_crit = damage_anisoBrittle_critLoad(f,instance)* &
-                     damageState(phase)%state0(index_d,constituent)* &
+                     damageState(phase)%state(index_d,constituent)* &
-                     damageState(phase)%state0(index_d,constituent)
+                     damageState(phase)%state(index_d,constituent)
                     
      udotd = &
        damage_anisoBrittle_sdot_0(instance)* &
@@ -384,11 +384,13 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el)
      damageState(phase)%dotState(index_o,constituent) = &
        (udotd + udott + udotn)/damage_anisoBrittle_critDisp(f,instance)
 
+     drivingForce = max(0.0_pReal, &
+                        2.0_pReal*damageState(phase)%state(index_d,constituent)*  &
+                        damageState(phase)%state(index_o,constituent) - &
+                        nonlocalFactor)
      damageState(phase)%dotState(index_d,constituent) = &
-      -damage_anisoBrittle_sdot_0(instance)* &
+       (exp(-drivingForce**damage_anisoBrittle_N(instance)) - 1.0_pReal)/ &
-       max(0.0_pReal,2.0_pReal*damageState(phase)%state(index_d,constituent)*  &
+       lattice_DamageMobility(phase)
-                     damageState(phase)%state(index_o,constituent) - &
-                     nonlocalFactor)**damage_anisoBrittle_N(instance)
 
      index_d = index_d + 1_pInt; index_o = index_o + 1_pInt
    enddo
@@ -450,8 +452,8 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
      traction_t    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
      traction_n    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
      traction_crit = damage_anisoBrittle_critLoad(f,instance)* &
-                     damageState(phase)%state0(index,constituent)* &
+                     damageState(phase)%state(index,constituent)* &
-                     damageState(phase)%state0(index,constituent)
+                     damageState(phase)%state(index,constituent)
      udotd = &
        sign(1.0_pReal,traction_d)* &
        damage_anisoBrittle_sdot_0(instance)* &
@@ -581,8 +583,8 @@ subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el)
      traction_t    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
      traction_n    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
      traction_crit = damage_anisoBrittle_critLoad(f,instance)* &
-                     damageState(phase)%state0(index,constituent)* &
+                     damageState(phase)%state(index,constituent)* &
-                     damageState(phase)%state0(index,constituent)
+                     damageState(phase)%state(index,constituent)
      udotd = &
        sign(1.0_pReal,traction_d)* &
        damage_anisoBrittle_sdot_0(instance)* &

code/damage_anisoDuctile.f90

@@ -33,7 +33,6 @@ module damage_anisoDuctile
    
  real(pReal),                         dimension(:),           allocatable,         private :: &
    damage_anisoDuctile_aTol_damage, &
-   damage_anisoDuctile_sdot_0, &
    damage_anisoDuctile_N
 
  real(pReal),                         dimension(:,:),         allocatable,         private :: &
@@ -140,7 +139,6 @@ subroutine damage_anisoDuctile_init(fileUnit)
  allocate(damage_anisoDuctile_Nslip(lattice_maxNslipFamily,maxNinstance),        source=0_pInt)
  allocate(damage_anisoDuctile_totalNslip(maxNinstance),                          source=0_pInt)
  allocate(damage_anisoDuctile_aTol_damage(maxNinstance),                         source=0.0_pReal) 
- allocate(damage_anisoDuctile_sdot_0(maxNinstance),                              source=0.0_pReal) 
  allocate(damage_anisoDuctile_N(maxNinstance),                                   source=0.0_pReal) 
 
  rewind(fileUnit)
@@ -177,9 +175,6 @@ subroutine damage_anisoDuctile_init(fileUnit)
        case ('atol_damage')
          damage_anisoDuctile_aTol_damage(instance) = IO_floatValue(line,positions,2_pInt)
          
-       case ('sdot_0')
-         damage_anisoDuctile_sdot_0(instance) = IO_floatValue(line,positions,2_pInt)
-         
        case ('rate_sensitivity_damage')
          damage_anisoDuctile_N(instance) = IO_floatValue(line,positions,2_pInt)
          
@@ -315,6 +310,7 @@ subroutine damage_anisoDuctile_dotState(nSlip, accumulatedSlip, ipc, ip, el)
    index, f, i
  real(pReal) :: &
    localDamage, &
+   drivingForce, &
    nonlocalFactor
 
  phase = mappingConstitutive(2,ipc,ip,el)
@@ -329,12 +325,13 @@ subroutine damage_anisoDuctile_dotState(nSlip, accumulatedSlip, ipc, ip, el)
  index = 1_pInt
  do f = 1_pInt,lattice_maxNslipFamily
    do i = 1_pInt,damage_anisoDuctile_Nslip(f,instance)                                            ! process each (active) slip system in family
+     drivingForce = max(0.0_pReal, &
+                        2.0_pReal*damageState(phase)%state(index+1,constituent)* &
+                        accumulatedSlip(index)/damage_anisoDuctile_critAccShear(f,instance) - &
+                        nonlocalFactor)
      damageState(phase)%dotState(index+1,constituent) = &
-      -damage_anisoDuctile_sdot_0(instance)* &
+       (exp(-drivingForce**damage_anisoDuctile_N(instance)) - 1.0_pReal)/ &
-       max(0.0_pReal, &
+       lattice_DamageMobility(phase)
-           2.0_pReal*damageState(phase)%state(index+1,constituent)* &
-           accumulatedSlip(index)/damage_anisoDuctile_critAccShear(f,instance) - &
-           nonlocalFactor)**damage_anisoDuctile_N(instance)
      index = index + 1_pInt
    enddo
  enddo

code/damage_isoBrittle.f90

@@ -27,6 +27,7 @@ module damage_isoBrittle
 
  real(pReal),                         dimension(:),     allocatable,         private :: &
    damage_isoBrittle_aTol, &
+   damage_isoBrittle_N, &
    damage_isoBrittle_critStrainEnergy
 
  enum, bind(c) 
@@ -43,11 +44,11 @@ module damage_isoBrittle
    damage_isoBrittle_stateInit, &
    damage_isoBrittle_aTolState, &
    damage_isoBrittle_dotState, &
-   damage_isoBrittle_microstructure, &
    damage_isoBrittle_getDamage, &
    damage_isoBrittle_putLocalDamage, &
    damage_isoBrittle_getLocalDamage, &
    damage_isoBrittle_getDamageDiffusion33, &
+   damage_isoBrittle_getDamagedC66, &
    damage_isoBrittle_postResults
 
 contains
@@ -124,6 +125,7 @@ subroutine damage_isoBrittle_init(fileUnit)
           damage_isoBrittle_output = ''
  allocate(damage_isoBrittle_outputID(maxval(phase_Noutput),maxNinstance),      source=undefined_ID)
  allocate(damage_isoBrittle_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(damage_isoBrittle_N(maxNinstance),                                   source=0.0_pReal) 
  allocate(damage_isoBrittle_critStrainEnergy(maxNinstance),                    source=0.0_pReal) 
  allocate(damage_isoBrittle_aTol(maxNinstance),                                source=0.0_pReal) 
 
@@ -158,6 +160,9 @@ subroutine damage_isoBrittle_init(fileUnit)
                                                        IO_lc(IO_stringValue(line,positions,2_pInt))
           end select
 
+       case ('rate_sensitivity_damage')
+         damage_isoBrittle_N(instance) = IO_floatValue(line,positions,2_pInt)
+
        case ('critical_strain_energy')
          damage_isoBrittle_critStrainEnergy(instance) = IO_floatValue(line,positions,2_pInt)
 
@@ -187,7 +192,7 @@ subroutine damage_isoBrittle_init(fileUnit)
        endif
      enddo outputsLoop
 ! Determine size of state array
-     sizeDotState              =   1_pInt
+     sizeDotState              =   2_pInt
      sizeState                 =   2_pInt
                 
      damageState(phase)%sizeState = sizeState
@@ -231,8 +236,7 @@ subroutine damage_isoBrittle_stateInit(phase)
 
  real(pReal), dimension(damageState(phase)%sizeState) :: tempState
 
- tempState(1) = 1.0_pReal
+ tempState = 1.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
@@ -258,11 +262,17 @@ end subroutine damage_isoBrittle_aTolState
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-subroutine damage_isoBrittle_dotState(ipc, ip, el)
+subroutine damage_isoBrittle_dotState(C, Fe, ipc, ip, el)
  use material, only: &
    mappingConstitutive, &
    phase_damageInstance, &
    damageState
+ use math, only : &
+   math_mul33x33, &
+   math_mul66x6, &
+   math_Mandel33to6, &
+   math_transpose33, &
+   math_I3
  use lattice, only: &
    lattice_DamageMobility
 
@@ -271,57 +281,38 @@ subroutine damage_isoBrittle_dotState(ipc, ip, el)
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
    el                                                                                               !< element
- integer(pInt) :: &
-   phase, constituent
-
- phase = mappingConstitutive(2,ipc,ip,el)
- constituent = mappingConstitutive(1,ipc,ip,el)
- 
- damageState(phase)%dotState(1,constituent) = &
-   (1.0_pReal/lattice_DamageMobility(phase))* &
-   (damageState(phase)%state(2,constituent) - &
-    damageState(phase)%state(1,constituent))
-  
-end subroutine damage_isoBrittle_dotState
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief calculates derived quantities from state
-!--------------------------------------------------------------------------------------------------
-subroutine damage_isoBrittle_microstructure(Tstar_v, Fe, ipc, ip, el)
- use material, only: &
-   mappingConstitutive, &
-   phase_damageInstance, &
-   damageState
- use math, only: &
-   math_Mandel6to33, &
-   math_mul33x33, &
-   math_transpose33, &
-   math_I3
-
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
- real(pReal),  intent(in), dimension(6) :: &
-   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor (Mandel)
  real(pReal),  intent(in), dimension(3,3) :: &
    Fe
  integer(pInt) :: &
    phase, constituent, instance
  real(pReal) :: &
-   strain(3,3)
+   drivingForce, &
+   strain(6), &
+   stress(6), &
+   C(6,6)
 
  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)
+ strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
- damageState(phase)%state(2,constituent) = min(damageState(phase)%state0(2,constituent), &
+ stress = math_mul66x6(C,strain) 
-                                               damage_isoBrittle_critStrainEnergy(instance)/ &
-                                               sum(abs(math_Mandel6to33(Tstar_v)*strain)))
  
-end subroutine damage_isoBrittle_microstructure
+ drivingForce = max(0.0_pReal, &
+                    2.0_pReal*damageState(phase)%state(2,constituent)* &
+                    sum(abs(stress*strain))/damage_isoBrittle_critStrainEnergy(instance) - &
+                    damage_isoBrittle_getDamage(ipc, ip, el) + &
+                    damageState(phase)%state(2,constituent))
+ 
+ damageState(phase)%dotState(1,constituent) = &
+   (damageState(phase)%state(2,constituent) - damageState(phase)%state(1,constituent))/ &
+   lattice_DamageMobility(phase)
+
+ damageState(phase)%dotState(2,constituent) = &
+   (exp(-drivingForce**damage_isoBrittle_N(instance)) - 1.0_pReal)/ &
+   lattice_DamageMobility(phase)
+
+end subroutine damage_isoBrittle_dotState
  
 !--------------------------------------------------------------------------------------------------
 !> @brief returns damage
@@ -417,11 +408,39 @@ function damage_isoBrittle_getDamageDiffusion33(ipc, ip, el)
  phase = mappingConstitutive(2,ipc,ip,el)
  constituent = mappingConstitutive(1,ipc,ip,el)
  damage_isoBrittle_getDamageDiffusion33 = &
-   damageState(phase)%state(2,constituent)* &
    lattice_DamageDiffusion33(1:3,1:3,phase)
     
 end function damage_isoBrittle_getDamageDiffusion33
 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns brittle damaged stiffness tensor 
+!--------------------------------------------------------------------------------------------------
+function damage_isoBrittle_getDamagedC66(C, ipc, ip, el)
+ use material, only: &
+   mappingConstitutive, &
+   damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),  intent(in), dimension(6,6) :: &
+   C
+ real(pReal),              dimension(6,6) :: &
+   damage_isoBrittle_getDamagedC66
+ integer(pInt) :: &
+   phase, constituent
+ 
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ damage_isoBrittle_getDamagedC66 = &
+   damageState(phase)%state(2,constituent)* &
+   damageState(phase)%state(2,constituent)* &
+   C
+    
+end function damage_isoBrittle_getDamagedC66
+
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of constitutive results
 !--------------------------------------------------------------------------------------------------

code/damage_isoDuctile.f90

@@ -27,6 +27,7 @@ module damage_isoDuctile
 
  real(pReal),                         dimension(:),           allocatable,         private :: &
    damage_isoDuctile_aTol, &
+   damage_isoDuctile_N, &
    damage_isoDuctile_critpStrain
 
  enum, bind(c) 
@@ -43,7 +44,6 @@ module damage_isoDuctile
    damage_isoDuctile_stateInit, &
    damage_isoDuctile_aTolState, &
    damage_isoDuctile_dotState, &
-   damage_isoDuctile_microstructure, &
    damage_isoDuctile_getDamage, &
    damage_isoDuctile_getSlipDamage, &
    damage_isoDuctile_putLocalDamage, &
@@ -124,6 +124,7 @@ subroutine damage_isoDuctile_init(fileUnit)
           damage_isoDuctile_output = ''
  allocate(damage_isoDuctile_outputID(maxval(phase_Noutput),maxNinstance),      source=undefined_ID)
  allocate(damage_isoDuctile_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(damage_isoDuctile_N(maxNinstance),                                   source=0.0_pReal) 
  allocate(damage_isoDuctile_critpStrain(maxNinstance),                         source=0.0_pReal) 
  allocate(damage_isoDuctile_aTol(maxNinstance),                                source=0.0_pReal) 
 
@@ -158,6 +159,9 @@ subroutine damage_isoDuctile_init(fileUnit)
                                                        IO_lc(IO_stringValue(line,positions,2_pInt))
           end select
 
+       case ('rate_sensitivity_damage')
+         damage_isoDuctile_N(instance) = IO_floatValue(line,positions,2_pInt)
+
        case ('critical_plastic_strain')
          damage_isoDuctile_critpStrain(instance) = IO_floatValue(line,positions,2_pInt)
 
@@ -187,7 +191,7 @@ subroutine damage_isoDuctile_init(fileUnit)
        endif
      enddo outputsLoop
 ! Determine size of state array
-     sizeDotState              =   1_pInt
+     sizeDotState              =   2_pInt
      sizeState                 =   2_pInt
                 
      damageState(phase)%sizeState = sizeState
@@ -231,10 +235,7 @@ subroutine damage_isoDuctile_stateInit(phase)
 
  real(pReal), dimension(damageState(phase)%sizeState) :: tempState
 
- tempState(1) = 1.0_pReal
+ tempState = 1.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
@@ -260,8 +261,9 @@ end subroutine damage_isoDuctile_aTolState
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-subroutine damage_isoDuctile_dotState(ipc, ip, el)
+subroutine damage_isoDuctile_dotState(nSlip,accumulatedSlip,ipc, ip, el)
  use material, only: &
+   phase_damageInstance, &
    mappingConstitutive, &
    damageState
  use math, only: &
@@ -269,39 +271,6 @@ subroutine damage_isoDuctile_dotState(ipc, ip, el)
  use lattice, only: &
    lattice_DamageMobility
 
- implicit none
- integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
- integer(pInt) :: &
-   phase, constituent
-
- phase = mappingConstitutive(2,ipc,ip,el)
- constituent = mappingConstitutive(1,ipc,ip,el)
- 
- damageState(phase)%dotState(1,constituent) = &
-   (1.0_pReal/lattice_DamageMobility(phase))* &
-   (damageState(phase)%state(2,constituent) - &
-    damageState(phase)%state(1,constituent))
-  
-end subroutine damage_isoDuctile_dotState
- 
-!--------------------------------------------------------------------------------------------------
-!> @brief calculates derived quantities from state
-!--------------------------------------------------------------------------------------------------
-subroutine damage_isoDuctile_microstructure(nSlip,accumulatedSlip,ipc, ip, el)
- use material, only: &
-   mappingConstitutive, &
-   phase_damageInstance, &
-   damageState
- use math, only: &
-   math_Mandel6to33, &
-   math_mul33x33, &
-   math_transpose33, &
-   math_I3, &
-   math_norm33
-
  implicit none
  integer(pInt), intent(in) :: &
    nSlip, &
@@ -312,15 +281,28 @@ subroutine damage_isoDuctile_microstructure(nSlip,accumulatedSlip,ipc, ip, el)
    accumulatedSlip
  integer(pInt) :: &
    phase, constituent, instance
+ real(pReal) :: &
+   drivingForce
 
  phase = mappingConstitutive(2,ipc,ip,el)
  constituent = mappingConstitutive(1,ipc,ip,el)
  instance = phase_damageInstance(phase)
- damageState(phase)%state(2,constituent) = min(damageState(phase)%state(2,constituent), &
-                                                     damage_isoDuctile_critpStrain(instance)/ &
-                                                     sum(accumulatedSlip))                         !< akin to damage surface
 
-end subroutine damage_isoDuctile_microstructure
+ drivingForce = max(0.0_pReal, &
+                    2.0_pReal*damageState(phase)%state(2,constituent)* &
+                    sum(accumulatedSlip)/damage_isoDuctile_critpStrain(instance) - &
+                    damage_isoDuctile_getDamage(ipc, ip, el) + &
+                    damageState(phase)%state(2,constituent))
+ 
+ damageState(phase)%dotState(1,constituent) = &
+   (damageState(phase)%state(2,constituent) - damageState(phase)%state(1,constituent))/ &
+   lattice_DamageMobility(phase)
+
+ damageState(phase)%dotState(2,constituent) = &
+   (exp(-drivingForce**damage_isoDuctile_N(instance)) - 1.0_pReal)/ &
+   lattice_DamageMobility(phase)
+  
+end subroutine damage_isoDuctile_dotState
  
 !--------------------------------------------------------------------------------------------------
 !> @brief returns damage