diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 42c0c624d..5ab7f7ee0 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -551,7 +551,7 @@ end function constitutive_damagedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el) +subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el) use prec, only: & pReal use material, only: & @@ -600,6 +600,8 @@ 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), intent(in) :: & + subdt !< timestep real(pReal), dimension(:), allocatable :: & accumulatedSlip integer(pInt) :: & @@ -620,22 +622,22 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el) select case (phase_damage(material_phase(ipc,ip,el))) case (LOCAL_DAMAGE_isoBrittle_ID) - call damage_isoBrittle_microstructure(constitutive_homogenizedC(ipc,ip,el), Fe, & + call damage_isoBrittle_microstructure(constitutive_homogenizedC(ipc,ip,el), Fe, subdt, & 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) + call damage_isoDuctile_microstructure(nSlip, accumulatedSlip, subdt, ipc, ip, el) case (LOCAL_DAMAGE_anisoBrittle_ID) - call damage_anisoBrittle_microstructure(ipc, ip, el) + call damage_anisoBrittle_microstructure(Tstar_v, subdt, 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) + call damage_anisoDuctile_microstructure(nSlip, accumulatedSlip, subdt, ipc, ip, el) case (LOCAL_DAMAGE_gurson_ID) call damage_gurson_microstructure(ipc, ip, el) case (LOCAL_DAMAGE_phaseField_ID) call damage_phaseField_microstructure(constitutive_homogenizedC(ipc,ip,el), Fe, & constitutive_getVacancyConcentration(ipc, ip, el), & - ipc, ip, el) + subdt, ipc, ip, el) end select @@ -1081,12 +1083,8 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su PLASTICITY_dislokmc_ID, & PLASTICITY_titanmod_ID, & PLASTICITY_nonlocal_ID, & - LOCAL_DAMAGE_isoBrittle_ID, & - LOCAL_DAMAGE_isoDuctile_ID, & - LOCAL_DAMAGE_anisoDuctile_ID, & LOCAL_DAMAGE_anisoBrittle_ID, & LOCAL_DAMAGE_gurson_ID, & - LOCAL_DAMAGE_phaseField_ID, & LOCAL_VACANCY_generation_ID use constitutive_j2, only: & constitutive_j2_dotState @@ -1100,18 +1098,10 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su constitutive_titanmod_dotState use constitutive_nonlocal, only: & constitutive_nonlocal_dotState - use damage_isoBrittle, only: & - damage_isoBrittle_dotState - use damage_isoDuctile, only: & - damage_isoDuctile_dotState use damage_anisoBrittle, only: & damage_anisoBrittle_dotState - use damage_anisoDuctile, only: & - damage_anisoDuctile_dotState use damage_gurson, only: & damage_gurson_dotState - use damage_phaseField, only: & - damage_phaseField_dotState use vacancy_generation, only: & vacancy_generation_dotState @@ -1160,18 +1150,10 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su end select select case (phase_damage(material_phase(ipc,ip,el))) - case (LOCAL_DAMAGE_isoBrittle_ID) - call damage_isoBrittle_dotState(ipc, ip, el) - case (LOCAL_DAMAGE_isoDuctile_ID) - 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 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) - call damage_phaseField_dotState(constitutive_getVacancyConcentration(ipc, ip, el), ipc, ip, el) end select select case (phase_vacancy(material_phase(ipc,ip,el))) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index b8f19e0c8..9dc0b7c36 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -416,7 +416,8 @@ subroutine crystallite_init call constitutive_microstructure( & crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e),g,i,e) ! update dependent state variables to be consistent with basic states + crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_subdt(g,i,e), g,i,e) ! update dependent state variables to be consistent with basic states enddo enddo enddo @@ -1594,7 +1595,8 @@ subroutine crystallite_integrateStateRK4() if (crystallite_todo(g,i,e)) & call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -1911,7 +1913,8 @@ subroutine crystallite_integrateStateRKCK45() if (crystallite_todo(g,i,e)) & call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2159,7 +2162,8 @@ subroutine crystallite_integrateStateRKCK45() if (crystallite_todo(g,i,e)) & call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2392,7 +2396,8 @@ subroutine crystallite_integrateStateAdaptiveEuler() if (crystallite_todo(g,i,e)) & call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO !$OMP END PARALLEL @@ -2728,7 +2733,8 @@ eIter = FEsolving_execElem(1:2) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) & call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO !$OMP END PARALLEL @@ -2975,7 +2981,8 @@ subroutine crystallite_integrateStateFPI() if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) & call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states p = mappingConstitutive(2,g,i,e) c = mappingConstitutive(1,g,i,e) plasticState(p)%previousDotState2(:,c) = plasticState(p)%previousDotState(:,c) diff --git a/code/damage_anisoDuctile.f90 b/code/damage_anisoDuctile.f90 index 054c3f9f9..9af9837c9 100644 --- a/code/damage_anisoDuctile.f90 +++ b/code/damage_anisoDuctile.f90 @@ -50,7 +50,6 @@ module damage_anisoDuctile damage_anisoDuctile_init, & damage_anisoDuctile_stateInit, & damage_anisoDuctile_aTolState, & - damage_anisoDuctile_dotState, & damage_anisoDuctile_microstructure, & damage_anisoDuctile_getDamage, & damage_anisoDuctile_putLocalDamage, & @@ -222,8 +221,9 @@ subroutine damage_anisoDuctile_init(fileUnit) endif enddo outputsLoop ! Determine size of state array - sizeDotState = 1_pInt ! non-local damage + sizeDotState = 0_pInt sizeState = sizeDotState + & + 1_pInt + & ! non-local damage damage_anisoDuctile_totalNslip(instance) damageState(phase)%sizeState = sizeState @@ -293,47 +293,15 @@ end subroutine damage_anisoDuctile_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine damage_anisoDuctile_dotState(ipc, ip, el) - use material, only: & - mappingConstitutive, & - phase_damageInstance, & - 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, & - instance - real(pReal) :: & - localDamage - - phase = mappingConstitutive(2,ipc,ip,el) - constituent = mappingConstitutive(1,ipc,ip,el) - instance = phase_damageInstance(phase) - - localDamage = minval(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) +subroutine damage_anisoDuctile_microstructure(nSlip, accumulatedSlip, subdt, ipc, ip, el) use material, only: & mappingConstitutive, & phase_damageInstance, & damageState use lattice, only: & lattice_maxNslipFamily + use lattice, only: & + lattice_DamageMobility implicit none integer(pInt), intent(in) :: & @@ -343,6 +311,8 @@ subroutine damage_anisoDuctile_microstructure(nSlip, accumulatedSlip, ipc, ip, e el !< element real(pReal), dimension(nSlip), intent(in) :: & accumulatedSlip + real(pReal), intent(in) :: & + subdt integer(pInt) :: & phase, & constituent, & @@ -377,6 +347,13 @@ subroutine damage_anisoDuctile_microstructure(nSlip, accumulatedSlip, ipc, ip, e enddo enddo + localDamage = minval(damageState(phase)%state(2:1+damage_anisoDuctile_totalNslip(instance),constituent)) + + damageState(phase)%state(1,constituent) = & + localDamage + & + (damageState(phase)%subState0(1,constituent) - localDamage)* & + exp(-subdt/lattice_DamageMobility(phase)) + end subroutine damage_anisoDuctile_microstructure !-------------------------------------------------------------------------------------------------- @@ -459,6 +436,8 @@ end function damage_anisoDuctile_getLocalDamage !> @brief returns slip system damage !-------------------------------------------------------------------------------------------------- function damage_anisoDuctile_getSlipDamage(ipc, ip, el) + use numerics, only: & + residualStiffness use material, only: & mappingConstitutive, & phase_damageInstance, & @@ -482,7 +461,8 @@ function damage_anisoDuctile_getSlipDamage(ipc, ip, el) instance = phase_damageInstance(phase) damage_anisoDuctile_getSlipDamage = & - damageState(phase)%state0(2:1+damage_anisoDuctile_totalNslip(instance),constituent) + damageState(phase)%state0(2:1+damage_anisoDuctile_totalNslip(instance),constituent) + & + residualStiffness end function damage_anisoDuctile_getSlipDamage diff --git a/code/damage_isoBrittle.f90 b/code/damage_isoBrittle.f90 index 248f0de05..52687a12e 100644 --- a/code/damage_isoBrittle.f90 +++ b/code/damage_isoBrittle.f90 @@ -42,7 +42,6 @@ module damage_isoBrittle damage_isoBrittle_init, & damage_isoBrittle_stateInit, & damage_isoBrittle_aTolState, & - damage_isoBrittle_dotState, & damage_isoBrittle_microstructure, & damage_isoBrittle_getDamage, & damage_isoBrittle_putLocalDamage, & @@ -196,7 +195,7 @@ subroutine damage_isoBrittle_init(fileUnit) endif enddo outputsLoop ! Determine size of state array - sizeDotState = 1_pInt + sizeDotState = 0_pInt sizeState = 2_pInt damageState(phase)%sizeState = sizeState @@ -266,34 +265,7 @@ end subroutine damage_isoBrittle_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -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) +subroutine damage_isoBrittle_microstructure(C, Fe, subdt, ipc, ip, el) use material, only: & mappingConstitutive, & phase_damageInstance, & @@ -304,6 +276,8 @@ subroutine damage_isoBrittle_microstructure(C, Fe, ipc, ip, el) math_Mandel33to6, & math_transpose33, & math_I3 + use lattice, only: & + lattice_DamageMobility implicit none integer(pInt), intent(in) :: & @@ -312,6 +286,8 @@ subroutine damage_isoBrittle_microstructure(C, Fe, ipc, ip, el) el !< element real(pReal), intent(in), dimension(3,3) :: & Fe + real(pReal), intent(in) :: & + subdt integer(pInt) :: & phase, constituent, instance real(pReal) :: & @@ -329,7 +305,12 @@ subroutine damage_isoBrittle_microstructure(C, Fe, ipc, ip, el) damageState(phase)%state(2,constituent) = & min(damageState(phase)%state0(2,constituent), & damage_isoBrittle_critStrainEnergy(instance)/sum(abs(stress*strain))) - + + damageState(phase)%state(1,constituent) = & + damageState(phase)%state(2,constituent) + & + (damageState(phase)%subState0(1,constituent) - damageState(phase)%state(2,constituent))* & + exp(-subdt/lattice_DamageMobility(phase)) + end subroutine damage_isoBrittle_microstructure !-------------------------------------------------------------------------------------------------- diff --git a/code/damage_isoDuctile.f90 b/code/damage_isoDuctile.f90 index 2c87c0da2..12ea9c7f5 100644 --- a/code/damage_isoDuctile.f90 +++ b/code/damage_isoDuctile.f90 @@ -42,7 +42,6 @@ module damage_isoDuctile damage_isoDuctile_init, & damage_isoDuctile_stateInit, & damage_isoDuctile_aTolState, & - damage_isoDuctile_dotState, & damage_isoDuctile_microstructure, & damage_isoDuctile_getDamage, & damage_isoDuctile_getSlipDamage, & @@ -195,7 +194,7 @@ subroutine damage_isoDuctile_init(fileUnit) endif enddo outputsLoop ! Determine size of state array - sizeDotState = 1_pInt + sizeDotState = 0_pInt sizeState = 2_pInt damageState(phase)%sizeState = sizeState @@ -265,40 +264,15 @@ end subroutine damage_isoDuctile_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -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) +subroutine damage_isoDuctile_microstructure(nSlip,accumulatedSlip,subdt,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) :: & @@ -308,6 +282,8 @@ subroutine damage_isoDuctile_microstructure(nSlip,accumulatedSlip,ipc, ip, el) el !< element real(pReal), dimension(nSlip), intent(in) :: & accumulatedSlip + real(pReal), intent(in) :: & + subdt integer(pInt) :: & phase, constituent, instance @@ -319,6 +295,11 @@ subroutine damage_isoDuctile_microstructure(nSlip,accumulatedSlip,ipc, ip, el) min(damageState(phase)%state0(2,constituent), & damage_isoDuctile_critPlasticStrain(instance)/sum(accumulatedSlip)) + damageState(phase)%state(1,constituent) = & + damageState(phase)%state(2,constituent) + & + (damageState(phase)%subState0(1,constituent) - damageState(phase)%state(2,constituent))* & + exp(-subdt/lattice_DamageMobility(phase)) + end subroutine damage_isoDuctile_microstructure !-------------------------------------------------------------------------------------------------- diff --git a/code/damage_phaseField.f90 b/code/damage_phaseField.f90 index f1c43f352..60bd39d56 100644 --- a/code/damage_phaseField.f90 +++ b/code/damage_phaseField.f90 @@ -43,7 +43,6 @@ module damage_phaseField damage_phaseField_init, & damage_phaseField_stateInit, & damage_phaseField_aTolState, & - damage_phaseField_dotState, & damage_phaseField_microstructure, & damage_phaseField_getDamage, & damage_phaseField_putLocalDamage, & @@ -201,7 +200,7 @@ subroutine damage_phaseField_init(fileUnit) endif enddo outputsLoop ! Determine size of state array - sizeDotState = 1_pInt + sizeDotState = 0_pInt sizeState = 2_pInt damageState(phase)%sizeState = sizeState @@ -271,37 +270,7 @@ end subroutine damage_phaseField_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine damage_phaseField_dotState(Cv, 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 - real(pReal), intent(in) :: & - Cv - 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)*(1.0_pReal - Cv) - & - damageState(phase)%state(1,constituent))/ & - lattice_DamageMobility(phase) - -end subroutine damage_phaseField_dotState - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates derived quantities from state -!-------------------------------------------------------------------------------------------------- -subroutine damage_phaseField_microstructure(C, Fe, Cv, ipc, ip, el) +subroutine damage_phaseField_microstructure(C, Fe, Cv, subdt, ipc, ip, el) use material, only: & mappingConstitutive, & phase_damageInstance, & @@ -312,6 +281,8 @@ 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) :: & @@ -324,6 +295,8 @@ subroutine damage_phaseField_microstructure(C, Fe, Cv, ipc, ip, el) C real(pReal), intent(in) :: & Cv + real(pReal), intent(in) :: & + subdt integer(pInt) :: & phase, constituent, instance real(pReal) :: & @@ -341,6 +314,11 @@ subroutine damage_phaseField_microstructure(C, Fe, Cv, ipc, ip, el) min(damageState(phase)%state0(2,constituent), & damage_phaseField_surfaceEnergy(instance)/ & (2.0_pReal*(sum(abs(stress*strain)) + Cv*damage_phaseField_vacancyFormationEnergy(instance)))) + + damageState(phase)%state(1,constituent) = & + damageState(phase)%state(2,constituent)*(1.0_pReal - Cv) + & + (damageState(phase)%subState0(1,constituent) - damageState(phase)%state(2,constituent)*(1.0_pReal - Cv))* & + exp(-subdt/lattice_DamageMobility(phase)) end subroutine damage_phaseField_microstructure