all source mechanisms covered

This commit is contained in:
Sharan Roongta 2020-09-23 01:08:13 +02:00
parent de3e13df88
commit 190b90d3d4
4 changed files with 44 additions and 44 deletions

View File

@ -12,11 +12,11 @@ submodule (constitutive:constitutive_damage) source_damage_anisoBrittle
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
sdot_0, & !< opening rate of cleavage planes dot_o, & !< opening rate of cleavage planes
n !< damage rate sensitivity q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
critDisp, & !< critical displacement s_crit, & !< critical displacement
critLoad !< critical load g_crit !< critical load
real(pReal), dimension(:,:,:,:), allocatable :: & real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems cleavage_systems
integer :: & integer :: &
@ -75,18 +75,18 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray) N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray)
prm%sum_N_cl = sum(abs(N_cl)) prm%sum_N_cl = sum(abs(N_cl))
prm%n = src%get_asFloat('q') prm%q = src%get_asFloat('q')
prm%sdot_0 = src%get_asFloat('dot_o') prm%dot_o = src%get_asFloat('dot_o')
prm%critDisp = src%get_asFloats('s_crit', requiredSize=size(N_cl)) prm%s_crit = src%get_asFloats('s_crit', requiredSize=size(N_cl))
prm%critLoad = src%get_asFloats('g_crit', requiredSize=size(N_cl)) prm%g_crit = src%get_asFloats('g_crit', requiredSize=size(N_cl))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),& prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase%get_asFloat('c/a',defaultVal=0.0_pReal))
! expand: family => system ! expand: family => system
prm%critDisp = math_expand(prm%critDisp,N_cl) prm%s_crit = math_expand(prm%s_crit,N_cl)
prm%critLoad = math_expand(prm%critLoad,N_cl) prm%g_crit = math_expand(prm%g_crit,N_cl)
#if defined (__GFORTRAN__) #if defined (__GFORTRAN__)
prm%output = output_asStrings(src) prm%output = output_asStrings(src)
@ -95,10 +95,10 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
#endif #endif
! sanity checks ! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o' if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o'
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -152,14 +152,14 @@ module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
traction_crit = prm%critLoad(i)*damage(homog)%p(damageOffset)**2.0_pReal traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
= sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
+ prm%sdot_0 / prm%critDisp(i) & + prm%dot_o / prm%s_crit(i) &
* ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%n + & * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%n + & (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%n) (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q)
enddo enddo
end associate end associate

View File

@ -12,9 +12,9 @@ submodule(constitutive:constitutive_damage) source_damage_anisoDuctile
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
n !< damage rate sensitivity q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
critPlasticStrain !< critical plastic strain per slip system gamma_crit !< critical plastic strain per slip system
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
end type tParameters end type tParameters
@ -68,11 +68,11 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray)
prm%n = src%get_asFloat('q') prm%q = src%get_asFloat('q')
prm%critPlasticStrain = src%get_asFloats('gamma_crit',requiredSize=size(N_sl)) prm%gamma_crit = src%get_asFloats('gamma_crit',requiredSize=size(N_sl))
! expand: family => system ! expand: family => system
prm%critPlasticStrain = math_expand(prm%critPlasticStrain,N_sl) prm%gamma_crit = math_expand(prm%gamma_crit,N_sl)
#if defined (__GFORTRAN__) #if defined (__GFORTRAN__)
prm%output = output_asStrings(src) prm%output = output_asStrings(src)
@ -81,8 +81,8 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
#endif #endif
! sanity checks ! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -127,7 +127,7 @@ module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
associate(prm => param(source_damage_anisoDuctile_instance(phase))) associate(prm => param(source_damage_anisoDuctile_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
= sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%n)/prm%critPlasticStrain) = sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit)
end associate end associate
end subroutine source_damage_anisoDuctile_dotState end subroutine source_damage_anisoDuctile_dotState

View File

@ -12,8 +12,8 @@ submodule (constitutive:constitutive_damage) source_damage_isoDuctile
type:: tParameters !< container type for internal constitutive parameters type:: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
critPlasticStrain, & !< critical plastic strain gamma_crit, & !< critical plastic strain
N q
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
end type tParameters end type tParameters
@ -64,8 +64,8 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
associate(prm => param(source_damage_isoDuctile_instance(p))) associate(prm => param(source_damage_isoDuctile_instance(p)))
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
prm%N = src%get_asFloat('q') prm%q = src%get_asFloat('q')
prm%critPlasticStrain = src%get_asFloat('gamma_crit') prm%gamma_crit = src%get_asFloat('gamma_crit')
#if defined (__GFORTRAN__) #if defined (__GFORTRAN__)
prm%output = output_asStrings(src) prm%output = output_asStrings(src)
@ -74,8 +74,8 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
#endif #endif
! sanity checks ! sanity checks
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -120,7 +120,7 @@ module subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
associate(prm => param(source_damage_isoDuctile_instance(phase))) associate(prm => param(source_damage_isoDuctile_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit
end associate end associate
end subroutine source_damage_isoDuctile_dotState end subroutine source_damage_isoDuctile_dotState

View File

@ -13,8 +13,8 @@ submodule(constitutive:constitutive_thermal) source_thermal_externalheat
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
time, & t_n, &
heat_rate f_T
integer :: & integer :: &
nIntervals nIntervals
end type tParameters end type tParameters
@ -64,10 +64,10 @@ module function source_thermal_externalheat_init(source_length) result(mySources
associate(prm => param(source_thermal_externalheat_instance(p))) associate(prm => param(source_thermal_externalheat_instance(p)))
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
prm%time = src%get_asFloats('t_n') prm%t_n = src%get_asFloats('t_n')
prm%nIntervals = size(prm%time) - 1 prm%nIntervals = size(prm%t_n) - 1
prm%heat_rate = src%get_asFloats('f_T',requiredSize = size(prm%time)) prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -121,13 +121,13 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d
associate(prm => param(source_thermal_externalheat_instance(phase))) associate(prm => param(source_thermal_externalheat_instance(phase)))
do interval = 1, prm%nIntervals ! scan through all rate segments do interval = 1, prm%nIntervals ! scan through all rate segments
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%time(interval)) & frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) &
/ (prm%time(interval+1) - prm%time(interval)) ! fractional time within segment / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment
if ( (frac_time < 0.0_pReal .and. interval == 1) & if ( (frac_time < 0.0_pReal .and. interval == 1) &
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
.or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
TDot = prm%heat_rate(interval ) * (1.0_pReal - frac_time) + & TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + &
prm%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries... prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
! ...or extrapolate if outside of bounds ! ...or extrapolate if outside of bounds
enddo enddo
dTDot_dT = 0.0 dTDot_dT = 0.0