223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

updated damage models: 

- coupling to plasticity handled within damage module instead of plasticity module
- anisotropic models more stable
This commit is contained in:
Pratheek Shanthraj 2015-01-29 13:56:09 +00:00
parent 2f68bb0681
commit 0b59519a2a
11 changed files with 165 additions and 417 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

130
code/constitutive.f90
View File

@ -38,7 +38,6 @@ module constitutive
constitutive_getLocalDamage, &
constitutive_putLocalDamage, &
constitutive_getDamage, &
constitutive_getSlipDamage, &
constitutive_getDamageDiffusion33, &
constitutive_getAdiabaticTemperature, &
constitutive_putAdiabaticTemperature, &
@ -701,23 +700,17 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el)
use plastic_j2, only: &
plastic_j2_LpAndItsTangent
use plastic_phenopowerlaw, only: &
plastic_phenopowerlaw_LpAndItsTangent, &
plastic_phenopowerlaw_totalNslip
plastic_phenopowerlaw_LpAndItsTangent
use plastic_dislotwin, only: &
plastic_dislotwin_LpAndItsTangent, &
plastic_dislotwin_totalNslip
plastic_dislotwin_LpAndItsTangent
use plastic_dislokmc, only: &
plastic_dislokmc_LpAndItsTangent, &
plastic_dislokmc_totalNslip
plastic_dislokmc_LpAndItsTangent
use plastic_disloucla, only: &
plastic_disloucla_LpAndItsTangent, &
plastic_disloucla_totalNslip
plastic_disloucla_LpAndItsTangent
use plastic_titanmod, only: &
plastic_titanmod_LpAndItsTangent, &
plastic_titanmod_totalNslip
plastic_titanmod_LpAndItsTangent
use plastic_nonlocal, only: &
plastic_nonlocal_LpAndItsTangent, &
totalNslip
plastic_nonlocal_LpAndItsTangent
implicit none
integer(pInt), intent(in) :: &
@ -730,8 +723,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el)
Lp !< plastic velocity gradient
real(pReal), intent(out), dimension(9,9) :: &
dLp_dTstar !< derivative of Lp with respect to Tstar (4th-order tensor)
integer(pInt) :: &
nSlip
select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -739,44 +730,28 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el)
Lp = 0.0_pReal
dLp_dTstar = 0.0_pReal
case (PLASTICITY_J2_ID)
nSlip = 1_pInt
call plastic_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), &
ipc,ip,el)
call plastic_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID)
nSlip = plastic_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))
call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), &
ipc,ip,el)
call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID)
nSlip = totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))
call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
constitutive_getTemperature(ipc,ip,el), &
nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), &
ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID)
nSlip = plastic_dislotwin_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))
call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
constitutive_getTemperature(ipc,ip,el), &
nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), &
ipc,ip,el)
case (PLASTICITY_DISLOKMC_ID)
nSlip = plastic_dislokmc_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))
call plastic_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
constitutive_getTemperature(ipc,ip,el), &
nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), &
ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID)
nSlip = plastic_disloucla_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))
call plastic_disloucla_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
constitutive_getTemperature(ipc,ip,el), &
nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), &
ipc,ip,el)
case (PLASTICITY_TITANMOD_ID)
nSlip = plastic_titanmod_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))
call plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
constitutive_getTemperature(ipc,ip,el), &
nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), &
ipc,ip,el)
end select
@ -1163,23 +1138,17 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
use plastic_j2, only: &
plastic_j2_dotState
use plastic_phenopowerlaw, only: &
plastic_phenopowerlaw_dotState, &
plastic_phenopowerlaw_totalNslip
plastic_phenopowerlaw_dotState
use plastic_dislotwin, only: &
plastic_dislotwin_dotState, &
plastic_dislotwin_totalNslip
plastic_dislotwin_dotState
use plastic_dislokmc, only: &
plastic_dislokmc_dotState, &
plastic_dislokmc_totalNslip
plastic_dislokmc_dotState
use plastic_disloucla, only: &
plastic_disloucla_dotState, &
plastic_disloucla_totalNslip
plastic_disloucla_dotState
use plastic_titanmod, only: &
plastic_titanmod_dotState, &
plastic_titanmod_totalNslip
plastic_titanmod_dotState
use plastic_nonlocal, only: &
plastic_nonlocal_dotState, &
totalNslip
plastic_nonlocal_dotState
use damage_gurson, only: &
damage_gurson_dotState
@ -1203,38 +1172,30 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
tick, tock, &
tickrate, &
maxticks
integer(pInt) :: &
nSlip
if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) &
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_J2_ID)
nSlip = 1_pInt
call plastic_j2_dotState (Tstar_v,nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el)
call plastic_j2_dotState (Tstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID)
nSlip = plastic_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))
call plastic_phenopowerlaw_dotState(Tstar_v,nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el)
call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID)
nSlip = plastic_dislotwin_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))
call plastic_dislotwin_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), &
nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el)
ipc,ip,el)
case (PLASTICITY_DISLOKMC_ID)
nSlip = plastic_dislokmc_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))
call plastic_dislokmc_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), &
nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el)
ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID)
nSlip = plastic_disloucla_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))
call plastic_disloucla_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), &
nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el)
call plastic_disloucla_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), &
ipc,ip,el)
case (PLASTICITY_TITANMOD_ID)
call plastic_titanmod_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), &
ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID)
nSlip = totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))
call plastic_nonlocal_dotState (Tstar_v,FeArray,FpArray,constitutive_getTemperature(ipc,ip,el), &
nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),subdt,subfracArray,ip,el)
subdt,subfracArray,ip,el)
end select
select case (phase_damage(material_phase(ipc,ip,el)))
@ -1497,52 +1458,6 @@ function constitutive_getDamage(ipc, ip, el)
end function constitutive_getDamage
!--------------------------------------------------------------------------------------------------
!> @brief Returns the damage on each slip system
!--------------------------------------------------------------------------------------------------
function constitutive_getSlipDamage(nSlip, Tstar_v, ipc, ip, el)
use prec, only: &
pReal
use material, only: &
material_phase, &
LOCAL_DAMAGE_isoDuctile_ID, &
LOCAL_DAMAGE_anisoDuctile_ID, &
LOCAL_DAMAGE_gurson_ID, &
phase_damage
use damage_isoDuctile, only: &
damage_isoDuctile_getSlipDamage
use damage_anisoDuctile, only: &
damage_anisoDuctile_getSlipDamage
use damage_gurson, only: &
damage_gurson_getSlipDamage
implicit none
integer(pInt), intent(in) :: &
nSlip, &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress
real(pReal), dimension(nSlip) :: &
constitutive_getSlipDamage
select case (phase_damage(material_phase(ipc,ip,el)))
case (LOCAL_DAMAGE_isoDuctile_ID)
constitutive_getSlipDamage = damage_isoDuctile_getSlipDamage(ipc, ip, el)
case (LOCAL_DAMAGE_anisoDuctile_ID)
constitutive_getSlipDamage = damage_anisoDuctile_getSlipDamage(ipc, ip, el)
case (LOCAL_DAMAGE_gurson_ID)
constitutive_getSlipDamage = damage_gurson_getSlipDamage(Tstar_v, ipc, ip, el)
case default
constitutive_getSlipDamage = 1.0_pReal
end select
end function constitutive_getSlipDamage
!--------------------------------------------------------------------------------------------------
!> @brief returns damage diffusion tensor
!--------------------------------------------------------------------------------------------------
@ -1569,12 +1484,13 @@ function constitutive_getDamageDiffusion33(ipc, ip, el)
real(pReal), dimension(3,3) :: &
constitutive_getDamageDiffusion33
constitutive_getDamageDiffusion33 = lattice_DamageDiffusion33(1:3,1:3,material_phase(ipc,ip,el))
select case(phase_damage(material_phase(ipc,ip,el)))
case (LOCAL_DAMAGE_isoBrittle_ID)
constitutive_getDamageDiffusion33 = damage_isoBrittle_getDamageDiffusion33(ipc, ip, el)
case (LOCAL_DAMAGE_phaseField_ID)
constitutive_getDamageDiffusion33 = damage_phaseField_getDamageDiffusion33(ipc, ip, el)
case default
constitutive_getDamageDiffusion33 = lattice_DamageDiffusion33(1:3,1:3,material_phase(ipc,ip,el))
end select

134
code/damage_anisoBrittle.f90
View File

@ -257,9 +257,8 @@ subroutine damage_anisoBrittle_init(fileUnit)
sizeDotState = 0_pInt
sizeState = sizeDotState + &
1_pInt + & ! non-local damage
damage_anisoBrittle_totalNcleavage(instance) + & ! opening on each damage system
damage_anisoBrittle_totalNcleavage(instance) + & ! damage on each damage system
9_pInt ! Fd
1_pInt + & ! opening on each damage system
9_pInt ! Fd
damageState(phase)%sizeState = sizeState
damageState(phase)%sizeDotState = sizeDotState
@ -304,16 +303,12 @@ subroutine damage_anisoBrittle_stateInit(phase,instance)
real(pReal), dimension(damageState(phase)%sizeState) :: tempState
tempState(1) = 1.0_pReal
tempState(2 : &
1 + damage_anisoBrittle_totalNcleavage(instance)) = 0.0_pReal
tempState(2 + damage_anisoBrittle_totalNcleavage(instance): &
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,:)))
damageState(phase)%state0 = damageState(phase)%state
damageState(phase)%partionedState0 = damageState(phase)%state
tempState(1) = 1.0_pReal
tempState(2) = 0.0_pReal
tempState(3:11) = reshape(math_I3, shape=[9])
damageState(phase)%state0 = spread(tempState,2,size(damageState(phase)%state(1,:)))
end subroutine damage_anisoBrittle_stateInit
!--------------------------------------------------------------------------------------------------
@ -329,20 +324,20 @@ subroutine damage_anisoBrittle_aTolState(phase,instance)
instance ! number specifying the current instance of the damage
real(pReal), dimension(damageState(phase)%sizeState) :: tempTol
tempTol(1) = damage_anisoBrittle_aTol_damage(instance)
tempTol(2 : &
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_damage(instance)
tempTol(2 +2*damage_anisoBrittle_totalNcleavage(instance): &
10+2*damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_damage(instance)
tempTol(1) = damage_anisoBrittle_aTol_damage(instance)
tempTol(2) = damage_anisoBrittle_aTol_disp (instance)
tempTol(3:11) = damage_anisoBrittle_aTol_damage(instance)
damageState(phase)%aTolState = tempTol
end subroutine damage_anisoBrittle_aTolState
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine damage_anisoBrittle_microstructure(Tstar_v, subdt, ipc, ip, el)
use numerics, only: &
residualStiffness
use material, only: &
mappingConstitutive, &
phase_damageInstance, &
@ -367,13 +362,9 @@ subroutine damage_anisoBrittle_microstructure(Tstar_v, subdt, ipc, ip, el)
phase, &
constituent, &
instance, &
f, i, index_d, index_o, index_myFamily
f, i, index_myFamily
real(pReal) :: &
localDamage, &
resForce, &
deltaState, &
available, &
predicted
localDamage
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit
@ -381,55 +372,29 @@ subroutine damage_anisoBrittle_microstructure(Tstar_v, subdt, ipc, ip, el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase)
localDamage = minval(damageState(phase)%state(2+ damage_anisoBrittle_totalNcleavage(instance): &
1+2*damage_anisoBrittle_totalNcleavage(instance),constituent))
index_o = 2_pInt
index_d = 2_pInt + damage_anisoBrittle_totalNcleavage(instance)
damageState(phase)%state(2,constituent) = damageState(phase)%subState0(2,constituent)
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
traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
traction_crit = damage_anisoBrittle_critLoad(f,instance)
deltaState = subdt*damage_anisoBrittle_sdot_0(instance)* &
( (abs(traction_d)/traction_crit)**damage_anisoBrittle_N(instance) + &
(abs(traction_t)/traction_crit)**damage_anisoBrittle_N(instance) + &
(max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance))
if (damageState(phase)%state(index_o,constituent) <= 1.0_pReal) then
damageState(phase)%state(index_o,constituent) = &
damageState(phase)%subState0(index_o,constituent) + &
deltaState
else
if (damageState(phase)%state(index_d,constituent) == localDamage) then
resForce = damage_anisoBrittle_getDamage(ipc, ip, el) - localDamage
else
resForce = 0.0_pReal
endif
available = (damageState(phase)%subState0(index_o,constituent) - resForce)** &
(1.0_pReal - damage_anisoBrittle_N(instance))
predicted = -(1.0_pReal - damage_anisoBrittle_N(instance))*deltaState
if (available > predicted) then
damageState(phase)%state(index_o,constituent) = resForce + &
(available - predicted)**(1.0_pReal/(1.0_pReal - damage_anisoBrittle_N(instance)))
else
damageState(phase)%state(index_o,constituent) = huge(1.0_pReal)
endif
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) - resForce)))
endif
traction_crit = damage_anisoBrittle_critLoad(f,instance)* &
damage_anisoBrittle_getDamage(ipc, ip, el)
damageState(phase)%state(2,constituent) = &
damageState(phase)%state(2,constituent) + &
subdt*damage_anisoBrittle_sdot_0(instance)* &
( (abs(traction_d)/traction_crit)**damage_anisoBrittle_N(instance) + &
(abs(traction_t)/traction_crit)**damage_anisoBrittle_N(instance) + &
(max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance))/ &
damage_anisoBrittle_critDisp(f,instance)
index_d = index_d + 1_pInt; index_o = index_o + 1_pInt
enddo
enddo
localDamage = minval(damageState(phase)%state(2+ damage_anisoBrittle_totalNcleavage(instance): &
1+2*damage_anisoBrittle_totalNcleavage(instance),constituent))
localDamage = min(1.0_pReal,max(residualStiffness,1.0/damageState(phase)%state(2,constituent)))
damageState(phase)%state(1,constituent) = &
localDamage + &
(damageState(phase)%subState0(1,constituent) - localDamage)* &
@ -441,8 +406,6 @@ end subroutine damage_anisoBrittle_microstructure
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el)
use numerics, only: &
residualStiffness
use material, only: &
mappingConstitutive, &
phase_damageInstance, &
@ -472,7 +435,7 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
phase, &
constituent, &
instance, &
f, i, index, index_myFamily, k, l, m, n
f, i, index_myFamily, k, l, m, n
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
@ -483,7 +446,6 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
Ld = 0.0_pReal
dLd_dTstar3333 = 0.0_pReal
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
@ -491,7 +453,7 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
traction_crit = damage_anisoBrittle_critLoad(f,instance)* &
(damageState(phase)%state0(index,constituent) + residualStiffness)
damage_anisoBrittle_getDamage(ipc, ip, el)
udotd = &
sign(1.0_pReal,traction_d)* &
damage_anisoBrittle_sdot_0(instance)* &
@ -530,7 +492,6 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
lattice_Scleavage(m,n,3,index_myFamily+i,phase)
endif
index = index + 1_pInt
enddo
enddo
dLd_dTstar = math_Plain3333to99(dLd_dTstar3333)
@ -562,11 +523,7 @@ pure function damage_anisoBrittle_getFd(ipc, ip, el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase)
damage_anisoBrittle_getFd = &
reshape(damageState(phase)% &
state(2*damage_anisoBrittle_totalNcleavage(instance)+2: &
2*damage_anisoBrittle_totalNcleavage(instance)+10,constituent), &
shape=[3,3])
damage_anisoBrittle_getFd = reshape(damageState(phase)%state(3:11,constituent),shape=[3,3])
end function damage_anisoBrittle_getFd
@ -574,8 +531,6 @@ end function damage_anisoBrittle_getFd
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el)
use numerics, only: &
residualStiffness
use material, only: &
mappingConstitutive, &
phase_damageInstance, &
@ -605,7 +560,7 @@ subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el)
phase, &
constituent, &
instance, &
f, i, index, index_myFamily
f, i, index_myFamily
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, udott, udotn
@ -615,7 +570,6 @@ subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el)
instance = phase_damageInstance(phase)
Ld = 0.0_pReal
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
@ -623,7 +577,7 @@ subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el)
traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
traction_crit = damage_anisoBrittle_critLoad(f,instance)* &
(damageState(phase)%state0(index,constituent) + residualStiffness)
damage_anisoBrittle_getDamage(ipc, ip, el)
udotd = &
sign(1.0_pReal,traction_d)* &
damage_anisoBrittle_sdot_0(instance)* &
@ -641,11 +595,9 @@ subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el)
(max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance)
Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
index = index + 1_pInt
enddo
enddo
damageState(phase)%state(2*damage_anisoBrittle_totalNcleavage(instance)+2: &
2*damage_anisoBrittle_totalNcleavage(instance)+10,constituent) = &
damageState(phase)%state(3:11,constituent) = &
reshape(math_mul33x33(math_inv33(math_I3 - dt*Ld), &
damage_anisoBrittle_getFd0(ipc, ip, el)), shape=[9])
@ -676,11 +628,7 @@ pure function damage_anisoBrittle_getFd0(ipc, ip, el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase)
damage_anisoBrittle_getFd0 = &
reshape(damageState(phase)% &
subState0(2*damage_anisoBrittle_totalNcleavage(instance)+2: &
2*damage_anisoBrittle_totalNcleavage(instance)+10,constituent), &
shape=[3,3])
damage_anisoBrittle_getFd0 = reshape(damageState(phase)%subState0(3:11,constituent),shape=[3,3])
end function damage_anisoBrittle_getFd0
@ -710,10 +658,7 @@ pure function damage_anisoBrittle_getPartionedFd0(ipc, ip, el)
instance = phase_damageInstance(phase)
damage_anisoBrittle_getPartionedFd0 = &
reshape(damageState(phase)% &
partionedState0(2*damage_anisoBrittle_totalNcleavage(instance)+2: &
2*damage_anisoBrittle_totalNcleavage(instance)+10,constituent), &
shape=[3,3])
reshape(damageState(phase)%partionedState0(3:11,constituent),shape=[3,3])
end function damage_anisoBrittle_getPartionedFd0
@ -741,7 +686,7 @@ function damage_anisoBrittle_getDamage(ipc, ip, el)
select case(field_damage_type(material_homog(ip,el)))
case (FIELD_DAMAGE_LOCAL_ID)
damage_anisoBrittle_getDamage = damageState(mappingConstitutive(2,ipc,ip,el))% &
state0(1,mappingConstitutive(1,ipc,ip,el))
state(1,mappingConstitutive(1,ipc,ip,el))
case (FIELD_DAMAGE_NONLOCAL_ID)
damage_anisoBrittle_getDamage = fieldDamage(material_homog(ip,el))% &
@ -799,7 +744,8 @@ end function damage_anisoBrittle_getLocalDamage
function damage_anisoBrittle_postResults(ipc,ip,el)
use material, only: &
mappingConstitutive, &
phase_damageInstance
phase_damageInstance, &
damageState
implicit none
integer(pInt), intent(in) :: &

157
code/damage_anisoDuctile.f90
View File

@ -66,7 +66,6 @@ module damage_anisoDuctile
damage_anisoDuctile_getDamage, &
damage_anisoDuctile_putLocalDamage, &
damage_anisoDuctile_getLocalDamage, &
damage_anisoDuctile_getSlipDamage, &
damage_anisoDuctile_postResults
contains
@ -254,7 +253,7 @@ subroutine damage_anisoDuctile_init(fileUnit)
sizeDotState = 0_pInt
sizeState = sizeDotState + &
1_pInt + & ! time regularised damage
damage_anisoDuctile_totalNslip(instance) + & ! slip system damages
1_pInt + & ! damaged plasticity
9 ! Fd
damageState(phase)%sizeState = sizeState
damageState(phase)%sizeDotState = sizeDotState
@ -297,16 +296,11 @@ subroutine damage_anisoDuctile_stateInit(phase, instance)
real(pReal), dimension(damageState(phase)%sizeState) :: tempState
tempState(1) = 1.0_pReal
tempState(2 : &
1 + damage_anisoDuctile_totalNslip(instance)) = 1.0_pReal
tempState(damage_anisoDuctile_totalNslip(instance)+2: &
damage_anisoDuctile_totalNslip(instance)+10) = reshape(math_I3, shape=[9])
tempState(1) = 1.0_pReal
tempState(2) = 0.0_pReal
tempState(3:11) = reshape(math_I3, shape=[9])
damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:)))
damageState(phase)%state0 = damageState(phase)%state
damageState(phase)%partionedState0 = damageState(phase)%state
damageState(phase)%state0 = spread(tempState,2,size(damageState(phase)%state(1,:)))
end subroutine damage_anisoDuctile_stateInit
@ -325,20 +319,22 @@ subroutine damage_anisoDuctile_aTolState(phase,instance)
tempTol = damage_anisoDuctile_aTol_damage(instance)
damageState(phase)%aTolState = tempTol
end subroutine damage_anisoDuctile_aTolState
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine damage_anisoDuctile_microstructure(subdt, ipc, ip, el)
use numerics, only: &
residualStiffness
use material, only: &
mappingConstitutive, &
phase_damageInstance, &
plasticState, &
damageState
use lattice, only: &
lattice_maxNslipFamily
use lattice, only: &
lattice_maxNslipFamily, &
lattice_DamageMobility
implicit none
@ -354,36 +350,30 @@ subroutine damage_anisoDuctile_microstructure(subdt, ipc, ip, el)
instance, &
index, f, i
real(pReal) :: &
localDamage, &
drivingForce
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))
index = 1_pInt
damageState(phase)%state(2,constituent) = damageState(phase)%subState0(2,constituent)
do f = 1_pInt,lattice_maxNslipFamily
do i = 1_pInt,damage_anisoDuctile_Nslip(f,instance) ! process each (active) slip system in family
if (localDamage == damageState(phase)%state(index+1,constituent)) then
drivingForce = plasticState(phase)%accumulatedSlip(index,constituent)/damage_anisoDuctile_critPlasticStrain(f,instance) - &
damage_anisoDuctile_getDamage(ipc, ip, el)
damageState(phase)%state(index+1,constituent) = &
min(damageState(phase)%state0(index+1,constituent), &
(sqrt(drivingForce*drivingForce + 4.0_pReal) - drivingForce)/2.0_pReal)
else
drivingForce = plasticState(phase)%accumulatedSlip(index,constituent)/damage_anisoDuctile_critPlasticStrain(f,instance)
damageState(phase)%state(index+1,constituent) = &
min(damageState(phase)%state0(index+1,constituent), &
1.0_pReal/drivingForce) ! irreversibility
endif
do i = 1_pInt,damage_anisoDuctile_Nslip(f,instance) ! process each (active) slip system in family
damageState(phase)%state(2,constituent) = &
damageState(phase)%state(2,constituent) + &
subdt* &
plasticState(phase)%slipRate(index,constituent)/ &
(damage_anisoDuctile_getDamage(ipc, ip, el)**damage_anisoDuctile_N(instance))/ &
damage_anisoDuctile_critPlasticStrain(f,instance)
index = index + 1_pInt
enddo
enddo
localDamage = minval(damageState(phase)%state(2:1+damage_anisoDuctile_totalNslip(instance),constituent))
localDamage = &
max(residualStiffness,min(1.0_pReal,1.0_pReal/damageState(phase)%state(2,constituent)))
damageState(phase)%state(1,constituent) = &
localDamage + &
(damageState(phase)%subState0(1,constituent) - localDamage)* &
@ -395,8 +385,6 @@ end subroutine damage_anisoDuctile_microstructure
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el)
use numerics, only: &
residualStiffness
use lattice, only: &
lattice_maxNslipFamily, &
lattice_NslipSystem, &
@ -436,7 +424,7 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
phase, &
constituent, &
instance, &
f, i, index, index_myFamily, k, l, m, n
f, i, index_myFamily, k, l, m, n
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
@ -447,13 +435,9 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
Ld = 0.0_pReal
dLd_dTstar3333 = 0.0_pReal
index = 2_pInt
do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,damage_anisoDuctile_Nslip(f,instance) ! process each (active) slip system in family
projection_d = math_tensorproduct(lattice_sd(1:3,index_myFamily+i,phase),&
lattice_sn(1:3,index_myFamily+i,phase))
projection_t = math_tensorproduct(lattice_st(1:3,index_myFamily+i,phase),&
@ -470,12 +454,13 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
traction_n = dot_product(Tstar_v,projection_n_v(1:6))
traction_crit = damage_anisoDuctile_critLoad(f,instance)* &
(damageState(phase)%state0(index,constituent) + residualStiffness) ! degrading critical load carrying capacity by damage
damage_anisoDuctile_getDamage(ipc, ip, el) ! degrading critical load carrying capacity by damage
udotd = &
sign(1.0_pReal,traction_d)* &
damage_anisoDuctile_sdot_0(instance)* &
(abs(traction_d)/traction_crit)**damage_anisoDuctile_N(instance)
(abs(traction_d)/traction_crit - &
abs(traction_d)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance)
if (udotd /= 0.0_pReal) then
Ld = Ld + udotd*projection_d
dudotd_dt = udotd*damage_anisoDuctile_N(instance)/traction_d
@ -487,7 +472,8 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
udott = &
sign(1.0_pReal,traction_t)* &
damage_anisoDuctile_sdot_0(instance)* &
(abs(traction_t)/traction_crit)**damage_anisoDuctile_N(instance)
(abs(traction_t)/traction_crit - &
abs(traction_t)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance)
if (udott /= 0.0_pReal) then
Ld = Ld + udott*projection_t
dudott_dt = udott*damage_anisoDuctile_N(instance)/traction_t
@ -497,15 +483,15 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
endif
udotn = &
damage_anisoDuctile_sdot_0(instance)* &
(max(0.0_pReal,traction_n)/traction_crit)**damage_anisoDuctile_N(instance)
(max(0.0_pReal,traction_n)/traction_crit - &
max(0.0_pReal,traction_n)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance)
if (udotn /= 0.0_pReal) then
Ld = Ld + udotn*projection_n(1:3,1:3)
dudotn_dt = udotn*damage_anisoDuctile_N(instance)/traction_n
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudotn_dt*projection_n(k,l)* projection_n(m,n)
endif
index = index + 1_pInt
endif
enddo
enddo
@ -539,10 +525,7 @@ pure function damage_anisoDuctile_getFd(ipc, ip, el)
instance = phase_damageInstance(phase)
damage_anisoDuctile_getFd = &
reshape(damageState(phase)% &
state(damage_anisoDuctile_totalNslip(instance)+2: &
damage_anisoDuctile_totalNslip(instance)+10,constituent), &
shape=[3,3])
reshape(damageState(phase)%state(3:11,constituent),shape=[3,3])
end function damage_anisoDuctile_getFd
@ -550,8 +533,6 @@ end function damage_anisoDuctile_getFd
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el)
use numerics, only: &
residualStiffness
use material, only: &
mappingConstitutive, &
phase_damageInstance, &
@ -559,9 +540,7 @@ subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el)
use math, only: &
math_mul33x33, &
math_inv33, &
math_Plain3333to99, &
math_I3, &
math_identity4th, &
math_symmetric33, &
math_Mandel33to6, &
math_tensorproduct
@ -591,7 +570,7 @@ subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el)
phase, &
constituent, &
instance, &
f, i, index, index_myFamily
f, i, index_myFamily
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, udott, udotn
@ -599,7 +578,7 @@ subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el)
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase)
index = 2_pInt
Ld = 0.0_pReal
do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family
@ -621,12 +600,13 @@ subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el)
traction_n = dot_product(Tstar_v,projection_n_v(1:6))
traction_crit = damage_anisoDuctile_critLoad(f,instance)* &
(damageState(phase)%state0(index,constituent) + residualStiffness) ! degrading critical load carrying capacity by damage
damage_anisoDuctile_getDamage(ipc, ip, el) ! degrading critical load carrying capacity by damage
udotd = &
sign(1.0_pReal,traction_d)* &
damage_anisoDuctile_sdot_0(instance)* &
(abs(traction_d)/traction_crit)**damage_anisoDuctile_N(instance)
(abs(traction_d)/traction_crit - &
abs(traction_d)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance)
if (udotd /= 0.0_pReal) then
Ld = Ld + udotd*projection_d
endif
@ -634,26 +614,25 @@ subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el)
udott = &
sign(1.0_pReal,traction_t)* &
damage_anisoDuctile_sdot_0(instance)* &
(abs(traction_t)/traction_crit)**damage_anisoDuctile_N(instance)
(abs(traction_t)/traction_crit) - &
abs(traction_t)/damage_anisoDuctile_critLoad(f,instance)**damage_anisoDuctile_N(instance)
if (udott /= 0.0_pReal) then
Ld = Ld + udott*projection_t
endif
udotn = &
damage_anisoDuctile_sdot_0(instance)* &
(max(0.0_pReal,traction_n)/traction_crit)**damage_anisoDuctile_N(instance)
(max(0.0_pReal,traction_n)/traction_crit - &
max(0.0_pReal,traction_n)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance)
if (udotn /= 0.0_pReal) then
Ld = Ld + udotn*projection_n(1:3,1:3)
endif
index = index + 1_pInt
enddo
enddo
damageState(phase)%state(damage_anisoDuctile_totalNslip(instance)+2: &
damage_anisoDuctile_totalNslip(instance)+10,constituent) = &
reshape(math_mul33x33(math_inv33(math_I3 - dt*Ld), &
damage_anisoDuctile_getFd0(ipc, ip, el)), shape=[9])
damageState(phase)%state(3:11,constituent) = reshape(math_mul33x33(math_inv33(math_I3 - dt*Ld), &
damage_anisoDuctile_getFd0(ipc, ip, el)), &
shape=[9])
end subroutine damage_anisoDuctile_putFd
@ -683,10 +662,7 @@ pure function damage_anisoDuctile_getFd0(ipc, ip, el)
instance = phase_damageInstance(phase)
damage_anisoDuctile_getFd0 = &
reshape(damageState(phase)% &
subState0(damage_anisoDuctile_totalNslip(instance)+2: &
damage_anisoDuctile_totalNslip(instance)+10,constituent), &
shape=[3,3])
reshape(damageState(phase)%subState0(3:11,constituent),shape=[3,3])
end function damage_anisoDuctile_getFd0
@ -716,10 +692,7 @@ pure function damage_anisoDuctile_getPartionedFd0(ipc, ip, el)
instance = phase_damageInstance(phase)
damage_anisoDuctile_getPartionedFd0 = &
reshape(damageState(phase)% &
partionedState0(damage_anisoDuctile_totalNslip(instance)+2: &
damage_anisoDuctile_totalNslip(instance)+10,constituent), &
shape=[3,3])
reshape(damageState(phase)%partionedState0(3:11,constituent),shape=[3,3])
end function damage_anisoDuctile_getPartionedFd0
@ -747,7 +720,7 @@ function damage_anisoDuctile_getDamage(ipc, ip, el)
select case(field_damage_type(material_homog(ip,el)))
case (FIELD_DAMAGE_LOCAL_ID)
damage_anisoDuctile_getDamage = damageState(mappingConstitutive(2,ipc,ip,el))% &
state0(1,mappingConstitutive(1,ipc,ip,el))
state(1,mappingConstitutive(1,ipc,ip,el))
case (FIELD_DAMAGE_NONLOCAL_ID)
damage_anisoDuctile_getDamage = fieldDamage(material_homog(ip,el))% &
@ -799,40 +772,6 @@ function damage_anisoDuctile_getLocalDamage(ipc, ip, el)
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, &
damageState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal) :: &
damage_anisoDuctile_getSlipDamage(damage_anisoDuctile_totalNslip( &
phase_damageInstance(mappingConstitutive(2,ipc,ip,el))))
integer(pInt) :: &
phase, &
constituent, &
instance
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase)
damage_anisoDuctile_getSlipDamage = &
damageState(phase)%state0(2:1+damage_anisoDuctile_totalNslip(instance),constituent) + &
residualStiffness
end function damage_anisoDuctile_getSlipDamage
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------

55
code/damage_isoDuctile.f90
View File

@ -27,7 +27,8 @@ module damage_isoDuctile
real(pReal), dimension(:), allocatable, private :: &
damage_isoDuctile_aTol, &
damage_isoDuctile_critPlasticStrain
damage_isoDuctile_critPlasticStrain, &
damage_isoDuctile_N
enum, bind(c)
enumerator :: undefined_ID, &
@ -44,7 +45,6 @@ module damage_isoDuctile
damage_isoDuctile_aTolState, &
damage_isoDuctile_microstructure, &
damage_isoDuctile_getDamage, &
damage_isoDuctile_getSlipDamage, &
damage_isoDuctile_putLocalDamage, &
damage_isoDuctile_getLocalDamage, &
damage_isoDuctile_getDamagedC66, &
@ -121,6 +121,7 @@ subroutine damage_isoDuctile_init(fileUnit)
allocate(damage_isoDuctile_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
allocate(damage_isoDuctile_Noutput(maxNinstance), source=0_pInt)
allocate(damage_isoDuctile_critPlasticStrain(maxNinstance), source=0.0_pReal)
allocate(damage_isoDuctile_N(maxNinstance), source=0.0_pReal)
allocate(damage_isoDuctile_aTol(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
@ -157,6 +158,9 @@ subroutine damage_isoDuctile_init(fileUnit)
case ('criticalplasticstrain')
damage_isoDuctile_critPlasticStrain(instance) = IO_floatValue(line,positions,2_pInt)
case ('damageratesensitivity')
damage_isoDuctile_N(instance) = IO_floatValue(line,positions,2_pInt)
case ('atol_damage')
damage_isoDuctile_aTol(instance) = IO_floatValue(line,positions,2_pInt)
@ -239,10 +243,11 @@ subroutine damage_isoDuctile_stateInit(phase)
real(pReal), dimension(damageState(phase)%sizeState) :: tempState
tempState = 1.0_pReal
damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:)))
damageState(phase)%state0 = damageState(phase)%state
damageState(phase)%partionedState0 = damageState(phase)%state
tempState(1) = 1.0_pReal
tempState(2) = 0.0_pReal
damageState(phase)%state0 = spread(tempState,2,size(damageState(phase)%state(1,:)))
end subroutine damage_isoDuctile_stateInit
!--------------------------------------------------------------------------------------------------
@ -260,6 +265,7 @@ subroutine damage_isoDuctile_aTolState(phase,instance)
tempTol = damage_isoDuctile_aTol(instance)
damageState(phase)%aTolState = tempTol
end subroutine damage_isoDuctile_aTolState
!--------------------------------------------------------------------------------------------------
@ -281,24 +287,30 @@ subroutine damage_isoDuctile_microstructure(subdt,ipc, ip, el)
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in) :: &
real(pReal), intent(in) :: &
subdt
integer(pInt) :: &
phase, constituent, instance
real(pReal) :: &
localDamage
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase)
damageState(phase)%state(2,constituent) = &
max(residualStiffness, &
min(damageState(phase)%state0(2,constituent), &
damage_isoDuctile_critPlasticStrain(instance)/ &
sum(plasticState(phase)%accumulatedSlip(:,constituent))))
damageState(phase)%subState0(2,constituent) + &
subdt* &
sum(plasticState(phase)%slipRate(:,constituent))/ &
(damage_isoDuctile_getDamage(ipc, ip, el)**damage_isoDuctile_N(instance))/ &
damage_isoDuctile_critPlasticStrain(instance)
localDamage = &
max(residualStiffness,min(1.0_pReal, 1.0_pReal/damageState(phase)%state(2,constituent)))
damageState(phase)%state(1,constituent) = &
damageState(phase)%state(2,constituent) + &
(damageState(phase)%subState0(1,constituent) - damageState(phase)%state(2,constituent))* &
localDamage + &
(damageState(phase)%subState0(1,constituent) - localDamage)* &
exp(-subdt/lattice_DamageMobility(phase))
end subroutine damage_isoDuctile_microstructure
@ -337,23 +349,6 @@ function damage_isoDuctile_getDamage(ipc, ip, el)
end function damage_isoDuctile_getDamage
!--------------------------------------------------------------------------------------------------
!> @brief returns slip damage
!--------------------------------------------------------------------------------------------------
function damage_isoDuctile_getSlipDamage(ipc, ip, el)
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal) :: damage_isoDuctile_getSlipDamage, damage
damage = damage_isoDuctile_getDamage(ipc, ip, el)
damage_isoDuctile_getSlipDamage = damage*damage
end function damage_isoDuctile_getSlipDamage
!--------------------------------------------------------------------------------------------------
!> @brief puts local damage
!--------------------------------------------------------------------------------------------------

16
code/plastic_disloKMC.f90
View File

@ -1112,8 +1112,7 @@ end subroutine plastic_disloKMC_microstructure
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_disloKMC_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,nSlipDamage,slipDamage, &
ipc,ip,el)
subroutine plastic_disloKMC_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el)
use prec, only: &
tol_math_check
use math, only: &
@ -1145,10 +1144,9 @@ subroutine plastic_disloKMC_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,
LATTICE_fcc_ID
implicit none
integer(pInt), intent(in) :: nSlipDamage,ipc,ip,el
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
real(pReal), dimension(nSlipDamage), intent(in) :: slipDamage
real(pReal), dimension(3,3), intent(out) :: Lp
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99
@ -1206,9 +1204,6 @@ subroutine plastic_disloKMC_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + plastic_disloKMC_nonSchmidCoeff(k,instance)*&
lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph)
enddo nonSchmidSystems
!* Applying damage to slip system
tau_slip_pos = tau_slip_pos/slipDamage(j)
tau_slip_neg = tau_slip_neg/slipDamage(j)
significantPostitiveStress: if((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then
!* Stress ratios
@ -1354,7 +1349,7 @@ end subroutine plastic_disloKMC_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_disloKMC_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage,ipc,ip,el)
subroutine plastic_disloKMC_dotState(Tstar_v,Temperature,ipc,ip,el)
use prec, only: &
tol_math_check
use math, only: &
@ -1385,12 +1380,9 @@ subroutine plastic_disloKMC_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage,
real(pReal), intent(in) :: &
temperature !< temperature at integration point
integer(pInt), intent(in) :: &
nSlipDamage, &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(nSlipDamage), intent(in) :: &
slipDamage
integer(pInt) :: instance,ns,nt,f,i,j,k,index_myFamily,s1,s2, &
ph, &
@ -1457,8 +1449,6 @@ subroutine plastic_disloKMC_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage,
tau_slip_neg = tau_slip_neg + plastic_disloKMC_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
enddo nonSchmidSystems
tau_slip_pos = tau_slip_pos/slipDamage(j)
tau_slip_neg = tau_slip_pos/slipDamage(j)
significantPositiveStress: if((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then
!* Stress ratios

16
code/plastic_disloUCLA.f90
View File

@ -1130,8 +1130,7 @@ end subroutine plastic_disloUCLA_microstructure
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,nSlipDamage,slipDamage, &
ipc,ip,el)
subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el)
use prec, only: &
tol_math_check
use math, only: &
@ -1163,10 +1162,9 @@ subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
LATTICE_fcc_ID
implicit none
integer(pInt), intent(in) :: nSlipDamage,ipc,ip,el
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
real(pReal), dimension(nSlipDamage), intent(in) :: slipDamage
real(pReal), dimension(3,3), intent(out) :: Lp
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99
@ -1224,9 +1222,6 @@ subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + plastic_disloUCLA_nonSchmidCoeff(k,instance)*&
lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph)
enddo nonSchmidSystems
!* Applying damage to slip system
tau_slip_pos = tau_slip_pos/slipDamage(j)
tau_slip_neg = tau_slip_neg/slipDamage(j)
significantPostitiveStress: if((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then
!* Stress ratio
@ -1437,7 +1432,7 @@ end subroutine plastic_disloUCLA_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage,ipc,ip,el)
subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el)
use prec, only: &
tol_math_check
use math, only: &
@ -1468,12 +1463,9 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage
real(pReal), intent(in) :: &
temperature !< temperature at integration point
integer(pInt), intent(in) :: &
nSlipDamage, &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(nSlipDamage), intent(in) :: &
slipDamage
integer(pInt) :: instance,ns,nt,f,i,j,k,index_myFamily,s1,s2, &
ph, &
@ -1539,8 +1531,6 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage
tau_slip_neg = tau_slip_neg + plastic_disloUCLA_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
enddo nonSchmidSystems
tau_slip_pos = tau_slip_pos/slipDamage(j)
tau_slip_neg = tau_slip_pos/slipDamage(j)
significantPositiveStress: if((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then
!* Stress ratios

14
code/plastic_dislotwin.f90
View File

@ -1363,7 +1363,7 @@ end subroutine plastic_dislotwin_microstructure
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,nSlipDamage,slipDamage,ipc,ip,el)
subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el)
use prec, only: &
tol_math_check
use math, only: &
@ -1399,10 +1399,9 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
LATTICE_fcc_ID
implicit none
integer(pInt), intent(in) :: nSlipDamage,ipc,ip,el
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
real(pReal), dimension(nSlipDamage), intent(in) :: slipDamage
real(pReal), dimension(3,3), intent(out) :: Lp
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99
@ -1461,7 +1460,7 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
!* Calculation of Lp
!* Resolved shear stress on slip system
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))/slipDamage(j)
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))
if((abs(tau_slip(j))-plasticState(ph)%state(7*ns+4*nt+2*nr+j, of)) > tol_math_check) then
!* Stress ratios
@ -1668,7 +1667,7 @@ end subroutine plastic_dislotwin_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage,ipc,ip,el)
subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
use prec, only: &
tol_math_check
use math, only: &
@ -1706,12 +1705,9 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage
real(pReal), intent(in) :: &
temperature !< temperature at integration point
integer(pInt), intent(in) :: &
nSlipDamage, &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(nSlipDamage), intent(in) :: &
slipDamage
integer(pInt) :: instance,ns,nt,nr,f,i,j,index_myFamily,s1,s2,a,b,sa,sb,ssa,ssb, &
ph, &
@ -1751,7 +1747,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage
j = j+1_pInt
!* Resolved shear stress on slip system
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))/slipDamage(j)
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))
if((abs(tau_slip(j))-plasticState(ph)%state(7*ns+4*nt+2*nr+j, of)) > tol_math_check) then
!* Stress ratios

14
code/plastic_j2.f90
View File

@ -362,7 +362,7 @@ end subroutine plastic_j2_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDamage,slipDamage,ipc,ip,el)
subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6, &
math_Mandel6to33, &
@ -388,12 +388,9 @@ subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDamage,slipDa
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
nSlipDamage, &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(nSlipDamage), intent(in) :: &
slipDamage
real(pReal), dimension(3,3) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
@ -417,7 +414,7 @@ subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDamage,slipDa
dLp_dTstar99 = 0.0_pReal
else
gamma_dot = plastic_j2_gdot0(instance) &
* (sqrt(1.5_pReal) * norm_Tstar_dev / (slipDamage(1)*plastic_j2_fTaylor(instance) * &
* (sqrt(1.5_pReal) * norm_Tstar_dev / (plastic_j2_fTaylor(instance) * &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el)))) &
**plastic_j2_n(instance)
@ -441,7 +438,7 @@ end subroutine plastic_j2_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_j2_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip,el)
subroutine plastic_j2_dotState(Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6
use mesh, only: &
@ -458,12 +455,9 @@ subroutine plastic_j2_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip,el)
real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
nSlipDamage, &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(nSlipDamage), intent(in) :: &
slipDamage
real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: &
@ -490,7 +484,7 @@ subroutine plastic_j2_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip,el)
! strain rate
gamma_dot = plastic_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!-----------------------------------------------------------------------------------
(slipDamage(1)*plastic_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**plastic_j2_n(instance)
(plastic_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**plastic_j2_n(instance)
!--------------------------------------------------------------------------------------------------
! hardening coefficient

19
code/plastic_nonlocal.f90
View File

@ -2021,8 +2021,7 @@ end subroutine plastic_nonlocal_kinetics
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, nSlipDamage, slipDamage, &
ipc, ip, el)
subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, ipc, ip, el)
use math, only: math_Plain3333to99, &
math_mul6x6, &
@ -2048,14 +2047,11 @@ use mesh, only: mesh_ipVolume
implicit none
!*** input variables
integer(pInt), intent(in) :: nSlipDamage, &
ipc, &
integer(pInt), intent(in) :: ipc, &
ip, & !< current integration point
el !< current element number
real(pReal), intent(in) :: Temperature !< temperature
real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(nSlipDamage), intent(in) :: &
slipDamage
!*** output variables
@ -2121,7 +2117,7 @@ tauThreshold = plasticState(ph)%state(iTauF(1:ns,instance),of)
do s = 1_pInt,ns
sLattice = slipSystemLattice(s,instance)
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph))/slipDamage(s)
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph))
tauNS(s,1) = tau(s)
tauNS(s,2) = tau(s)
if (tau(s) > 0.0_pReal) then
@ -2408,7 +2404,7 @@ end subroutine plastic_nonlocal_deltaState
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature,nSlipDamage,slipDamage, &
subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, &
timestep,subfrac, ip,el)
use prec, only: DAMASK_NaN
@ -2463,8 +2459,7 @@ implicit none
!*** input variables
integer(pInt), intent(in) :: ip, & !< current integration point
el, & !< current element number
nSlipDamage
el !< current element number
real(pReal), intent(in) :: Temperature, & !< temperature
timestep !< substepped crystallite time increment
real(pReal), dimension(6), intent(in) :: Tstar_v !< current 2nd Piola-Kirchhoff stress in Mandel notation
@ -2473,8 +2468,6 @@ real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), in
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe, & !< elastic deformation gradient
Fp !< plastic deformation gradient
real(pReal), dimension(nSlipDamage), intent(in) :: &
slipDamage
!*** local variables
@ -2639,7 +2632,7 @@ forall (t = 1_pInt:4_pInt) &
do s = 1_pInt,ns ! loop over slip systems
sLattice = slipSystemLattice(s,instance)
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph))/slipDamage(s) + tauBack(s)
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph)) + tauBack(s)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo

19
code/plastic_phenopowerlaw.f90
View File

@ -717,8 +717,7 @@ end subroutine plastic_phenopowerlaw_aTolState
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDamage,slipDamage, &
ipc,ip,el)
subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use math, only: &
math_Plain3333to99, &
math_Mandel6to33
@ -745,14 +744,11 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDa
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
integer(pInt), intent(in) :: &
nSlipDamage, &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), dimension(nSlipDamage), intent(in) :: &
slipDamage
integer(pInt) :: &
instance, &
@ -807,11 +803,11 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,nSlipDa
lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph)
enddo
gdot_slip_pos = 0.5_pReal*plastic_phenopowerlaw_gdot0_slip(instance)* &
((abs(tau_slip_pos)/(slipDamage(j)*plasticState(ph)%state(j,of))) &
((abs(tau_slip_pos)/(plasticState(ph)%state(j,of))) &
**plastic_phenopowerlaw_n_slip(instance))*sign(1.0_pReal,tau_slip_pos)
gdot_slip_neg = 0.5_pReal*plastic_phenopowerlaw_gdot0_slip(instance)* &
((abs(tau_slip_neg)/(slipDamage(j)*plasticState(ph)%state(j,of))) &
((abs(tau_slip_neg)/(plasticState(ph)%state(j,of))) &
**plastic_phenopowerlaw_n_slip(instance))*sign(1.0_pReal,tau_slip_neg)
Lp = Lp + (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
@ -871,7 +867,7 @@ end subroutine plastic_phenopowerlaw_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip,el)
subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
use lattice, only: &
lattice_Sslip_v, &
lattice_Stwin_v, &
@ -891,12 +887,9 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip,
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
nSlipDamage, &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element !< microstructure state
real(pReal), dimension(nSlipDamage), intent(in) :: &
slipDamage
integer(pInt) :: &
instance,ph, &
@ -968,8 +961,8 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,nSlipDamage,slipDamage,ipc,ip,
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
enddo nonSchmidSystems
gdot_slip(j) = plastic_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos)/(slipDamage(j)*plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg)/(slipDamage(j)*plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance))&
((abs(tau_slip_pos)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos)
enddo slipSystems1
enddo slipFamilies1

8
code/plastic_titanmod.f90
View File

@ -1332,8 +1332,7 @@ end subroutine plastic_titanmod_microstructure
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,temperature,nSlipDamage,slipDamage, &
ipc,ip,el)
subroutine plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,temperature,ipc,ip,el)
use math, only: &
math_Plain3333to99, &
math_Mandel6to33
@ -1365,7 +1364,6 @@ subroutine plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,temperature,
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
integer(pInt), intent(in) :: &
nSlipDamage, &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
@ -1373,8 +1371,6 @@ subroutine plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,temperature,
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature !< temperature at IP
real(pReal), dimension(nSlipDamage), intent(in) :: &
slipDamage
integer(pInt) :: &
index_myFamily, instance, &
ns,nt, &
@ -1432,7 +1428,7 @@ subroutine plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,temperature,
!* Calculation of Lp
!* Resolved shear stress on slip system
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))/slipDamage(j)
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))
if(lattice_structure(ph)==LATTICE_hex_ID) then ! only for prismatic and pyr <a> systems in hex
screwvelocity_prefactor=plastic_titanmod_debyefrequency(instance)* &
plasticState(ph)%state(4_pInt*ns+nt+j, of)*(plastic_titanmod_burgersPerSlipSys(j,instance)/ &