Merge branch 'temperature-direct' into 'development'
use data from other physics directly See merge request damask/DAMASK!485
This commit is contained in:
Sharan Roongta
commit
b4f3ac4577
|
|
@ -155,7 +155,7 @@ module phase
|
|||
real(pReal), dimension(3,3) :: P
|
||||
end function phase_P
|
||||
|
||||
module function thermal_T(ph,en) result(T)
|
||||
pure module function thermal_T(ph,en) result(T)
|
||||
integer, intent(in) :: ph,en
|
||||
real(pReal) :: T
|
||||
end function thermal_T
|
||||
|
|
|
|||
|
|
@ -73,47 +73,37 @@ submodule(phase:mechanical) plastic
|
|||
en
|
||||
end subroutine kinehardening_LpAndItsTangent
|
||||
|
||||
module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
|
||||
module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
|
||||
real(pReal), dimension(3,3), intent(out) :: &
|
||||
Lp
|
||||
real(pReal), dimension(3,3,3,3), intent(out) :: &
|
||||
dLp_dMp
|
||||
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
Mp
|
||||
real(pReal), intent(in) :: &
|
||||
T
|
||||
integer, intent(in) :: &
|
||||
ph, &
|
||||
en
|
||||
end subroutine dislotwin_LpAndItsTangent
|
||||
|
||||
pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
|
||||
pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
|
||||
real(pReal), dimension(3,3), intent(out) :: &
|
||||
Lp
|
||||
real(pReal), dimension(3,3,3,3), intent(out) :: &
|
||||
dLp_dMp
|
||||
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
Mp
|
||||
real(pReal), intent(in) :: &
|
||||
T
|
||||
integer, intent(in) :: &
|
||||
ph, &
|
||||
en
|
||||
end subroutine dislotungsten_LpAndItsTangent
|
||||
|
||||
module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
|
||||
Mp,Temperature,ph,en)
|
||||
module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
|
||||
real(pReal), dimension(3,3), intent(out) :: &
|
||||
Lp
|
||||
real(pReal), dimension(3,3,3,3), intent(out) :: &
|
||||
dLp_dMp
|
||||
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
Mp !< Mandel stress
|
||||
real(pReal), intent(in) :: &
|
||||
Temperature
|
||||
integer, intent(in) :: &
|
||||
ph, &
|
||||
en
|
||||
|
|
@ -282,13 +272,13 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
|
|||
call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
|
||||
|
||||
case (PLASTIC_NONLOCAL_ID) plasticType
|
||||
call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en)
|
||||
call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
|
||||
|
||||
case (PLASTIC_DISLOTWIN_ID) plasticType
|
||||
call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en)
|
||||
call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
|
||||
|
||||
case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
|
||||
call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en)
|
||||
call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
|
||||
|
||||
end select plasticType
|
||||
|
||||
|
|
|
|||
|
|
@ -257,27 +257,27 @@ end function plastic_dislotungsten_init
|
|||
!> @brief Calculate plastic velocity gradient and its tangent.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
|
||||
Mp,T,ph,en)
|
||||
Mp,ph,en)
|
||||
real(pReal), dimension(3,3), intent(out) :: &
|
||||
Lp !< plastic velocity gradient
|
||||
real(pReal), dimension(3,3,3,3), intent(out) :: &
|
||||
dLp_dMp !< derivative of Lp with respect to the Mandel stress
|
||||
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
Mp !< Mandel stress
|
||||
real(pReal), intent(in) :: &
|
||||
T !< temperature
|
||||
integer, intent(in) :: &
|
||||
ph, &
|
||||
en
|
||||
|
||||
integer :: &
|
||||
i,k,l,m,n
|
||||
real(pReal) :: &
|
||||
T !< temperature
|
||||
real(pReal), dimension(param(ph)%sum_N_sl) :: &
|
||||
dot_gamma_pos,dot_gamma_neg, &
|
||||
ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
|
||||
|
||||
|
||||
T = thermal_T(ph,en)
|
||||
Lp = 0.0_pReal
|
||||
dLp_dMp = 0.0_pReal
|
||||
|
||||
|
|
|
|||
|
|
@ -476,18 +476,18 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
|
|||
C66_tw, &
|
||||
C66_tr
|
||||
integer :: i
|
||||
real(pReal) :: f_unrotated
|
||||
real(pReal) :: f_matrix
|
||||
|
||||
|
||||
C = elastic_C66(ph,en)
|
||||
|
||||
associate(prm => param(ph), stt => state(ph))
|
||||
|
||||
f_unrotated = 1.0_pReal &
|
||||
- sum(stt%f_tw(1:prm%sum_N_tw,en)) &
|
||||
- sum(stt%f_tr(1:prm%sum_N_tr,en))
|
||||
f_matrix = 1.0_pReal &
|
||||
- sum(stt%f_tw(1:prm%sum_N_tw,en)) &
|
||||
- sum(stt%f_tr(1:prm%sum_N_tr,en))
|
||||
|
||||
homogenizedC = f_unrotated * C
|
||||
homogenizedC = f_matrix * C
|
||||
|
||||
twinActive: if (prm%sum_N_tw > 0) then
|
||||
C66_tw = lattice_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph))
|
||||
|
|
@ -513,20 +513,20 @@ end function plastic_dislotwin_homogenizedC
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Calculate plastic velocity gradient and its tangent.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
|
||||
module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
|
||||
|
||||
real(pReal), dimension(3,3), intent(out) :: Lp
|
||||
real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp
|
||||
real(pReal), dimension(3,3), intent(in) :: Mp
|
||||
integer, intent(in) :: ph,en
|
||||
real(pReal), intent(in) :: T
|
||||
|
||||
integer :: i,k,l,m,n
|
||||
real(pReal) :: &
|
||||
f_unrotated,StressRatio_p,&
|
||||
E_kB_T, &
|
||||
ddot_gamma_dtau, &
|
||||
tau
|
||||
f_matrix,StressRatio_p,&
|
||||
E_kB_T, &
|
||||
ddot_gamma_dtau, &
|
||||
tau, &
|
||||
T
|
||||
real(pReal), dimension(param(ph)%sum_N_sl) :: &
|
||||
dot_gamma_sl,ddot_gamma_dtau_sl
|
||||
real(pReal), dimension(param(ph)%sum_N_tw) :: &
|
||||
|
|
@ -556,69 +556,71 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
|
|||
0, 1, 1 &
|
||||
],pReal),[ 3,6])
|
||||
|
||||
associate(prm => param(ph), stt => state(ph))
|
||||
|
||||
f_unrotated = 1.0_pReal &
|
||||
- sum(stt%f_tw(1:prm%sum_N_tw,en)) &
|
||||
- sum(stt%f_tr(1:prm%sum_N_tr,en))
|
||||
|
||||
T = thermal_T(ph,en)
|
||||
Lp = 0.0_pReal
|
||||
dLp_dMp = 0.0_pReal
|
||||
|
||||
call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl)
|
||||
slipContribution: do i = 1, prm%sum_N_sl
|
||||
Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i)
|
||||
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
|
||||
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
|
||||
+ ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
|
||||
end do slipContribution
|
||||
associate(prm => param(ph), stt => state(ph))
|
||||
|
||||
call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
|
||||
twinContibution: do i = 1, prm%sum_N_tw
|
||||
Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
|
||||
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
|
||||
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
|
||||
+ ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
|
||||
end do twinContibution
|
||||
f_matrix = 1.0_pReal &
|
||||
- sum(stt%f_tw(1:prm%sum_N_tw,en)) &
|
||||
- sum(stt%f_tr(1:prm%sum_N_tr,en))
|
||||
|
||||
call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
|
||||
transContibution: do i = 1, prm%sum_N_tr
|
||||
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
|
||||
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
|
||||
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
|
||||
+ ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
|
||||
end do transContibution
|
||||
call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl)
|
||||
slipContribution: do i = 1, prm%sum_N_sl
|
||||
Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i)
|
||||
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
|
||||
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
|
||||
+ ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
|
||||
end do slipContribution
|
||||
|
||||
Lp = Lp * f_unrotated
|
||||
dLp_dMp = dLp_dMp * f_unrotated
|
||||
call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
|
||||
twinContibution: do i = 1, prm%sum_N_tw
|
||||
Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
|
||||
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
|
||||
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
|
||||
+ ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
|
||||
end do twinContibution
|
||||
|
||||
shearBandingContribution: if (dNeq0(prm%v_sb)) then
|
||||
call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
|
||||
transContibution: do i = 1, prm%sum_N_tr
|
||||
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
|
||||
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
|
||||
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
|
||||
+ ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
|
||||
end do transContibution
|
||||
|
||||
E_kB_T = prm%E_sb/(K_B*T)
|
||||
call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design?
|
||||
Lp = Lp * f_matrix
|
||||
dLp_dMp = dLp_dMp * f_matrix
|
||||
|
||||
do i = 1,6
|
||||
P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),&
|
||||
matmul(eigVectors,sb_mComposition(1:3,i)))
|
||||
tau = math_tensordot(Mp,P_sb)
|
||||
shearBandingContribution: if (dNeq0(prm%v_sb)) then
|
||||
|
||||
significantShearBandStress: if (abs(tau) > tol_math_check) then
|
||||
StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
|
||||
dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
|
||||
ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb &
|
||||
* (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
|
||||
* (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
|
||||
E_kB_T = prm%E_sb/(K_B*T)
|
||||
call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design?
|
||||
|
||||
Lp = Lp + dot_gamma_sb * P_sb
|
||||
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
|
||||
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
|
||||
+ ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n)
|
||||
end if significantShearBandStress
|
||||
end do
|
||||
do i = 1,6
|
||||
P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),&
|
||||
matmul(eigVectors,sb_mComposition(1:3,i)))
|
||||
tau = math_tensordot(Mp,P_sb)
|
||||
|
||||
end if shearBandingContribution
|
||||
significantShearBandStress: if (abs(tau) > tol_math_check) then
|
||||
StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
|
||||
dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
|
||||
ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb &
|
||||
* (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
|
||||
* (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
|
||||
|
||||
end associate
|
||||
Lp = Lp + dot_gamma_sb * P_sb
|
||||
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
|
||||
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
|
||||
+ ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n)
|
||||
end if significantShearBandStress
|
||||
end do
|
||||
|
||||
end if shearBandingContribution
|
||||
|
||||
end associate
|
||||
|
||||
end subroutine dislotwin_LpAndItsTangent
|
||||
|
||||
|
|
@ -638,7 +640,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
|
|||
|
||||
integer :: i
|
||||
real(pReal) :: &
|
||||
f_unrotated, &
|
||||
f_matrix, &
|
||||
d_hat, &
|
||||
v_cl, & !< climb velocity
|
||||
tau, &
|
||||
|
|
@ -661,9 +663,9 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
|
|||
mu = elastic_mu(ph,en)
|
||||
nu = elastic_nu(ph,en)
|
||||
|
||||
f_unrotated = 1.0_pReal &
|
||||
- sum(stt%f_tw(1:prm%sum_N_tw,en)) &
|
||||
- sum(stt%f_tr(1:prm%sum_N_tr,en))
|
||||
f_matrix = 1.0_pReal &
|
||||
- sum(stt%f_tw(1:prm%sum_N_tw,en)) &
|
||||
- sum(stt%f_tr(1:prm%sum_N_tr,en))
|
||||
|
||||
call kinetics_sl(Mp,T,ph,en,dot_gamma_sl)
|
||||
dot%gamma_sl(:,en) = abs(dot_gamma_sl)
|
||||
|
|
@ -709,10 +711,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
|
|||
- dot_rho_dip_climb
|
||||
|
||||
call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw)
|
||||
dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char
|
||||
dot%f_tw(:,en) = f_matrix*dot_gamma_tw/prm%gamma_char
|
||||
|
||||
call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr)
|
||||
dot%f_tr(:,en) = f_unrotated*dot_gamma_tr
|
||||
dot%f_tr(:,en) = f_matrix*dot_gamma_tr
|
||||
|
||||
end associate
|
||||
|
||||
|
|
|
|||
|
|
@ -741,7 +741,7 @@ end subroutine nonlocal_dependentState
|
|||
!> @brief calculates plastic velocity gradient and its tangent
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
|
||||
Mp,Temperature,ph,en)
|
||||
Mp,ph,en)
|
||||
real(pReal), dimension(3,3), intent(out) :: &
|
||||
Lp !< plastic velocity gradient
|
||||
real(pReal), dimension(3,3,3,3), intent(out) :: &
|
||||
|
|
@ -749,9 +749,6 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
|
|||
integer, intent(in) :: &
|
||||
ph, &
|
||||
en
|
||||
real(pReal), intent(in) :: &
|
||||
Temperature !< temperature
|
||||
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
Mp
|
||||
!< derivative of Lp with respect to Mp
|
||||
|
|
@ -771,67 +768,72 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
|
|||
real(pReal), dimension(param(ph)%sum_N_sl) :: &
|
||||
tau, & !< resolved shear stress including backstress terms
|
||||
dot_gamma !< shear rate
|
||||
real(pReal) :: &
|
||||
Temperature !< temperature
|
||||
|
||||
|
||||
Temperature = thermal_T(ph,en)
|
||||
Lp = 0.0_pReal
|
||||
dLp_dMp = 0.0_pReal
|
||||
|
||||
associate(prm => param(ph),dst=>dependentState(ph),stt=>state(ph))
|
||||
|
||||
!*** shortcut to state variables
|
||||
rho = getRho(ph,en)
|
||||
rhoSgl = rho(:,sgl)
|
||||
!*** shortcut to state variables
|
||||
rho = getRho(ph,en)
|
||||
rhoSgl = rho(:,sgl)
|
||||
|
||||
do s = 1,prm%sum_N_sl
|
||||
tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s))
|
||||
tauNS(s,1) = tau(s)
|
||||
tauNS(s,2) = tau(s)
|
||||
if (tau(s) > 0.0_pReal) then
|
||||
tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s))
|
||||
tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s))
|
||||
else
|
||||
tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s))
|
||||
tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s))
|
||||
end if
|
||||
end do
|
||||
tauNS = tauNS + spread(dst%tau_back(:,en),2,4)
|
||||
tau = tau + dst%tau_back(:,en)
|
||||
|
||||
! edges
|
||||
call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
|
||||
tau, tauNS(:,1), dst%tau_pass(:,en),1,Temperature, ph)
|
||||
v(:,2) = v(:,1)
|
||||
dv_dtau(:,2) = dv_dtau(:,1)
|
||||
dv_dtauNS(:,2) = dv_dtauNS(:,1)
|
||||
|
||||
!screws
|
||||
if (prm%nonSchmidActive) then
|
||||
do t = 3,4
|
||||
call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
|
||||
tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph)
|
||||
do s = 1,prm%sum_N_sl
|
||||
tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s))
|
||||
tauNS(s,1) = tau(s)
|
||||
tauNS(s,2) = tau(s)
|
||||
if (tau(s) > 0.0_pReal) then
|
||||
tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s))
|
||||
tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s))
|
||||
else
|
||||
tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s))
|
||||
tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s))
|
||||
end if
|
||||
end do
|
||||
else
|
||||
v(:,3:4) = spread(v(:,1),2,2)
|
||||
dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2)
|
||||
dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2)
|
||||
end if
|
||||
tauNS = tauNS + spread(dst%tau_back(:,en),2,4)
|
||||
tau = tau + dst%tau_back(:,en)
|
||||
|
||||
stt%v(:,en) = pack(v,.true.)
|
||||
! edges
|
||||
call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
|
||||
tau, tauNS(:,1), dst%tau_pass(:,en),1,Temperature, ph)
|
||||
v(:,2) = v(:,1)
|
||||
dv_dtau(:,2) = dv_dtau(:,1)
|
||||
dv_dtauNS(:,2) = dv_dtauNS(:,1)
|
||||
|
||||
!*** Bauschinger effect
|
||||
forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
|
||||
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
|
||||
!screws
|
||||
if (prm%nonSchmidActive) then
|
||||
do t = 3,4
|
||||
call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
|
||||
tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph)
|
||||
end do
|
||||
else
|
||||
v(:,3:4) = spread(v(:,1),2,2)
|
||||
dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2)
|
||||
dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2)
|
||||
end if
|
||||
|
||||
dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
|
||||
stt%v(:,en) = pack(v,.true.)
|
||||
|
||||
Lp = 0.0_pReal
|
||||
dLp_dMp = 0.0_pReal
|
||||
do s = 1,prm%sum_N_sl
|
||||
Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s)
|
||||
forall (i=1:3,j=1:3,k=1:3,l=1:3) &
|
||||
dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
|
||||
+ prm%P_sl(i,j,s) * prm%P_sl(k,l,s) &
|
||||
* sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) &
|
||||
+ prm%P_sl(i,j,s) &
|
||||
* (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
|
||||
- prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s)
|
||||
end do
|
||||
!*** Bauschinger effect
|
||||
forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
|
||||
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
|
||||
|
||||
dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
|
||||
|
||||
do s = 1,prm%sum_N_sl
|
||||
Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s)
|
||||
forall (i=1:3,j=1:3,k=1:3,l=1:3) &
|
||||
dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
|
||||
+ prm%P_sl(i,j,s) * prm%P_sl(k,l,s) &
|
||||
* sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) &
|
||||
+ prm%P_sl(i,j,s) &
|
||||
* (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
|
||||
- prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s)
|
||||
end do
|
||||
|
||||
end associate
|
||||
|
||||
|
|
@ -870,7 +872,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
|
|||
dUpper, & ! current maximum stable dipole distance for edges and screws
|
||||
dUpperOld, & ! old maximum stable dipole distance for edges and screws
|
||||
deltaDUpper ! change in maximum stable dipole distance for edges and screws
|
||||
|
||||
|
||||
|
||||
associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph))
|
||||
|
||||
mu = elastic_mu(ph,en)
|
||||
|
|
@ -1394,78 +1397,79 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
|
|||
belowThreshold
|
||||
type(rotation) :: mis
|
||||
|
||||
|
||||
associate(prm => param(ph))
|
||||
ns = prm%sum_N_sl
|
||||
ns = prm%sum_N_sl
|
||||
|
||||
en = material_phaseMemberAt(1,i,e)
|
||||
!*** start out fully compatible
|
||||
my_compatibility = 0.0_pReal
|
||||
forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal
|
||||
en = material_phaseMemberAt(1,i,e)
|
||||
!*** start out fully compatible
|
||||
my_compatibility = 0.0_pReal
|
||||
forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal
|
||||
|
||||
neighbors: do n = 1,nIPneighbors
|
||||
neighbor_e = IPneighborhood(1,n,i,e)
|
||||
neighbor_i = IPneighborhood(2,n,i,e)
|
||||
neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e)
|
||||
neighbor_phase = material_phaseAt(1,neighbor_e)
|
||||
neighbors: do n = 1,nIPneighbors
|
||||
neighbor_e = IPneighborhood(1,n,i,e)
|
||||
neighbor_i = IPneighborhood(2,n,i,e)
|
||||
neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e)
|
||||
neighbor_phase = material_phaseAt(1,neighbor_e)
|
||||
|
||||
if (neighbor_e <= 0 .or. neighbor_i <= 0) then
|
||||
!* FREE SURFACE
|
||||
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface)
|
||||
elseif (neighbor_phase /= ph) then
|
||||
!* PHASE BOUNDARY
|
||||
if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) &
|
||||
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal
|
||||
elseif (prm%chi_GB >= 0.0_pReal) then
|
||||
!* GRAIN BOUNDARY
|
||||
if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), &
|
||||
phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. &
|
||||
plasticState(neighbor_phase)%nonlocal) &
|
||||
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
|
||||
else
|
||||
!* GRAIN BOUNDARY ?
|
||||
!* Compatibility defined by relative orientation of slip systems:
|
||||
!* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection.
|
||||
!* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection.
|
||||
!* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one),
|
||||
!* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that
|
||||
!* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one.
|
||||
!* Finally the smallest compatibility value is decreased until the sum is exactly equal to one.
|
||||
!* All values below the threshold are set to zero.
|
||||
mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me))
|
||||
mySlipSystems: do s1 = 1,ns
|
||||
neighborSlipSystems: do s2 = 1,ns
|
||||
my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), &
|
||||
mis%rotate(prm%slip_normal(1:3,s2))) &
|
||||
* abs(math_inner(prm%slip_direction(1:3,s1), &
|
||||
mis%rotate(prm%slip_direction(1:3,s2))))
|
||||
my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), &
|
||||
mis%rotate(prm%slip_normal(1:3,s2)))) &
|
||||
* abs(math_inner(prm%slip_direction(1:3,s1), &
|
||||
mis%rotate(prm%slip_direction(1:3,s2))))
|
||||
end do neighborSlipSystems
|
||||
if (neighbor_e <= 0 .or. neighbor_i <= 0) then
|
||||
!* FREE SURFACE
|
||||
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface)
|
||||
elseif (neighbor_phase /= ph) then
|
||||
!* PHASE BOUNDARY
|
||||
if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) &
|
||||
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal
|
||||
elseif (prm%chi_GB >= 0.0_pReal) then
|
||||
!* GRAIN BOUNDARY
|
||||
if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), &
|
||||
phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. &
|
||||
plasticState(neighbor_phase)%nonlocal) &
|
||||
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
|
||||
else
|
||||
!* GRAIN BOUNDARY ?
|
||||
!* Compatibility defined by relative orientation of slip systems:
|
||||
!* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection.
|
||||
!* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection.
|
||||
!* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one),
|
||||
!* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that
|
||||
!* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one.
|
||||
!* Finally the smallest compatibility value is decreased until the sum is exactly equal to one.
|
||||
!* All values below the threshold are set to zero.
|
||||
mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me))
|
||||
mySlipSystems: do s1 = 1,ns
|
||||
neighborSlipSystems: do s2 = 1,ns
|
||||
my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), &
|
||||
mis%rotate(prm%slip_normal(1:3,s2))) &
|
||||
* abs(math_inner(prm%slip_direction(1:3,s1), &
|
||||
mis%rotate(prm%slip_direction(1:3,s2))))
|
||||
my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), &
|
||||
mis%rotate(prm%slip_normal(1:3,s2)))) &
|
||||
* abs(math_inner(prm%slip_direction(1:3,s1), &
|
||||
mis%rotate(prm%slip_direction(1:3,s2))))
|
||||
end do neighborSlipSystems
|
||||
|
||||
my_compatibilitySum = 0.0_pReal
|
||||
belowThreshold = .true.
|
||||
do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold))
|
||||
thresholdValue = maxval(my_compatibility(2,:,s1,n), belowThreshold) ! screws always positive
|
||||
nThresholdValues = real(count(my_compatibility(2,:,s1,n) >= thresholdValue),pReal)
|
||||
where (my_compatibility(2,:,s1,n) >= thresholdValue) belowThreshold = .false.
|
||||
if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
|
||||
where (abs(my_compatibility(:,:,s1,n)) >= thresholdValue) &
|
||||
my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,&
|
||||
my_compatibility(:,:,s1,n))
|
||||
my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue
|
||||
end do
|
||||
my_compatibilitySum = 0.0_pReal
|
||||
belowThreshold = .true.
|
||||
do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold))
|
||||
thresholdValue = maxval(my_compatibility(2,:,s1,n), belowThreshold) ! screws always positive
|
||||
nThresholdValues = real(count(my_compatibility(2,:,s1,n) >= thresholdValue),pReal)
|
||||
where (my_compatibility(2,:,s1,n) >= thresholdValue) belowThreshold = .false.
|
||||
if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
|
||||
where (abs(my_compatibility(:,:,s1,n)) >= thresholdValue) &
|
||||
my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,&
|
||||
my_compatibility(:,:,s1,n))
|
||||
my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue
|
||||
end do
|
||||
|
||||
where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal
|
||||
where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal
|
||||
where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal
|
||||
where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal
|
||||
|
||||
end do mySlipSystems
|
||||
end if
|
||||
end do mySlipSystems
|
||||
end if
|
||||
|
||||
end do neighbors
|
||||
end do neighbors
|
||||
|
||||
compatibility(:,:,:,:,i,e) = my_compatibility
|
||||
compatibility(:,:,:,:,i,e) = my_compatibility
|
||||
|
||||
end associate
|
||||
|
||||
|
|
@ -1646,53 +1650,55 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T,
|
|||
criticalStress_P, & !< maximum obstacle strength
|
||||
criticalStress_S !< maximum obstacle strength
|
||||
|
||||
associate(prm => param(ph))
|
||||
|
||||
v = 0.0_pReal
|
||||
dv_dtau = 0.0_pReal
|
||||
dv_dtauNS = 0.0_pReal
|
||||
|
||||
do s = 1,prm%sum_N_sl
|
||||
if (abs(tau(s)) > tauThreshold(s)) then
|
||||
associate(prm => param(ph))
|
||||
|
||||
!* Peierls contribution
|
||||
tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s))
|
||||
lambda_P = prm%b_sl(s)
|
||||
activationVolume_P = prm%w *prm%b_sl(s)**3
|
||||
criticalStress_P = prm%peierlsStress(s,c)
|
||||
activationEnergy_P = criticalStress_P * activationVolume_P
|
||||
tauRel_P = min(1.0_pReal, tauEff / criticalStress_P)
|
||||
tPeierls = 1.0_pReal / prm%nu_a &
|
||||
* exp(activationEnergy_P / (K_B * T) &
|
||||
* (1.0_pReal - tauRel_P**prm%p)**prm%q)
|
||||
dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) &
|
||||
* (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), &
|
||||
0.0_pReal, &
|
||||
tauEff < criticalStress_P)
|
||||
do s = 1,prm%sum_N_sl
|
||||
if (abs(tau(s)) > tauThreshold(s)) then
|
||||
|
||||
! Contribution from solid solution strengthening
|
||||
tauEff = abs(tau(s)) - tauThreshold(s)
|
||||
lambda_S = prm%b_sl(s) / sqrt(prm%c_sol)
|
||||
activationVolume_S = prm%f_sol * prm%b_sl(s)**3 / sqrt(prm%c_sol)
|
||||
criticalStress_S = prm%Q_sol / activationVolume_S
|
||||
tauRel_S = min(1.0_pReal, tauEff / criticalStress_S)
|
||||
tSolidSolution = 1.0_pReal / prm%nu_a &
|
||||
* exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
|
||||
dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) &
|
||||
* (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), &
|
||||
0.0_pReal, &
|
||||
tauEff < criticalStress_S)
|
||||
!* Peierls contribution
|
||||
tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s))
|
||||
lambda_P = prm%b_sl(s)
|
||||
activationVolume_P = prm%w *prm%b_sl(s)**3
|
||||
criticalStress_P = prm%peierlsStress(s,c)
|
||||
activationEnergy_P = criticalStress_P * activationVolume_P
|
||||
tauRel_P = min(1.0_pReal, tauEff / criticalStress_P)
|
||||
tPeierls = 1.0_pReal / prm%nu_a &
|
||||
* exp(activationEnergy_P / (K_B * T) &
|
||||
* (1.0_pReal - tauRel_P**prm%p)**prm%q)
|
||||
dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) &
|
||||
* (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), &
|
||||
0.0_pReal, &
|
||||
tauEff < criticalStress_P)
|
||||
|
||||
!* viscous glide velocity
|
||||
tauEff = abs(tau(s)) - tauThreshold(s)
|
||||
! Contribution from solid solution strengthening
|
||||
tauEff = abs(tau(s)) - tauThreshold(s)
|
||||
lambda_S = prm%b_sl(s) / sqrt(prm%c_sol)
|
||||
activationVolume_S = prm%f_sol * prm%b_sl(s)**3 / sqrt(prm%c_sol)
|
||||
criticalStress_S = prm%Q_sol / activationVolume_S
|
||||
tauRel_S = min(1.0_pReal, tauEff / criticalStress_S)
|
||||
tSolidSolution = 1.0_pReal / prm%nu_a &
|
||||
* exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
|
||||
dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) &
|
||||
* (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), &
|
||||
0.0_pReal, &
|
||||
tauEff < criticalStress_S)
|
||||
|
||||
!* viscous glide velocity
|
||||
tauEff = abs(tau(s)) - tauThreshold(s)
|
||||
|
||||
|
||||
v(s) = sign(1.0_pReal,tau(s)) &
|
||||
/ (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff))
|
||||
dv_dtau(s) = v(s)**2 * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2))
|
||||
dv_dtauNS(s) = v(s)**2 * dtPeierls_dtau / lambda_P
|
||||
v(s) = sign(1.0_pReal,tau(s)) &
|
||||
/ (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff))
|
||||
dv_dtau(s) = v(s)**2 * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2))
|
||||
dv_dtauNS(s) = v(s)**2 * dtPeierls_dtau / lambda_P
|
||||
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
end associate
|
||||
|
||||
|
|
|
|||
|
|
@ -271,7 +271,7 @@ end subroutine thermal_forward
|
|||
!----------------------------------------------------------------------------------------------
|
||||
!< @brief Get temperature (for use by non-thermal physics)
|
||||
!----------------------------------------------------------------------------------------------
|
||||
module function thermal_T(ph,en) result(T)
|
||||
pure module function thermal_T(ph,en) result(T)
|
||||
|
||||
integer, intent(in) :: ph, en
|
||||
real(pReal) :: T
|
||||
|
|
|
|||
Loading…
Reference in New Issue