From 36dc59b09fb2089da3788fddb001d9e9ebe2e675 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Tue, 25 Nov 2014 23:56:52 +0000 Subject: [PATCH] reverted back to rate independent damage evolution. updated example config files to recent changes --- code/constitutive.f90 | 36 +++++++++--- code/damage_anisoBrittle.f90 | 103 +++++++++++++++++++++++------------ code/damage_anisoDuctile.f90 | 96 ++++++++++++++++++-------------- code/damage_isoBrittle.f90 | 76 ++++++++++++++------------ code/damage_isoDuctile.f90 | 75 +++++++++++++------------ code/damage_phaseField.f90 | 13 +---- 6 files changed, 233 insertions(+), 166 deletions(-) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 9dc3da19f..d548c79aa 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -564,6 +564,8 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el) PLASTICITY_nonlocal_ID, & LOCAL_DAMAGE_isoBrittle_ID, & LOCAL_DAMAGE_isoDuctile_ID, & + LOCAL_DAMAGE_anisoBrittle_ID, & + LOCAL_DAMAGE_anisoDuctile_ID, & LOCAL_DAMAGE_gurson_ID, & LOCAL_DAMAGE_phaseField_ID @@ -575,8 +577,14 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el) constitutive_dislotwin_microstructure use constitutive_dislokmc, only: & constitutive_dislokmc_microstructure - use damage_gurson, only: & - damage_gurson_microstructure + use damage_isoBrittle, only: & + damage_isoBrittle_microstructure + use damage_isoDuctile, only: & + damage_isoDuctile_microstructure + use damage_anisoBrittle, only: & + damage_anisoBrittle_microstructure + use damage_anisoDuctile, only: & + damage_anisoDuctile_microstructure use damage_gurson, only: & damage_gurson_microstructure use damage_phaseField, only: & @@ -592,6 +600,10 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el) real(pReal), intent(in), dimension(3,3) :: & Fe, & !< elastic deformation gradient Fp !< plastic deformation gradient + real(pReal), dimension(:), allocatable :: & + accumulatedSlip + integer(pInt) :: & + nSlip select case (phase_plasticity(material_phase(ipc,ip,el))) @@ -607,6 +619,17 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el) end select select case (phase_damage(material_phase(ipc,ip,el))) + case (LOCAL_DAMAGE_isoBrittle_ID) + call damage_isoBrittle_microstructure(constitutive_homogenizedC(ipc,ip,el), Fe, & + ipc, ip, el) + case (LOCAL_DAMAGE_isoDuctile_ID) + call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,ipc, ip, el) + call damage_isoDuctile_microstructure(nSlip, accumulatedSlip, ipc, ip, el) + case (LOCAL_DAMAGE_anisoBrittle_ID) + call damage_anisoBrittle_microstructure(ipc, ip, el) + case (LOCAL_DAMAGE_anisoDuctile_ID) + call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,ipc, ip, el) + call damage_anisoDuctile_microstructure(nSlip, accumulatedSlip, ipc, ip, el) case (LOCAL_DAMAGE_gurson_ID) call damage_gurson_microstructure(ipc, ip, el) case (LOCAL_DAMAGE_phaseField_ID) @@ -1138,16 +1161,13 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su select case (phase_damage(material_phase(ipc,ip,el))) case (LOCAL_DAMAGE_isoBrittle_ID) - call damage_isoBrittle_dotState(constitutive_homogenizedC(ipc,ip,el), & - FeArray(1:3,1:3,ipc,ip,el), ipc, ip, el) + call damage_isoBrittle_dotState(ipc, ip, el) case (LOCAL_DAMAGE_isoDuctile_ID) - call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,ipc, ip, el) - call damage_isoDuctile_dotState(nSlip,accumulatedSlip,ipc, ip, el) + call damage_isoDuctile_dotState(ipc, ip, el) case (LOCAL_DAMAGE_anisoBrittle_ID) call damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el) case (LOCAL_DAMAGE_anisoDuctile_ID) - call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,ipc,ip,el) - call damage_anisoDuctile_dotState(nSlip, accumulatedSlip, ipc, ip, el) + call damage_anisoDuctile_dotState(ipc, ip, el) case (LOCAL_DAMAGE_gurson_ID) call damage_gurson_dotState(Tstar_v, Lp, ipc, ip, el) case (LOCAL_DAMAGE_phaseField_ID) diff --git a/code/damage_anisoBrittle.f90 b/code/damage_anisoBrittle.f90 index d1efa231f..12ce8018c 100644 --- a/code/damage_anisoBrittle.f90 +++ b/code/damage_anisoBrittle.f90 @@ -55,6 +55,7 @@ module damage_anisoBrittle damage_anisoBrittle_stateInit, & damage_anisoBrittle_aTolState, & damage_anisoBrittle_dotState, & + damage_anisoBrittle_microstructure, & damage_anisoBrittle_LdAndItsTangent, & damage_anisoBrittle_getFd, & damage_anisoBrittle_putFd, & @@ -78,9 +79,6 @@ subroutine damage_anisoBrittle_init(fileUnit) debug_level,& debug_constitutive,& debug_levelBasic - use mesh, only: & - mesh_maxNips, & - mesh_NcpElems use IO, only: & IO_read, & IO_lc, & @@ -95,7 +93,6 @@ subroutine damage_anisoBrittle_init(fileUnit) IO_timeStamp, & IO_EOF use material, only: & - homogenization_maxNgrains, & phase_damage, & phase_damageInstance, & phase_Noutput, & @@ -189,10 +186,10 @@ subroutine damage_anisoBrittle_init(fileUnit) case ('atol_disp') damage_anisoBrittle_aTol_disp(instance) = IO_floatValue(line,positions,2_pInt) - case ('sdot_0') + case ('sdot0') damage_anisoBrittle_sdot_0(instance) = IO_floatValue(line,positions,2_pInt) - case ('rate_sensitivity_damage') + case ('damageratesensitivity') damage_anisoBrittle_N(instance) = IO_floatValue(line,positions,2_pInt) case ('ncleavage') ! @@ -201,12 +198,12 @@ subroutine damage_anisoBrittle_init(fileUnit) damage_anisoBrittle_Ncleavage(j,instance) = IO_intValue(line,positions,1_pInt+j) enddo - case ('critical_displacement') + case ('criticaldisplacement') do j = 1_pInt, Nchunks_CleavageFamilies damage_anisoBrittle_critDisp(j,instance) = IO_floatValue(line,positions,1_pInt+j) enddo - case ('critical_load') + case ('criticalload') do j = 1_pInt, Nchunks_CleavageFamilies damage_anisoBrittle_critLoad(j,instance) = IO_floatValue(line,positions,1_pInt+j) enddo @@ -259,9 +256,9 @@ subroutine damage_anisoBrittle_init(fileUnit) enddo outputsLoop ! Determine size of state array sizeDotState = 1_pInt + & ! non-local damage - damage_anisoBrittle_totalNcleavage(instance) + & ! local damage on each damage system - damage_anisoBrittle_totalNcleavage(instance) ! local damage on each damage system + damage_anisoBrittle_totalNcleavage(instance) ! opening on each damage system sizeState = sizeDotState + & + damage_anisoBrittle_totalNcleavage(instance) + & ! local damage on each damage system 9_pInt ! Fd damageState(phase)%sizeState = sizeState @@ -309,9 +306,9 @@ subroutine damage_anisoBrittle_stateInit(phase,instance) tempState(1) = 1.0_pReal tempState(2 : & - 1 + damage_anisoBrittle_totalNcleavage(instance)) = 1.0_pReal + 1 + damage_anisoBrittle_totalNcleavage(instance)) = 0.0_pReal tempState(2 + damage_anisoBrittle_totalNcleavage(instance): & - 1 +2*damage_anisoBrittle_totalNcleavage(instance)) = 0.0_pReal + 1 +2*damage_anisoBrittle_totalNcleavage(instance)) = 1.0_pReal tempState(2 +2*damage_anisoBrittle_totalNcleavage(instance): & 10+2*damage_anisoBrittle_totalNcleavage(instance)) = reshape(math_I3, shape=[9]) damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:))) @@ -334,9 +331,9 @@ subroutine damage_anisoBrittle_aTolState(phase,instance) tempTol(1) = damage_anisoBrittle_aTol_damage(instance) tempTol(2 : & - 1 + damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_damage(instance) + 1 + damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_disp (instance) tempTol(2 + damage_anisoBrittle_totalNcleavage(instance): & - 1 +2*damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_disp (instance) + 1 +2*damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_damage(instance) tempTol(2 +2*damage_anisoBrittle_totalNcleavage(instance): & 10+2*damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_damage(instance) damageState(phase)%aTolState = tempTol @@ -351,7 +348,6 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el) phase_damageInstance, & damageState use lattice, only: & - lattice_Scleavage, & lattice_Scleavage_v, & lattice_maxNcleavageFamily, & lattice_NcleavageSystem, & @@ -372,7 +368,7 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el) real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, udott, udotn, & - drivingForce, nonlocalFactor, localDamage + nonlocalFactor, localDamage phase = mappingConstitutive(2,ipc,ip,el) constituent = mappingConstitutive(1,ipc,ip,el) @@ -381,10 +377,12 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el) localDamage = max(0.0_pReal, & 1.0_pReal - sum(1.0_pReal - damageState(phase)% & state(2:1+damage_anisoBrittle_totalNcleavage(instance),constituent))) + damageState(phase)%dotState(1,constituent) = & + (localDamage - damageState(phase)%state(1,constituent))/lattice_DamageMobility(phase) nonlocalFactor = damage_anisoBrittle_getDamage(ipc, ip, el) - localDamage - index_d = 2_pInt - index_o = 2_pInt + damage_anisoBrittle_totalNcleavage(instance) + index_o = 2_pInt + index_d = 2_pInt + damage_anisoBrittle_totalNcleavage(instance) do f = 1_pInt,lattice_maxNcleavageFamily index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family @@ -408,21 +406,61 @@ 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) = & - (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 + enddo + +end subroutine damage_anisoBrittle_dotState + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine damage_anisoBrittle_microstructure(ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + use lattice, only: & + lattice_maxNcleavageFamily, & + lattice_NcleavageSystem + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + integer(pInt) :: & + phase, & + constituent, & + instance, & + f, i, index_d, index_o, index_myFamily + real(pReal) :: & + nonlocalFactor, localDamage + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + localDamage = max(0.0_pReal, & + 1.0_pReal - sum(1.0_pReal - damageState(phase)% & + state(2:1+damage_anisoBrittle_totalNcleavage(instance),constituent))) + nonlocalFactor = damage_anisoBrittle_getDamage(ipc, ip, el) - localDamage + + index_o = 2_pInt + index_d = 2_pInt + damage_anisoBrittle_totalNcleavage(instance) + do f = 1_pInt,lattice_maxNcleavageFamily + index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family + do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family + damageState(phase)%state(index_d,constituent) = & + min(damageState(phase)%state0(index_d,constituent), & + 1.0_pReal/max(0.0_pReal,damageState(phase)%state(index_o,constituent) - & + nonlocalFactor)) index_d = index_d + 1_pInt; index_o = index_o + 1_pInt enddo enddo - damageState(phase)%dotState(1,constituent) = & - (localDamage - damageState(phase)%state(1,constituent))/lattice_DamageMobility(phase) -end subroutine damage_anisoBrittle_dotState +end subroutine damage_anisoBrittle_microstructure !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient @@ -468,7 +506,7 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, Ld = 0.0_pReal dLd_dTstar3333 = 0.0_pReal - index = 2_pInt + index = 2_pInt + damage_anisoBrittle_totalNcleavage(instance) do f = 1_pInt,lattice_maxNcleavageFamily index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family @@ -599,7 +637,7 @@ subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el) instance = phase_damageInstance(phase) Ld = 0.0_pReal - index = 2_pInt + index = 2_pInt + damage_anisoBrittle_totalNcleavage(instance) do f = 1_pInt,lattice_maxNcleavageFamily index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family @@ -739,7 +777,6 @@ end function damage_anisoBrittle_getDamage subroutine damage_anisoBrittle_putLocalDamage(ipc, ip, el, localDamage) use material, only: & mappingConstitutive, & - phase_damageInstance, & damageState implicit none @@ -761,7 +798,6 @@ end subroutine damage_anisoBrittle_putLocalDamage function damage_anisoBrittle_getLocalDamage(ipc, ip, el) use material, only: & mappingConstitutive, & - phase_damageInstance, & damageState implicit none @@ -783,8 +819,7 @@ end function damage_anisoBrittle_getLocalDamage function damage_anisoBrittle_postResults(ipc,ip,el) use material, only: & mappingConstitutive, & - phase_damageInstance,& - damageState + phase_damageInstance implicit none integer(pInt), intent(in) :: & diff --git a/code/damage_anisoDuctile.f90 b/code/damage_anisoDuctile.f90 index f518afa16..85acf526d 100644 --- a/code/damage_anisoDuctile.f90 +++ b/code/damage_anisoDuctile.f90 @@ -32,11 +32,10 @@ module damage_anisoDuctile damage_anisoDuctile_Nslip !< number of slip systems per family real(pReal), dimension(:), allocatable, private :: & - damage_anisoDuctile_aTol_damage, & - damage_anisoDuctile_N + damage_anisoDuctile_aTol_damage real(pReal), dimension(:,:), allocatable, private :: & - damage_anisoDuctile_critAccShear + damage_anisoDuctile_critPlasticStrain enum, bind(c) enumerator :: undefined_ID, & @@ -52,6 +51,7 @@ module damage_anisoDuctile damage_anisoDuctile_stateInit, & damage_anisoDuctile_aTolState, & damage_anisoDuctile_dotState, & + damage_anisoDuctile_microstructure, & damage_anisoDuctile_getDamage, & damage_anisoDuctile_putLocalDamage, & damage_anisoDuctile_getLocalDamage, & @@ -71,9 +71,6 @@ subroutine damage_anisoDuctile_init(fileUnit) debug_level,& debug_constitutive,& debug_levelBasic - use mesh, only: & - mesh_maxNips, & - mesh_NcpElems use IO, only: & IO_read, & IO_lc, & @@ -88,7 +85,6 @@ subroutine damage_anisoDuctile_init(fileUnit) IO_timeStamp, & IO_EOF use material, only: & - homogenization_maxNgrains, & phase_damage, & phase_damageInstance, & phase_Noutput, & @@ -136,11 +132,10 @@ subroutine damage_anisoDuctile_init(fileUnit) damage_anisoDuctile_output = '' allocate(damage_anisoDuctile_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) allocate(damage_anisoDuctile_Noutput(maxNinstance), source=0_pInt) - allocate(damage_anisoDuctile_critAccShear(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(damage_anisoDuctile_critPlasticStrain(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) 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_N(maxNinstance), source=0.0_pReal) rewind(fileUnit) phase = 0_pInt @@ -176,18 +171,15 @@ subroutine damage_anisoDuctile_init(fileUnit) case ('atol_damage') damage_anisoDuctile_aTol_damage(instance) = IO_floatValue(line,positions,2_pInt) - case ('rate_sensitivity_damage') - damage_anisoDuctile_N(instance) = IO_floatValue(line,positions,2_pInt) - case ('nslip') ! Nchunks_SlipFamilies = positions(1) - 1_pInt do j = 1_pInt, Nchunks_SlipFamilies damage_anisoDuctile_Nslip(j,instance) = IO_intValue(line,positions,1_pInt+j) enddo - case ('critical_accshear') + case ('criticalplasticstrain') do j = 1_pInt, Nchunks_SlipFamilies - damage_anisoDuctile_critAccShear(j,instance) = IO_floatValue(line,positions,1_pInt+j) + damage_anisoDuctile_critPlasticStrain(j,instance) = IO_floatValue(line,positions,1_pInt+j) enddo end select @@ -205,10 +197,8 @@ subroutine damage_anisoDuctile_init(fileUnit) damage_anisoDuctile_totalNslip(instance) = sum(damage_anisoDuctile_Nslip(:,instance)) if (damage_anisoDuctile_aTol_damage(instance) < 0.0_pReal) & damage_anisoDuctile_aTol_damage(instance) = 1.0e-3_pReal ! default absolute tolerance 1e-3 - if (any(damage_anisoDuctile_critAccShear(:,instance) < 0.0_pReal)) & - call IO_error(211_pInt,el=instance,ext_msg='critical_accshear ('//LOCAL_DAMAGE_anisoDuctile_LABEL//')') - if (damage_anisoDuctile_N(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity_damage ('//LOCAL_DAMAGE_anisoDuctile_LABEL//')') + if (any(damage_anisoDuctile_critPlasticStrain(:,instance) < 0.0_pReal)) & + call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//LOCAL_DAMAGE_anisoDuctile_LABEL//')') endif myPhase enddo sanityChecks @@ -232,9 +222,9 @@ subroutine damage_anisoDuctile_init(fileUnit) endif enddo outputsLoop ! Determine size of state array - sizeDotState = 1_pInt + & - damage_anisoDuctile_totalNslip(instance) ! non-local damage - sizeState = sizeDotState + sizeDotState = 1_pInt ! non-local damage + sizeState = sizeDotState + & + damage_anisoDuctile_totalNslip(instance) damageState(phase)%sizeState = sizeState damageState(phase)%sizeDotState = sizeDotState @@ -271,8 +261,6 @@ end subroutine damage_anisoDuctile_init subroutine damage_anisoDuctile_stateInit(phase) use material, only: & damageState - use math, only: & - math_I3 implicit none integer(pInt), intent(in) :: phase !< number specifying the phase of the damage @@ -305,15 +293,49 @@ end subroutine damage_anisoDuctile_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine damage_anisoDuctile_dotState(nSlip, accumulatedSlip, ipc, ip, el) +subroutine damage_anisoDuctile_dotState(ipc, ip, el) use material, only: & mappingConstitutive, & phase_damageInstance, & damageState use lattice, only: & - lattice_maxNslipFamily, & lattice_DamageMobility + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + integer(pInt) :: & + phase, & + constituent, & + instance + real(pReal) :: & + localDamage + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + localDamage = max(0.0_pReal, & + 1.0_pReal - sum(1.0_pReal - damageState(phase)% & + state(2:1+damage_anisoDuctile_totalNslip(instance),constituent))) + damageState(phase)%dotState(1,constituent) = & + (localDamage - damageState(phase)%state(1,constituent))/lattice_DamageMobility(phase) + +end subroutine damage_anisoDuctile_dotState + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine damage_anisoDuctile_microstructure(nSlip, accumulatedSlip, ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + use lattice, only: & + lattice_maxNslipFamily + implicit none integer(pInt), intent(in) :: & nSlip, & @@ -329,7 +351,6 @@ subroutine damage_anisoDuctile_dotState(nSlip, accumulatedSlip, ipc, ip, el) index, f, i real(pReal) :: & localDamage, & - drivingForce, & nonlocalFactor phase = mappingConstitutive(2,ipc,ip,el) @@ -344,20 +365,16 @@ 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) = & - (exp(-drivingForce**damage_anisoDuctile_N(instance)) - 1.0_pReal)/ & - lattice_DamageMobility(phase) + damageState(phase)%state(index+1,constituent) = & + min(damageState(phase)%state0(index+1,constituent), & + 1.0_pReal/max(0.0_pReal,accumulatedSlip(index)/ & + damage_anisoDuctile_critPlasticStrain(f,instance) - & + nonlocalFactor)) index = index + 1_pInt enddo enddo - damageState(phase)%dotState(1,constituent) = & - (localDamage - damageState(phase)%state(1,constituent))/lattice_DamageMobility(phase) -end subroutine damage_anisoDuctile_dotState +end subroutine damage_anisoDuctile_microstructure !-------------------------------------------------------------------------------------------------- !> @brief returns damage @@ -396,7 +413,6 @@ end function damage_anisoDuctile_getDamage subroutine damage_anisoDuctile_putLocalDamage(ipc, ip, el, localDamage) use material, only: & mappingConstitutive, & - phase_damageInstance, & damageState implicit none @@ -418,7 +434,6 @@ end subroutine damage_anisoDuctile_putLocalDamage function damage_anisoDuctile_getLocalDamage(ipc, ip, el) use material, only: & mappingConstitutive, & - phase_damageInstance, & damageState implicit none @@ -442,8 +457,6 @@ function damage_anisoDuctile_getSlipDamage(ipc, ip, el) mappingConstitutive, & phase_damageInstance, & damageState - use lattice, only: & - lattice_maxNslipFamily implicit none integer(pInt), intent(in) :: & @@ -474,8 +487,7 @@ end function damage_anisoDuctile_getSlipDamage function damage_anisoDuctile_postResults(ipc,ip,el) use material, only: & mappingConstitutive, & - phase_damageInstance,& - damageState + phase_damageInstance implicit none integer(pInt), intent(in) :: & diff --git a/code/damage_isoBrittle.f90 b/code/damage_isoBrittle.f90 index fa62f4503..67c4be0d0 100644 --- a/code/damage_isoBrittle.f90 +++ b/code/damage_isoBrittle.f90 @@ -27,7 +27,6 @@ module damage_isoBrittle real(pReal), dimension(:), allocatable, private :: & damage_isoBrittle_aTol, & - damage_isoBrittle_N, & damage_isoBrittle_critStrainEnergy enum, bind(c) @@ -44,6 +43,7 @@ module damage_isoBrittle damage_isoBrittle_stateInit, & damage_isoBrittle_aTolState, & damage_isoBrittle_dotState, & + damage_isoBrittle_microstructure, & damage_isoBrittle_getDamage, & damage_isoBrittle_putLocalDamage, & damage_isoBrittle_getLocalDamage, & @@ -64,9 +64,6 @@ subroutine damage_isoBrittle_init(fileUnit) debug_level,& debug_constitutive,& debug_levelBasic - use mesh, only: & - mesh_maxNips, & - mesh_NcpElems use IO, only: & IO_read, & IO_lc, & @@ -81,7 +78,6 @@ subroutine damage_isoBrittle_init(fileUnit) IO_timeStamp, & IO_EOF use material, only: & - homogenization_maxNgrains, & phase_damage, & phase_damageInstance, & phase_Noutput, & @@ -125,7 +121,6 @@ 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) @@ -160,10 +155,7 @@ 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') + case ('criticalstrainenergy') damage_isoBrittle_critStrainEnergy(instance) = IO_floatValue(line,positions,2_pInt) case ('atol_damage') @@ -182,9 +174,7 @@ subroutine damage_isoBrittle_init(fileUnit) if (damage_isoBrittle_aTol(instance) < 0.0_pReal) & damage_isoBrittle_aTol(instance) = 1.0e-3_pReal ! default absolute tolerance 1e-3 if (damage_isoBrittle_critStrainEnergy(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='critical_strain_energy ('//LOCAL_DAMAGE_isoBrittle_LABEL//')') - if (damage_isoBrittle_N(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity_damage ('//LOCAL_DAMAGE_isoBrittle_LABEL//')') + call IO_error(211_pInt,el=instance,ext_msg='criticalStrainEnergy ('//LOCAL_DAMAGE_isoBrittle_LABEL//')') endif myPhase enddo sanityChecks @@ -206,7 +196,7 @@ subroutine damage_isoBrittle_init(fileUnit) endif enddo outputsLoop ! Determine size of state array - sizeDotState = 2_pInt + sizeDotState = 1_pInt sizeState = 2_pInt damageState(phase)%sizeState = sizeState @@ -276,7 +266,34 @@ end subroutine damage_isoBrittle_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine damage_isoBrittle_dotState(C, Fe, ipc, ip, el) +subroutine damage_isoBrittle_dotState(ipc, ip, el) + use material, only: & + mappingConstitutive, & + damageState + 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) = & + (damageState(phase)%state(2,constituent) - damageState(phase)%state(1,constituent))/ & + lattice_DamageMobility(phase) + +end subroutine damage_isoBrittle_dotState + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine damage_isoBrittle_microstructure(C, Fe, ipc, ip, el) use material, only: & mappingConstitutive, & phase_damageInstance, & @@ -287,8 +304,6 @@ subroutine damage_isoBrittle_dotState(C, Fe, ipc, ip, el) math_Mandel33to6, & math_transpose33, & math_I3 - use lattice, only: & - lattice_DamageMobility implicit none integer(pInt), intent(in) :: & @@ -300,7 +315,6 @@ subroutine damage_isoBrittle_dotState(C, Fe, ipc, ip, el) integer(pInt) :: & phase, constituent, instance real(pReal) :: & - drivingForce, & strain(6), & stress(6), & C(6,6) @@ -312,21 +326,11 @@ subroutine damage_isoBrittle_dotState(C, Fe, ipc, ip, el) 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) + damageState(phase)%state(2,constituent) = & + min(damageState(phase)%state0(2,constituent), & + damage_isoBrittle_critStrainEnergy(instance)/sum(abs(stress*strain))) - damageState(phase)%dotState(2,constituent) = & - (exp(-drivingForce**damage_isoBrittle_N(instance)) - 1.0_pReal)/ & - lattice_DamageMobility(phase) - -end subroutine damage_isoBrittle_dotState +end subroutine damage_isoBrittle_microstructure !-------------------------------------------------------------------------------------------------- !> @brief returns damage @@ -449,8 +453,8 @@ function damage_isoBrittle_getDamagedC66(C, ipc, ip, el) 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)* & + damageState(phase)%state(1,constituent)* & + damageState(phase)%state(1,constituent)* & C end function damage_isoBrittle_getDamagedC66 @@ -461,8 +465,8 @@ end function damage_isoBrittle_getDamagedC66 function damage_isoBrittle_postResults(ipc,ip,el) use material, only: & mappingConstitutive, & - phase_damageInstance,& - damageState + damageState, & + phase_damageInstance implicit none integer(pInt), intent(in) :: & diff --git a/code/damage_isoDuctile.f90 b/code/damage_isoDuctile.f90 index 00130d0ba..25c04c299 100644 --- a/code/damage_isoDuctile.f90 +++ b/code/damage_isoDuctile.f90 @@ -27,8 +27,7 @@ module damage_isoDuctile real(pReal), dimension(:), allocatable, private :: & damage_isoDuctile_aTol, & - damage_isoDuctile_N, & - damage_isoDuctile_critpStrain + damage_isoDuctile_critPlasticStrain enum, bind(c) enumerator :: undefined_ID, & @@ -44,6 +43,7 @@ module damage_isoDuctile damage_isoDuctile_stateInit, & damage_isoDuctile_aTolState, & damage_isoDuctile_dotState, & + damage_isoDuctile_microstructure, & damage_isoDuctile_getDamage, & damage_isoDuctile_getSlipDamage, & damage_isoDuctile_putLocalDamage, & @@ -63,9 +63,6 @@ subroutine damage_isoDuctile_init(fileUnit) debug_level,& debug_constitutive,& debug_levelBasic - use mesh, only: & - mesh_maxNips, & - mesh_NcpElems use IO, only: & IO_read, & IO_lc, & @@ -80,7 +77,6 @@ subroutine damage_isoDuctile_init(fileUnit) IO_timeStamp, & IO_EOF use material, only: & - homogenization_maxNgrains, & phase_damage, & phase_damageInstance, & phase_Noutput, & @@ -124,8 +120,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_critPlasticStrain(maxNinstance), source=0.0_pReal) allocate(damage_isoDuctile_aTol(maxNinstance), source=0.0_pReal) rewind(fileUnit) @@ -159,11 +154,8 @@ 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) + case ('criticalplasticstrain') + damage_isoDuctile_critPlasticStrain(instance) = IO_floatValue(line,positions,2_pInt) case ('atol_damage') damage_isoDuctile_aTol(instance) = IO_floatValue(line,positions,2_pInt) @@ -179,10 +171,8 @@ subroutine damage_isoDuctile_init(fileUnit) ! sanity checks if (damage_isoDuctile_aTol(instance) < 0.0_pReal) & damage_isoDuctile_aTol(instance) = 1.0e-3_pReal ! default absolute tolerance 1e-3 - if (damage_isoDuctile_critpStrain(instance) <= 0.0_pReal) & + if (damage_isoDuctile_critPlasticStrain(instance) <= 0.0_pReal) & call IO_error(211_pInt,el=instance,ext_msg='critical_plastic_strain ('//LOCAL_DAMAGE_isoDuctile_LABEL//')') - if (damage_isoDuctile_N(instance) <= 0.0_pReal) & - call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity_damage ('//LOCAL_DAMAGE_isoDuctile_LABEL//')') endif myPhase enddo sanityChecks @@ -205,7 +195,7 @@ subroutine damage_isoDuctile_init(fileUnit) endif enddo outputsLoop ! Determine size of state array - sizeDotState = 2_pInt + sizeDotState = 1_pInt sizeState = 2_pInt damageState(phase)%sizeState = sizeState @@ -275,15 +265,40 @@ end subroutine damage_isoDuctile_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine damage_isoDuctile_dotState(nSlip,accumulatedSlip,ipc, ip, el) +subroutine damage_isoDuctile_dotState(ipc, ip, el) + use material, only: & + mappingConstitutive, & + damageState + 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) = & + (damageState(phase)%state(2,constituent) - damageState(phase)%state(1,constituent))/ & + lattice_DamageMobility(phase) + +end subroutine damage_isoDuctile_dotState + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine damage_isoDuctile_microstructure(nSlip,accumulatedSlip,ipc, ip, el) use material, only: & phase_damageInstance, & mappingConstitutive, & damageState use math, only: & math_norm33 - use lattice, only: & - lattice_DamageMobility implicit none integer(pInt), intent(in) :: & @@ -295,28 +310,16 @@ subroutine damage_isoDuctile_dotState(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) - 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)%state(2,constituent) = & + min(damageState(phase)%state0(2,constituent), & + damage_isoDuctile_critPlasticStrain(instance)/sum(accumulatedSlip)) - 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 +end subroutine damage_isoDuctile_microstructure !-------------------------------------------------------------------------------------------------- !> @brief returns damage diff --git a/code/damage_phaseField.f90 b/code/damage_phaseField.f90 index 298acfdae..ab9ae0078 100644 --- a/code/damage_phaseField.f90 +++ b/code/damage_phaseField.f90 @@ -1,5 +1,5 @@ !-------------------------------------------------------------------------------------------------- -! $Id: damage_phaseField.f90 3736 2014-11-21 13:12:54Z MPIE\l.sharma $ +! $Id$ !-------------------------------------------------------------------------------------------------- !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH @@ -65,9 +65,6 @@ subroutine damage_phaseField_init(fileUnit) debug_level,& debug_constitutive,& debug_levelBasic - use mesh, only: & - mesh_maxNips, & - mesh_NcpElems use IO, only: & IO_read, & IO_lc, & @@ -82,7 +79,6 @@ subroutine damage_phaseField_init(fileUnit) IO_timeStamp, & IO_EOF use material, only: & - homogenization_maxNgrains, & phase_damage, & phase_damageInstance, & phase_Noutput, & @@ -109,7 +105,7 @@ subroutine damage_phaseField_init(fileUnit) mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- damage_'//LOCAL_damage_phaseField_label//' init -+>>>' - write(6,'(a)') ' $Id: damage_phaseField.f90 3736 2014-11-21 13:12:54Z MPIE\l.sharma $' + write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" endif mainProcess @@ -316,8 +312,6 @@ subroutine damage_phaseField_microstructure(C, Fe, Cv, ipc, ip, el) math_Mandel33to6, & math_transpose33, & math_I3 - use lattice, only: & - lattice_DamageMobility implicit none integer(pInt), intent(in) :: & @@ -454,8 +448,7 @@ end function damage_phaseField_getDamageDiffusion33 !-------------------------------------------------------------------------------------------------- function damage_phaseField_getDamagedC66(C, ipc, ip, el) use material, only: & - mappingConstitutive, & - damageState + mappingConstitutive implicit none integer(pInt), intent(in) :: &