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

Merge branch 'dynamic-C' into 'development' 

preparing for temperature dependent C by calling it dynamically

See merge request damask/DAMASK!457
This commit is contained in:
Martin Diehl 2021-11-18 17:58:54 +00:00
parent 52d0ba4a74 eee455c0d3
commit 202af32429
5 changed files with 236 additions and 149 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
src/phase_mechanical.f90
View File

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

151
src/phase_mechanical_elastic.f90
View File

@ -1,18 +1,22 @@
submodule(phase:mechanical) elastic submodule(phase:mechanical) elastic
type :: tParameters type :: tParameters
real(pReal), dimension(6,6) :: &
C66 = 0.0_pReal !< Elastic constants in Voigt notation
real(pReal) :: & real(pReal) :: &
mu, & C_11 = 0.0_pReal, &
nu C_12 = 0.0_pReal, &
C_13 = 0.0_pReal, &
C_33 = 0.0_pReal, &
C_44 = 0.0_pReal, &
C_66 = 0.0_pReal
end type tParameters end type tParameters
type(tParameters), allocatable, dimension(:) :: param type(tParameters), allocatable, dimension(:) :: param
contains contains
!--------------------------------------------------------------------------------------------------
!> @brief Initialize elasticity.
!--------------------------------------------------------------------------------------------------
module subroutine elastic_init(phases) module subroutine elastic_init(phases)
class(tNode), pointer :: & class(tNode), pointer :: &
@ -41,22 +45,15 @@ module subroutine elastic_init(phases)
associate(prm => param(ph)) associate(prm => param(ph))
prm%C66(1,1) = elastic%get_asFloat('C_11') prm%C_11 = elastic%get_asFloat('C_11')
prm%C66(1,2) = elastic%get_asFloat('C_12') prm%C_12 = elastic%get_asFloat('C_12')
prm%C66(4,4) = elastic%get_asFloat('C_44') prm%C_44 = elastic%get_asFloat('C_44')
if (any(phase_lattice(ph) == ['hP','tI'])) then if (any(phase_lattice(ph) == ['hP','tI'])) then
prm%C66(1,3) = elastic%get_asFloat('C_13') prm%C_13 = elastic%get_asFloat('C_13')
prm%C66(3,3) = elastic%get_asFloat('C_33') prm%C_33 = elastic%get_asFloat('C_33')
end if end if
if (phase_lattice(ph) == 'tI') prm%C66(6,6) = elastic%get_asFloat('C_66') if (phase_lattice(ph) == 'tI') prm%C_66 = 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
end associate end associate
end do end do
@ -64,6 +61,93 @@ module subroutine elastic_init(phases)
end subroutine elastic_init end subroutine elastic_init
!--------------------------------------------------------------------------------------------------
!> @brief Return 6x6 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 Return 6x6 elasticity tensor.
!--------------------------------------------------------------------------------------------------
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 Return 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 Return 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 !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermediate deformation gradients using Hooke's law !> the elastic and intermediate deformation gradients using Hooke's law
@ -112,42 +196,15 @@ module function phase_homogenizedC(ph,en) result(C)
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
plasticType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_DISLOTWIN_ID) plasticType case (PLASTICITY_DISLOTWIN_ID) plasticType
C = plastic_dislotwin_homogenizedC(ph,en) C = plastic_dislotwin_homogenizedC(ph,en)
case default plasticType case default plasticType
C = param(ph)%C66 C = elastic_C66(ph,en)
end select plasticType end select plasticType
end function phase_homogenizedC 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 end submodule elastic

14
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -13,7 +13,6 @@ submodule(phase:plastic) dislotungsten
type :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: &
D = 1.0_pReal, & !< grain size D = 1.0_pReal, & !< grain size
mu = 1.0_pReal, & !< equivalent shear modulus
D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient
Q_cl = 1.0_pReal !< activation energy for dislocation climb Q_cl = 1.0_pReal !< activation energy for dislocation climb
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
@ -130,8 +129,6 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray) prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
#endif #endif
prm%mu = elastic_mu(ph)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) 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_formation, &
dot_rho_dip_climb, & dot_rho_dip_climb, &
d_hat d_hat
real(pReal) :: &
mu
associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
mu = elastic_mu(ph,en)
call kinetics(Mp,T,ph,en,& call kinetics(Mp,T,ph,en,&
dot_gamma_pos,dot_gamma_neg, & dot_gamma_pos,dot_gamma_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_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_formation = 0.0_pReal
dot_rho_dip_climb = 0.0_pReal dot_rho_dip_climb = 0.0_pReal
else where 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 prm%d_caron, & ! lower limit
dst%Lambda_sl(:,en)) ! upper 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, & 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, & 0.0_pReal, &
prm%dipoleformation) 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)) * (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? 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 end where
@ -376,7 +376,7 @@ module subroutine dislotungsten_dependentState(ph,en)
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) 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))) * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
Lambda_sl_inv = 1.0_pReal/prm%D & Lambda_sl_inv = 1.0_pReal/prm%D &

166
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -14,8 +14,6 @@ submodule(phase:plastic) dislotwin
type :: tParameters type :: tParameters
real(pReal) :: & 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 Q_cl = 1.0_pReal, & !< activation energy for dislocation climb
omega = 1.0_pReal, & !< frequency factor for dislocation climb omega = 1.0_pReal, & !< frequency factor for dislocation climb
D = 1.0_pReal, & !< grain size 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 delta_G = 1.0_pReal, & !< Free energy difference between austensite and martensite
i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation
h = 1.0_pReal, & !< Stack height of hex nucleus 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) :: & real(pReal), dimension(2) :: &
Gamma_sf = 0.0_pReal Gamma_sf = 0.0_pReal
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
@ -62,20 +62,22 @@ submodule(phase:plastic) dislotwin
h_sl_tr, & !< components of slip-trans interaction matrix h_sl_tr, & !< components of slip-trans interaction matrix
h_tr_tr, & !< components of trans-trans interaction matrix h_tr_tr, & !< components of trans-trans interaction matrix
n0_sl, & !< slip system normal n0_sl, & !< slip system normal
forestProjection, & forestProjection
C66
real(pReal), allocatable, dimension(:,:,:) :: & real(pReal), allocatable, dimension(:,:,:) :: &
P_sl, & P_sl, &
P_tw, & P_tw, &
P_tr, & P_tr
C66_tw, &
C66_tr
integer :: & integer :: &
sum_N_sl, & !< total number of active slip system sum_N_sl, & !< total number of active slip system
sum_N_tw, & !< total number of active twin system sum_N_tw, & !< total number of active twin system
sum_N_tr !< total number of active transformation system sum_N_tr !< total number of active transformation system
integer, allocatable, dimension(:) :: &
N_tw, &
N_tr
integer, allocatable, dimension(:,:) :: & integer, allocatable, dimension(:,:) :: &
fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans
character(len=:), allocatable :: &
lattice_tr
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
logical :: & logical :: &
@ -134,7 +136,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
sizeState, sizeDotState, & sizeState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
N_sl, N_tw, N_tr N_sl
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
rho_mob_0, & !< initial unipolar dislocation density per slip system rho_mob_0, & !< initial unipolar dislocation density per slip system
rho_dip_0 !< initial dipole 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) prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
#endif #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 ! slip related parameters
@ -260,35 +258,33 @@ module function plastic_dislotwin_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! twin related parameters ! twin related parameters
N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray) prm%N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(abs(N_tw)) prm%sum_N_tw = sum(abs(prm%N_tw))
twinActive: if (prm%sum_N_tw > 0) then twinActive: if (prm%sum_N_tw > 0) then
prm%systems_tw = lattice_labels_twin(N_tw,phase_lattice(ph)) prm%systems_tw = lattice_labels_twin(prm%N_tw,phase_lattice(ph))
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%P_tw = lattice_SchmidMatrix_twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'), & prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,pl%get_as1dFloat('h_tw-tw'), &
phase_lattice(ph)) phase_lattice(ph))
prm%b_tw = pl%get_as1dFloat('b_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(N_tw)) prm%t_tw = pl%get_as1dFloat('t_tw', requiredSize=size(prm%N_tw))
prm%r = pl%get_as1dFloat('p_tw', requiredSize=size(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%x_c_tw = pl%get_asFloat('x_c_tw')
prm%L_tw = pl%get_asFloat('L_tw') prm%L_tw = pl%get_asFloat('L_tw')
prm%i_tw = pl%get_asFloat('i_tw') prm%i_tw = pl%get_asFloat('i_tw')
prm%gamma_char= lattice_characteristicShear_Twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%gamma_char= lattice_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%C66_tw = lattice_C66_twin(N_tw,prm%C66,phase_lattice(ph),phase_cOverA(ph))
if (.not. prm%fccTwinTransNucleation) then if (.not. prm%fccTwinTransNucleation) then
prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw') 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)
end if endif
! expand: family => system ! expand: family => system
prm%b_tw = math_expand(prm%b_tw,N_tw) prm%b_tw = math_expand(prm%b_tw,prm%N_tw)
prm%t_tw = math_expand(prm%t_tw,N_tw) prm%t_tw = math_expand(prm%t_tw,prm%N_tw)
prm%r = math_expand(prm%r,N_tw) prm%r = math_expand(prm%r,prm%N_tw)
! sanity checks ! sanity checks
if ( prm%x_c_tw < 0.0_pReal) extmsg = trim(extmsg)//' x_c_tw' 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 ! transformation related parameters
N_tr = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray) prm%N_tr = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray)
prm%sum_N_tr = sum(abs(N_tr)) prm%sum_N_tr = sum(abs(prm%N_tr))
transActive: if (prm%sum_N_tr > 0) then transActive: if (prm%sum_N_tr > 0) then
prm%b_tr = pl%get_as1dFloat('b_tr') 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%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%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%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%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%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)) phase_lattice(ph))
prm%C66_tr = lattice_C66_trans(N_tr,prm%C66,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%P_tr = lattice_SchmidMatrix_trans(N_tr,pl%get_asString('lattice_tr'), &
0.0_pReal, & 0.0_pReal, &
pl%get_asFloat('a_cI', defaultVal=0.0_pReal), & prm%a_cI, &
pl%get_asFloat('a_cF', defaultVal=0.0_pReal)) prm%a_cF)
if (phase_lattice(ph) /= 'cF') then if (phase_lattice(ph) /= 'cF') then
prm%dot_N_0_tr = pl%get_as1dFloat('dot_N_0_tr') 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)
end if endif
prm%t_tr = pl%get_as1dFloat('t_tr') 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 = 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 ! sanity checks
if ( prm%x_c_tr < 0.0_pReal) extmsg = trim(extmsg)//' x_c_tr' if ( prm%x_c_tr < 0.0_pReal) extmsg = trim(extmsg)//' x_c_tr'
@ -386,16 +381,16 @@ module function plastic_dislotwin_init() result(myPlasticity)
end if end if
slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then 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)) 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'
end if slipAndTwinActive endif slipAndTwinActive
slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then 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)) 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'
end if slipAndTransActive endif slipAndTransActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
@ -478,27 +473,40 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
integer, intent(in) :: & integer, intent(in) :: &
ph, en ph, en
real(pReal), dimension(6,6) :: & real(pReal), dimension(6,6) :: &
homogenizedC homogenizedC, &
C66
real(pReal), dimension(:,:,:), allocatable :: &
C66_tw, &
C66_tr
integer :: i integer :: i
real(pReal) :: f_unrotated real(pReal) :: f_unrotated
associate(prm => param(ph), stt => state(ph)) associate(prm => param(ph), stt => state(ph))
C66 = elastic_C66(ph,en)
f_unrotated = 1.0_pReal & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,en)) & - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,en)) - sum(stt%f_tr(1:prm%sum_N_tr,en))
homogenizedC = f_unrotated * prm%C66 homogenizedC = f_unrotated * C66
do i=1,prm%sum_N_tw
homogenizedC = homogenizedC & twinActive: if (prm%sum_N_tw > 0) then
+ stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i) C66_tw = lattice_C66_twin(prm%N_tw,C66,phase_lattice(ph),phase_cOverA(ph))
end do do i=1,prm%sum_N_tw
do i=1,prm%sum_N_tr homogenizedC = homogenizedC &
homogenizedC = homogenizedC & + stt%f_tw(i,en)*C66_tw(1:6,1:6,i)
+ stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i) end do
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 end associate
@ -647,10 +655,15 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
dot_gamma_tw dot_gamma_tw
real(pReal), dimension(param(ph)%sum_N_tr) :: & real(pReal), dimension(param(ph)%sum_N_tr) :: &
dot_gamma_tr dot_gamma_tr
real(pReal) :: &
mu, &
nu
associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) 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 & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,en)) & - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,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_formation(i) = 0.0_pReal
dot_rho_dip_climb(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal
else significantSlipStress 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, right = dst%Lambda_sl(i,en))
d_hat = math_clip(d_hat, left = prm%d_caron(i)) d_hat = math_clip(d_hat, left = prm%d_caron(i))
@ -677,8 +690,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
else else
! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 ! 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))) 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) & b_d = merge(24.0_pReal*PI*(1.0_pReal - nu)/(2.0_pReal + nu) &
* (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (prm%mu*prm%b_sl(i)), & * (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (mu*prm%b_sl(i)), &
1.0_pReal, & 1.0_pReal, &
prm%ExtendedDislocations) prm%ExtendedDislocations)
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & 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 f_over_t_tr
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
x0 x0
real(pReal) :: &
mu, &
nu
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) 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_tw = sum(stt%f_tw(1:prm%sum_N_tw,en))
sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en)) sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
@ -760,24 +778,24 @@ 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) dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
!* threshold stress for dislocation motion !* 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 !* threshold stress for growing twin/martensite
dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) & 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) & 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) + 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_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 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 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) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) 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 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) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0) 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 end associate

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