preparing for temperature dependent C by calling it dynamically

This commit is contained in:
Sharan Roongta 2021-11-17 19:10:01 +01:00
parent 5a901a42ff
commit 55d6b1dd1a
5 changed files with 232 additions and 149 deletions

View File

@ -168,19 +168,19 @@ submodule(phase) mechanical
integer, intent(in) :: ph,en
end function plastic_dislotwin_homogenizedC
module function elastic_C66(ph) result(C66)
module function elastic_C66(ph,en) result(C66)
real(pReal), dimension(6,6) :: C66
integer, intent(in) :: ph
integer, intent(in) :: ph, en
end function elastic_C66
module function elastic_mu(ph) result(mu)
module function elastic_mu(ph,en) result(mu)
real(pReal) :: mu
integer, intent(in) :: ph
integer, intent(in) :: ph, en
end function elastic_mu
module function elastic_nu(ph) result(nu)
module function elastic_nu(ph,en) result(nu)
real(pReal) :: nu
integer, intent(in) :: ph
integer, intent(in) :: ph, en
end function elastic_nu
end interface

View File

@ -1,18 +1,22 @@
submodule(phase:mechanical) elastic
type :: tParameters
real(pReal), dimension(6,6) :: &
C66 = 0.0_pReal !< Elastic constants in Voigt notation
real(pReal) :: &
mu, &
nu
C_11, &
C_12, &
C_13, &
C_33, &
C_44, &
C_66
end type tParameters
type(tParameters), allocatable, dimension(:) :: param
contains
!--------------------------------------------------------------------------------------------------
!> @brief Initialize elasticity
!--------------------------------------------------------------------------------------------------
module subroutine elastic_init(phases)
class(tNode), pointer :: &
@ -41,28 +45,102 @@ module subroutine elastic_init(phases)
associate(prm => param(ph))
prm%C66(1,1) = elastic%get_asFloat('C_11')
prm%C66(1,2) = elastic%get_asFloat('C_12')
prm%C66(4,4) = elastic%get_asFloat('C_44')
prm%C_11 = elastic%get_asFloat('C_11')
prm%C_12 = elastic%get_asFloat('C_12')
prm%C_44 = elastic%get_asFloat('C_44')
if (any(phase_lattice(ph) == ['hP','tI'])) then
prm%C66(1,3) = elastic%get_asFloat('C_13')
prm%C66(3,3) = elastic%get_asFloat('C_33')
endif
if (phase_lattice(ph) == 'tI') prm%C66(6,6) = elastic%get_asFloat('C_66')
prm%C66 = lattice_symmetrize_C66(prm%C66,phase_lattice(ph))
prm%nu = lattice_equivalent_nu(prm%C66,'voigt')
prm%mu = lattice_equivalent_mu(prm%C66,'voigt')
prm%C66 = math_sym3333to66(math_Voigt66to3333(prm%C66)) ! Literature data is in Voigt notation
prm%C_13 = elastic%get_asFloat('C_13')
prm%C_33 = elastic%get_asFloat('C_33')
end if
if (phase_lattice(ph) == 'tI') prm%C_66 = elastic%get_asFloat('C_66')
end associate
enddo
end do
end subroutine elastic_init
!--------------------------------------------------------------------------------------------------
!> @brief returns 6x6 elasticity tensor
! internal function call to return dynamic values of the elasticity tensor
!--------------------------------------------------------------------------------------------------
function get_C66(ph,en)
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(6,6) :: get_C66
associate(prm => param(ph))
get_C66 = 0.0_pReal
get_C66(1,1) = prm%C_11
get_C66(1,2) = prm%C_12
get_C66(4,4) = prm%C_44
if (any(phase_lattice(ph) == ['hP','tI'])) then
get_C66(1,3) = prm%C_13
get_C66(3,3) = prm%C_33
end if
if (phase_lattice(ph) == 'tI') &
get_C66(6,6) = prm%C_66
get_C66 = lattice_symmetrize_C66(get_C66,phase_lattice(ph))
end associate
end function get_C66
!--------------------------------------------------------------------------------------------------
!> @brief returns 6x6 elasticity tensor in Voigt notation
!--------------------------------------------------------------------------------------------------
module function elastic_C66(ph,en) result(C66)
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(6,6) :: &
C66
associate(prm => param(ph))
C66 = get_C66(ph,en)
C66 = math_sym3333to66(math_Voigt66to3333(C66)) ! Literature data is in Voigt notation
end associate
end function elastic_C66
!--------------------------------------------------------------------------------------------------
!> @brief returns value of shear modulus
!--------------------------------------------------------------------------------------------------
module function elastic_mu(ph,en) result(mu)
integer, intent(in) :: &
ph, &
en
real(pReal) :: &
mu
mu = lattice_equivalent_mu(get_C66(ph,en),'voigt')
end function elastic_mu
!--------------------------------------------------------------------------------------------------
!> @brief returns value of poisson ratio
!--------------------------------------------------------------------------------------------------
module function elastic_nu(ph,en) result(nu)
integer, intent(in) :: &
ph, &
en
real(pReal) :: &
nu
nu = lattice_equivalent_nu(get_C66(ph,en),'voigt')
end function elastic_nu
!--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
@ -98,7 +176,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
do i =1, 3;do j=1,3
dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
enddo; enddo
end do; end do
end subroutine phase_hooke_SandItsTangents
@ -116,38 +194,10 @@ module function phase_homogenizedC(ph,en) result(C)
case (PLASTICITY_DISLOTWIN_ID) plasticType
C = plastic_dislotwin_homogenizedC(ph,en)
case default plasticType
C = param(ph)%C66
C = elastic_C66(ph,en)
end select plasticType
end function phase_homogenizedC
module function elastic_C66(ph) result(C66)
real(pReal), dimension(6,6) :: C66
integer, intent(in) :: ph
C66 = param(ph)%C66
end function elastic_C66
module function elastic_mu(ph) result(mu)
real(pReal) :: mu
integer, intent(in) :: ph
mu = param(ph)%mu
end function elastic_mu
module function elastic_nu(ph) result(nu)
real(pReal) :: nu
integer, intent(in) :: ph
nu = param(ph)%nu
end function elastic_nu
end submodule elastic

View File

@ -13,7 +13,6 @@ submodule(phase:plastic) dislotungsten
type :: tParameters
real(pReal) :: &
D = 1.0_pReal, & !< grain size
mu = 1.0_pReal, & !< equivalent shear modulus
D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient
Q_cl = 1.0_pReal !< activation energy for dislocation climb
real(pReal), allocatable, dimension(:) :: &
@ -130,8 +129,6 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
#endif
prm%mu = elastic_mu(ph)
!--------------------------------------------------------------------------------------------------
! slip related parameters
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
@ -324,10 +321,13 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
dot_rho_dip_formation, &
dot_rho_dip_climb, &
d_hat
real(pReal) :: &
mu
associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
mu = elastic_mu(ph,en)
call kinetics(Mp,T,ph,en,&
dot_gamma_pos,dot_gamma_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg)
@ -338,13 +338,13 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
dot_rho_dip_formation = 0.0_pReal
dot_rho_dip_climb = 0.0_pReal
else where
d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), &
d_hat = math_clip(3.0_pReal*mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), &
prm%d_caron, & ! lower limit
dst%Lambda_sl(:,en)) ! upper limit
dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, &
0.0_pReal, &
prm%dipoleformation)
v_cl = (3.0_pReal*prm%mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) &
v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) &
* (1.0_pReal/(d_hat+prm%d_caron))
dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency?
end where
@ -372,11 +372,12 @@ module subroutine dislotungsten_dependentState(ph,en)
real(pReal), dimension(param(ph)%sum_N_sl) :: &
Lambda_sl_inv
real(pReal) :: &
mu
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
dst%tau_pass(:,en) = prm%mu*prm%b_sl &
dst%tau_pass(:,en) = elastic_mu(ph,en)*prm%b_sl &
* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
Lambda_sl_inv = 1.0_pReal/prm%D &

View File

@ -14,8 +14,6 @@ submodule(phase:plastic) dislotwin
type :: tParameters
real(pReal) :: &
mu = 1.0_pReal, & !< equivalent shear modulus
nu = 1.0_pReal, & !< equivalent shear Poisson's ratio
Q_cl = 1.0_pReal, & !< activation energy for dislocation climb
omega = 1.0_pReal, & !< frequency factor for dislocation climb
D = 1.0_pReal, & !< grain size
@ -33,7 +31,9 @@ submodule(phase:plastic) dislotwin
delta_G = 1.0_pReal, & !< Free energy difference between austensite and martensite
i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation
h = 1.0_pReal, & !< Stack height of hex nucleus
T_ref = 0.0_pReal
T_ref = 0.0_pReal, &
a_cI = 1.0_pReal, &
a_cF = 1.0_pReal
real(pReal), dimension(2) :: &
Gamma_sf = 0.0_pReal
real(pReal), allocatable, dimension(:) :: &
@ -62,20 +62,22 @@ submodule(phase:plastic) dislotwin
h_sl_tr, & !< components of slip-trans interaction matrix
h_tr_tr, & !< components of trans-trans interaction matrix
n0_sl, & !< slip system normal
forestProjection, &
C66
forestProjection
real(pReal), allocatable, dimension(:,:,:) :: &
P_sl, &
P_tw, &
P_tr, &
C66_tw, &
C66_tr
P_tr
integer :: &
sum_N_sl, & !< total number of active slip system
sum_N_tw, & !< total number of active twin system
sum_N_tr !< total number of active transformation system
integer, allocatable, dimension(:) :: &
N_tw, &
N_tr
integer, allocatable, dimension(:,:) :: &
fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans
character(len=:), allocatable :: &
lattice_tr
character(len=pStringLen), allocatable, dimension(:) :: &
output
logical :: &
@ -134,7 +136,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
sizeState, sizeDotState, &
startIndex, endIndex
integer, dimension(:), allocatable :: &
N_sl, N_tw, N_tr
N_sl
real(pReal), allocatable, dimension(:) :: &
rho_mob_0, & !< initial unipolar dislocation density per slip system
rho_dip_0 !< initial dipole dislocation density per slip system
@ -185,10 +187,6 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
#endif
! This data is read in already in lattice
prm%mu = elastic_mu(ph)
prm%nu = elastic_nu(ph)
prm%C66 = elastic_C66(ph)
!--------------------------------------------------------------------------------------------------
! slip related parameters
@ -260,35 +258,33 @@ module function plastic_dislotwin_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! twin related parameters
N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(abs(N_tw))
prm%N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(abs(prm%N_tw))
twinActive: if (prm%sum_N_tw > 0) then
prm%systems_tw = lattice_labels_twin(N_tw,phase_lattice(ph))
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'), &
prm%systems_tw = lattice_labels_twin(prm%N_tw,phase_lattice(ph))
prm%P_tw = lattice_SchmidMatrix_twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,pl%get_as1dFloat('h_tw-tw'), &
phase_lattice(ph))
prm%b_tw = pl%get_as1dFloat('b_tw', requiredSize=size(N_tw))
prm%t_tw = pl%get_as1dFloat('t_tw', requiredSize=size(N_tw))
prm%r = pl%get_as1dFloat('p_tw', requiredSize=size(N_tw))
prm%b_tw = pl%get_as1dFloat('b_tw', requiredSize=size(prm%N_tw))
prm%t_tw = pl%get_as1dFloat('t_tw', requiredSize=size(prm%N_tw))
prm%r = pl%get_as1dFloat('p_tw', requiredSize=size(prm%N_tw))
prm%x_c_tw = pl%get_asFloat('x_c_tw')
prm%L_tw = pl%get_asFloat('L_tw')
prm%i_tw = pl%get_asFloat('i_tw')
prm%gamma_char= lattice_characteristicShear_Twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%C66_tw = lattice_C66_twin(N_tw,prm%C66,phase_lattice(ph),phase_cOverA(ph))
prm%gamma_char= lattice_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
if (.not. prm%fccTwinTransNucleation) then
prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw')
prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,N_tw)
prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,prm%N_tw)
endif
! expand: family => system
prm%b_tw = math_expand(prm%b_tw,N_tw)
prm%t_tw = math_expand(prm%t_tw,N_tw)
prm%r = math_expand(prm%r,N_tw)
prm%b_tw = math_expand(prm%b_tw,prm%N_tw)
prm%t_tw = math_expand(prm%t_tw,prm%N_tw)
prm%r = math_expand(prm%r,prm%N_tw)
! sanity checks
if ( prm%x_c_tw < 0.0_pReal) extmsg = trim(extmsg)//' x_c_tw'
@ -307,39 +303,38 @@ module function plastic_dislotwin_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! transformation related parameters
N_tr = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray)
prm%sum_N_tr = sum(abs(N_tr))
prm%N_tr = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray)
prm%sum_N_tr = sum(abs(prm%N_tr))
transActive: if (prm%sum_N_tr > 0) then
prm%b_tr = pl%get_as1dFloat('b_tr')
prm%b_tr = math_expand(prm%b_tr,N_tr)
prm%b_tr = math_expand(prm%b_tr,prm%N_tr)
prm%h = pl%get_asFloat('h', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%i_tr = pl%get_asFloat('i_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%delta_G = pl%get_asFloat('delta_G')
prm%x_c_tr = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%L_tr = pl%get_asFloat('L_tr')
prm%a_cI = pl%get_asFloat('a_cI', defaultVal=0.0_pReal)
prm%a_cF = pl%get_asFloat('a_cF', defaultVal=0.0_pReal)
prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_as1dFloat('h_tr-tr'),&
prm%lattice_tr = pl%get_asString('lattice_tr')
prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,pl%get_as1dFloat('h_tr-tr'),&
phase_lattice(ph))
prm%C66_tr = lattice_C66_trans(N_tr,prm%C66,pl%get_asString('lattice_tr'), &
0.0_pReal, &
pl%get_asFloat('a_cI', defaultVal=0.0_pReal), &
pl%get_asFloat('a_cF', defaultVal=0.0_pReal))
prm%P_tr = lattice_SchmidMatrix_trans(N_tr,pl%get_asString('lattice_tr'), &
prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr,prm%lattice_tr, &
0.0_pReal, &
pl%get_asFloat('a_cI', defaultVal=0.0_pReal), &
pl%get_asFloat('a_cF', defaultVal=0.0_pReal))
prm%a_cI, &
prm%a_cF)
if (phase_lattice(ph) /= 'cF') then
prm%dot_N_0_tr = pl%get_as1dFloat('dot_N_0_tr')
prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,N_tr)
prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr)
endif
prm%t_tr = pl%get_as1dFloat('t_tr')
prm%t_tr = math_expand(prm%t_tr,N_tr)
prm%t_tr = math_expand(prm%t_tr,prm%N_tr)
prm%s = pl%get_as1dFloat('p_tr',defaultVal=[0.0_pReal])
prm%s = math_expand(prm%s,N_tr)
prm%s = math_expand(prm%s,prm%N_tr)
! sanity checks
if ( prm%x_c_tr < 0.0_pReal) extmsg = trim(extmsg)//' x_c_tr'
@ -386,15 +381,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
endif
slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), &
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dFloat('h_sl-tw'), &
phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
if (prm%fccTwinTransNucleation .and. size(prm%N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
endif slipAndTwinActive
slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then
prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,pl%get_as1dFloat('h_sl-tr'), &
prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,prm%N_tr,pl%get_as1dFloat('h_sl-tr'), &
phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
if (prm%fccTwinTransNucleation .and. size(prm%N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
endif slipAndTransActive
!--------------------------------------------------------------------------------------------------
@ -478,27 +473,40 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
integer, intent(in) :: &
ph, en
real(pReal), dimension(6,6) :: &
homogenizedC
homogenizedC, &
C66
real(pReal), dimension(:,:,:), allocatable :: &
C66_tw, &
C66_tr
integer :: i
real(pReal) :: f_unrotated
associate(prm => param(ph), stt => state(ph))
C66 = elastic_C66(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))
homogenizedC = f_unrotated * prm%C66
do i=1,prm%sum_N_tw
homogenizedC = homogenizedC &
+ stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i)
enddo
do i=1,prm%sum_N_tr
homogenizedC = homogenizedC &
+ stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i)
enddo
homogenizedC = f_unrotated * C66
twinActive: if (prm%sum_N_tw > 0) then
C66_tw = lattice_C66_twin(prm%N_tw,C66,phase_lattice(ph),phase_cOverA(ph))
do i=1,prm%sum_N_tw
homogenizedC = homogenizedC &
+ stt%f_tw(i,en)*C66_tw(1:6,1:6,i)
end do
end if twinActive
transActive: if (prm%sum_N_tr > 0) then
C66_tr = lattice_C66_trans(prm%N_tr,C66,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF)
do i=1,prm%sum_N_tr
homogenizedC = homogenizedC &
+ stt%f_tr(i,en)*C66_tr(1:6,1:6,i)
end do
end if transActive
end associate
@ -647,10 +655,15 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
dot_gamma_tw
real(pReal), dimension(param(ph)%sum_N_tr) :: &
dot_gamma_tr
real(pReal) :: &
mu, &
nu
associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
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))
@ -665,7 +678,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
dot_rho_dip_formation(i) = 0.0_pReal
dot_rho_dip_climb(i) = 0.0_pReal
else significantSlipStress
d_hat = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
d_hat = 3.0_pReal*mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
d_hat = math_clip(d_hat, right = dst%Lambda_sl(i,en))
d_hat = math_clip(d_hat, left = prm%d_caron(i))
@ -677,8 +690,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
else
! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) &
* (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (prm%mu*prm%b_sl(i)), &
b_d = merge(24.0_pReal*PI*(1.0_pReal - nu)/(2.0_pReal + nu) &
* (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (mu*prm%b_sl(i)), &
1.0_pReal, &
prm%ExtendedDislocations)
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
@ -732,10 +745,15 @@ module subroutine dislotwin_dependentState(T,ph,en)
f_over_t_tr
real(pReal), dimension(:), allocatable :: &
x0
real(pReal) :: &
mu, &
nu
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
mu = elastic_mu(ph,en)
nu = elastic_nu(ph,en)
sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en))
sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
@ -760,24 +778,26 @@ module subroutine dislotwin_dependentState(T,ph,en)
dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
!* threshold stress for dislocation motion
dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
dst%tau_pass(:,en) = mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
!* threshold stress for growing twin/martensite
dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) &
+ 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_tw)
+ 3.0_pReal*prm%b_tw*mu/(prm%L_tw*prm%b_tw)
dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
+ 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_tr) &
+ 3.0_pReal*prm%b_tr*mu/(prm%L_tr*prm%b_tr) &
+ prm%h*prm%delta_G/(3.0_pReal*prm%b_tr)
dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal
dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal
x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
x0 = mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/&
(1.0_pReal-nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tw(:,en) = mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
x0 = mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/&
(1.0_pReal-nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tr(:,en) = mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
end associate

View File

@ -242,9 +242,6 @@ module function plastic_nonlocal_init() result(myPlasticity)
prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
prm%mu = elastic_mu(ph)
prm%nu = elastic_nu(ph)
ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(ini%N_sl))
slipActive: if (prm%sum_N_sl > 0) then
@ -570,7 +567,9 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
n
real(pReal) :: &
FVsize, &
nRealNeighbors ! number of really existing neighbors
nRealNeighbors, & ! number of really existing neighbors
mu, &
nu
integer, dimension(2) :: &
neighbors
real(pReal), dimension(2) :: &
@ -609,6 +608,8 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
associate(prm => param(ph),dst => dependentState(ph), stt => state(ph))
mu = elastic_mu(ph,en)
nu = elastic_nu(ph,en)
rho = getRho(ph,en)
stt%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
@ -627,7 +628,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
myInteractionMatrix = prm%h_sl_sl
endif
dst%tau_pass(:,en) = prm%mu * prm%b_sl &
dst%tau_pass(:,en) = mu * prm%b_sl &
* sqrt(matmul(myInteractionMatrix,sum(abs(rho),2)))
!*** calculate the dislocation stress of the neighboring excess dislocation densities
@ -728,8 +729,8 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal
! ... gives the local stress correction when multiplied with a factor
dst%tau_back(s,en) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) &
* ( rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) &
dst%tau_back(s,en) = - mu * prm%b_sl(s) / (2.0_pReal * PI) &
* ( rhoExcessGradient_over_rho(1) / (1.0_pReal - nu) &
+ rhoExcessGradient_over_rho(2))
enddo
endif
@ -855,6 +856,9 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
c, & ! character of dislocation
t, & ! type of dislocation
s ! index of my current slip system
real(pReal) :: &
mu, &
nu
real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
deltaRhoRemobilization, & ! density increment by remobilization
deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change)
@ -869,9 +873,12 @@ 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)
nu = elastic_nu(ph,en)
!*** shortcut to state variables
forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
forall (s = 1:prm%sum_N_sl, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en)
@ -901,8 +908,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo
dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
dUpper(:,1) = mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - nu) * abs(tau))
dUpper(:,2) = mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) &
dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1))
@ -975,7 +982,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
dLower, & !< minimum stable dipole distance for edges and screws
dUpper !< current maximum stable dipole distance for edges and screws
real(pReal) :: &
D_SD
D_SD, &
mu, &
nu
if (timestep <= 0.0_pReal) then
plasticState(ph)%dotState = 0.0_pReal
@ -984,6 +993,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph))
mu = elastic_mu(ph,en)
nu = elastic_nu(ph,en)
tau = 0.0_pReal
dot_gamma = 0.0_pReal
@ -1005,8 +1017,8 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
enddo
dLower = prm%minDipoleHeight
dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
dUpper(:,1) = mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - nu) * abs(tau))
dUpper(:,2) = mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) &
dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1))
@ -1083,8 +1095,8 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
! thermally activated annihilation of edge dipoles by climb
rhoDotThermalAnnihilation = 0.0_pReal
D_SD = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) ! eq. 3.53
v_climb = D_SD * prm%mu * prm%V_at &
/ (PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1)) * kB * Temperature) ! eq. 3.54
v_climb = D_SD * mu * prm%V_at &
/ (PI * (1.0_pReal-nu) * (dUpper(:,1) + dLower(:,1)) * kB * Temperature) ! eq. 3.54
forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) &
rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &