From 5ff3c882b5adb1e1cff7bebfc29e129505817236 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Tue, 11 Nov 2014 22:44:11 +0000 Subject: [PATCH] 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 | 24 ++++---- code/damage_anisoDuctile.f90 | 17 +++--- code/damage_isoBrittle.f90 | 111 ++++++++++++++++++++--------------- code/damage_isoDuctile.f90 | 72 +++++++++-------------- 4 files changed, 112 insertions(+), 112 deletions(-) diff --git a/code/damage_anisoBrittle.f90 b/code/damage_anisoBrittle.f90 index 51d421902..59c4633c1 100644 --- a/code/damage_anisoBrittle.f90 +++ b/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)%state0(index_d,constituent) + damageState(phase)%state(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)* & - max(0.0_pReal,2.0_pReal*damageState(phase)%state(index_d,constituent)* & - damageState(phase)%state(index_o,constituent) - & - nonlocalFactor)**damage_anisoBrittle_N(instance) + (exp(-drivingForce**damage_anisoBrittle_N(instance)) - 1.0_pReal)/ & + lattice_DamageMobility(phase) 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)%state0(index,constituent) + damageState(phase)%state(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)%state0(index,constituent) + damageState(phase)%state(index,constituent)* & + damageState(phase)%state(index,constituent) udotd = & sign(1.0_pReal,traction_d)* & damage_anisoBrittle_sdot_0(instance)* & diff --git a/code/damage_anisoDuctile.f90 b/code/damage_anisoDuctile.f90 index dace9560d..815ec47b7 100644 --- a/code/damage_anisoDuctile.f90 +++ b/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)* & - max(0.0_pReal, & - 2.0_pReal*damageState(phase)%state(index+1,constituent)* & - accumulatedSlip(index)/damage_anisoDuctile_critAccShear(f,instance) - & - nonlocalFactor)**damage_anisoDuctile_N(instance) + (exp(-drivingForce**damage_anisoDuctile_N(instance)) - 1.0_pReal)/ & + lattice_DamageMobility(phase) index = index + 1_pInt enddo enddo diff --git a/code/damage_isoBrittle.f90 b/code/damage_isoBrittle.f90 index c174c80e0..8da7541a9 100644 --- a/code/damage_isoBrittle.f90 +++ b/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(2) = 1.0_pReal + tempState = 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,58 +281,39 @@ 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) - damageState(phase)%state(2,constituent) = min(damageState(phase)%state0(2,constituent), & - damage_isoBrittle_critStrainEnergy(instance)/ & - sum(abs(math_Mandel6to33(Tstar_v)*strain))) + strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3) + stress = math_mul66x6(C,strain) + + 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) -end subroutine damage_isoBrittle_microstructure + 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 !-------------------------------------------------------------------------------------------------- diff --git a/code/damage_isoDuctile.f90 b/code/damage_isoDuctile.f90 index 5b2071837..5b2fe1808 100644 --- a/code/damage_isoDuctile.f90 +++ b/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(2) = 1.0_pReal - - + tempState = 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,16 +281,29 @@ 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 !--------------------------------------------------------------------------------------------------