diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index a74361198..db2faf914 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -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 diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 8b05adc31..9c9a159a0 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -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 = 0.0_pReal, & + 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 type(tParameters), allocatable, dimension(:) :: param contains - +!-------------------------------------------------------------------------------------------------- +!> @brief Initialize elasticity. +!-------------------------------------------------------------------------------------------------- module subroutine elastic_init(phases) class(tNode), pointer :: & @@ -41,22 +45,15 @@ 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') + prm%C_13 = elastic%get_asFloat('C_13') + prm%C_33 = elastic%get_asFloat('C_33') end if - 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 + if (phase_lattice(ph) == 'tI') prm%C_66 = elastic%get_asFloat('C_66') end associate end do @@ -64,6 +61,93 @@ module subroutine elastic_init(phases) 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 !> 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 integer, intent(in) :: ph, en + plasticType: select case (phase_plasticity(ph)) 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 diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 33bec98cc..102e009fe 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -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 @@ -376,7 +376,7 @@ module subroutine dislotungsten_dependentState(ph,en) 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 & diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index d7ddff80a..62119c819 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -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) - end if + 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) - end if + 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,16 +381,16 @@ module function plastic_dislotwin_init() result(myPlasticity) end if 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' - end if slipAndTwinActive + 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' - end if slipAndTransActive + if (prm%fccTwinTransNucleation .and. size(prm%N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation' + endif slipAndTransActive !-------------------------------------------------------------------------------------------------- ! allocate state arrays @@ -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) - end do - do i=1,prm%sum_N_tr - homogenizedC = homogenizedC & - + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i) - end do + 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,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) !* 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 diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 727278e18..f53a16042 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -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 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))) !*** 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)) end do end if @@ -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 end do - 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, & end do 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) &