damage driving force is a history variable and damage field is strictly bounded between 0 and 1

This commit is contained in:
Pratheek Shanthraj 2015-06-11 09:03:51 +00:00
parent e8ee5d6723
commit 40701cedc4
7 changed files with 34 additions and 27 deletions

View File

@ -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 &

View File

@ -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 &

View File

@ -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

View File

@ -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

View File

@ -263,30 +263,41 @@ 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))* &
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))* &
C
stiffness
end select
enddo
strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
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) = &
2.0_pReal*sum(strain*math_mul66x6(C,strain))/source_damage_isoBrittle_critStrainEnergy(instance) - &
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

View File

@ -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

View File

@ -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)