From 40701cedc476ef9ea385e55dc43ad1ca301fbf27 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Thu, 11 Jun 2015 09:03:51 +0000 Subject: [PATCH] damage driving force is a history variable and damage field is strictly bounded between 0 and 1 --- code/damage_local.f90 | 3 ++- code/numerics.f90 | 9 ++++---- code/source_damage_anisoBrittle.f90 | 5 ++--- code/source_damage_anisoDuctile.f90 | 4 ++-- code/source_damage_isoBrittle.f90 | 34 ++++++++++++++++++----------- code/source_damage_isoDuctile.f90 | 4 ++-- code/source_vacancy_thermalfluc.f90 | 2 +- 7 files changed, 34 insertions(+), 27 deletions(-) diff --git a/code/damage_local.f90 b/code/damage_local.f90 index afc3999a9..7875bdfe7 100644 --- a/code/damage_local.f90 +++ b/code/damage_local.f90 @@ -182,6 +182,7 @@ end subroutine damage_local_init !-------------------------------------------------------------------------------------------------- function damage_local_updateState(subdt, ip, el) use numerics, only: & + residualStiffness, & err_damage_tolAbs, & err_damage_tolRel use material, only: & @@ -206,7 +207,7 @@ function damage_local_updateState(subdt, ip, el) offset = mappingHomogenization(1,ip,el) phi = damageState(homog)%subState0(1,offset) call damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - phi = phi + subdt*phiDot + phi = max(residualStiffness,min(1.0_pReal,phi + subdt*phiDot)) damage_local_updateState = [ abs(phi - damageState(homog)%state(1,offset)) & <= err_damage_tolAbs & diff --git a/code/numerics.f90 b/code/numerics.f90 index 6f8b6542d..183487594 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -137,7 +137,7 @@ module numerics thermalOrder = 2_pInt, & !< order of temperature field shape functions damageOrder = 1_pInt, & !< order of damage field shape functions vacancyfluxOrder = 2_pInt, & !< order of vacancy concentration and chemical potential field shape functions - porosityOrder = 1_pInt, & !< order of porosity field shape functions + porosityOrder = 2_pInt, & !< order of porosity field shape functions hydrogenfluxOrder = 2_pInt !< order of hydrogen concentration and chemical potential field shape functions logical, protected, public :: & BBarStabilisation = .false. @@ -173,12 +173,11 @@ module numerics &-vacancy_pc_type ml & &-vacancy_mg_levels_ksp_type chebyshev & &-vacancy_mg_levels_pc_type sor & - &-porosity_snes_type vinewtonrsls & + &-porosity_snes_type newtonls & &-porosity_snes_atol 1e-8 & - &-porosity_ksp_type preonly & + &-porosity_ksp_type fgmres & &-porosity_ksp_max_it 25 & - &-porosity_pc_type cholesky & - &-porosity_pc_factor_mat_solver_package mumps & + &-porosity_pc_type hypre & &-hydrogen_snes_type newtonls & &-hydrogen_snes_linesearch_type cp & &-hydrogen_snes_atol 1e-9 & diff --git a/code/source_damage_anisoBrittle.f90 b/code/source_damage_anisoBrittle.f90 index 99289411e..2df036a64 100644 --- a/code/source_damage_anisoBrittle.f90 +++ b/code/source_damage_anisoBrittle.f90 @@ -376,10 +376,9 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph sourceOffset = source_damage_anisoBrittle_offset(phase) localphiDot = 1.0_pReal - & - max(1.0_pReal,sourceState(phase)%p(sourceOffset)%state(1,constituent))* & - phi + sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi - dLocalphiDot_dPhi = -max(1.0_pReal,sourceState(phase)%p(sourceOffset)%state(1,constituent)) + dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) end subroutine source_damage_anisobrittle_getRateAndItsTangent diff --git a/code/source_damage_anisoDuctile.f90 b/code/source_damage_anisoDuctile.f90 index 213eb8ebd..9a34c3bf3 100644 --- a/code/source_damage_anisoDuctile.f90 +++ b/code/source_damage_anisoDuctile.f90 @@ -365,10 +365,10 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph sourceOffset = source_damage_anisoDuctile_offset(phase) localphiDot = 1.0_pReal - & - max(1.0_pReal,sourceState(phase)%p(sourceOffset)%state(1,constituent))* & + sourceState(phase)%p(sourceOffset)%state(1,constituent)* & phi - dLocalphiDot_dPhi = -max(1.0_pReal,sourceState(phase)%p(sourceOffset)%state(1,constituent)) + dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) end subroutine source_damage_anisoDuctile_getRateAndItsTangent diff --git a/code/source_damage_isoBrittle.f90 b/code/source_damage_isoBrittle.f90 index 36227b707..1e918aa8f 100644 --- a/code/source_damage_isoBrittle.f90 +++ b/code/source_damage_isoBrittle.f90 @@ -263,31 +263,42 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) el !< element real(pReal), intent(in), dimension(3,3) :: & Fe + real(pReal), intent(in), dimension(6,6) :: & + C integer(pInt) :: & phase, constituent, instance, sourceOffset, mech real(pReal) :: & strain(6), & - C(6,6) + stiffness(6,6), & + strainenergy phase = mappingConstitutive(2,ipc,ip,el) constituent = mappingConstitutive(1,ipc,ip,el) instance = source_damage_isoBrittle_instance(phase) sourceOffset = source_damage_isoBrittle_offset(phase) + stiffness = C do mech = 1_pInt, phase_NstiffnessDegradations(phase) select case(phase_stiffnessDegradation(mech,phase)) case (STIFFNESS_DEGRADATION_porosity_ID) - C = porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* & - porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* & - C + stiffness = porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* & + porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* & + stiffness end select enddo strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3) - sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & - 2.0_pReal*sum(strain*math_mul66x6(C,strain))/source_damage_isoBrittle_critStrainEnergy(instance) - & - sourceState(phase)%p(sourceOffset)%state(1,constituent) - + strainenergy = 2.0_pReal*sum(strain*math_mul66x6(stiffness,strain))/ & + source_damage_isoBrittle_critStrainEnergy(instance) + if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then + sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & + strainenergy - sourceState(phase)%p(sourceOffset)%state(1,constituent) + else + sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & + sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - & + sourceState(phase)%p(sourceOffset)%state(1,constituent) + endif + end subroutine source_damage_isoBrittle_deltaState !-------------------------------------------------------------------------------------------------- @@ -315,11 +326,8 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD constituent = mappingConstitutive(1,ipc,ip,el) sourceOffset = source_damage_isoBrittle_offset(phase) - localphiDot = 1.0_pReal - & - max(1.0_pReal,sourceState(phase)%p(sourceOffset)%state(1,constituent))* & - phi - - dLocalphiDot_dPhi = -max(1.0_pReal,sourceState(phase)%p(sourceOffset)%state(1,constituent)) + localphiDot = 1.0_pReal - phi*sourceState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) end subroutine source_damage_isoBrittle_getRateAndItsTangent diff --git a/code/source_damage_isoDuctile.f90 b/code/source_damage_isoDuctile.f90 index 39d5d72bc..fcf82412d 100644 --- a/code/source_damage_isoDuctile.f90 +++ b/code/source_damage_isoDuctile.f90 @@ -301,10 +301,10 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD sourceOffset = source_damage_isoDuctile_offset(phase) localphiDot = 1.0_pReal - & - max(1.0_pReal,sourceState(phase)%p(sourceOffset)%state(1,constituent))* & + sourceState(phase)%p(sourceOffset)%state(1,constituent)* & phi - dLocalphiDot_dPhi = -max(1.0_pReal,sourceState(phase)%p(sourceOffset)%state(1,constituent)) + dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) end subroutine source_damage_isoDuctile_getRateAndItsTangent diff --git a/code/source_vacancy_thermalfluc.f90 b/code/source_vacancy_thermalfluc.f90 index 1150d1567..41a99b12a 100644 --- a/code/source_vacancy_thermalfluc.f90 +++ b/code/source_vacancy_thermalfluc.f90 @@ -159,7 +159,7 @@ subroutine source_vacancy_thermalfluc_init(fileUnit) sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_thermalfluc_sizePostResults(instance) - allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.0_pReal) + allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.1_pReal) allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)