From 5267794ff20d40f4dbfec2226f7293a903bb8f73 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 26 Jan 2022 07:24:29 +0100 Subject: [PATCH 01/17] second order expansion of temperature-dependent variables --- src/phase_mechanical_plastic_dislotwin.f90 | 59 ++++++++++++++-------- 1 file changed, 37 insertions(+), 22 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index bc6fb0ee6..f200f6091 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -25,14 +25,14 @@ submodule(phase:plastic) dislotwin xi_sb = 1.0_pReal, & !< value for shearband resistance v_sb = 1.0_pReal, & !< value for shearband velocity_0 E_sb = 1.0_pReal, & !< activation energy for shear bands - 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 = T_ROOM, & a_cI = 1.0_pReal, & a_cF = 1.0_pReal - real(pReal), dimension(2) :: & - Gamma_sf = 0.0_pReal + real(pReal), dimension(3) :: & + Gamma_sf = 0.0_pReal, & !< stacking fault energy + Delta_G = 0.0_pReal !< free energy difference between austensite and martensite real(pReal), allocatable, dimension(:) :: & b_sl, & !< absolute length of Burgers vector [m] for each slip system b_tw, & !< absolute length of Burgers vector [m] for each twin system @@ -303,7 +303,9 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%h = pl%get_asFloat('h', defaultVal=0.0_pReal) ! ToDo: This is not optional! prm%i_tr = pl%get_asFloat('i_tr', defaultVal=0.0_pReal) ! ToDo: This is not optional! - prm%delta_G = pl%get_asFloat('delta_G') + prm%Delta_G(1) = pl%get_asFloat('Delta_G') + prm%Delta_G(2) = pl%get_asFloat('Delta_G,T', defaultVal=0.0_pReal) + prm%Delta_G(3) = pl%get_asFloat('Delta_G,T^2',defaultVal=0.0_pReal) prm%x_c_tr = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: This is not optional! prm%L_tr = pl%get_asFloat('L_tr') prm%a_cI = pl%get_asFloat('a_cI', defaultVal=0.0_pReal) @@ -363,6 +365,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%T_ref = pl%get_asFloat('T_ref') prm%Gamma_sf(1) = pl%get_asFloat('Gamma_sf') prm%Gamma_sf(2) = pl%get_asFloat('Gamma_sf,T',defaultVal=0.0_pReal) + prm%Gamma_sf(3) = pl%get_asFloat('Gamma_sf,T^2',defaultVal=0.0_pReal) end if slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then @@ -559,7 +562,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) end do slipContribution - call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) + if (prm%sum_N_tw > 0) 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) & @@ -567,7 +570,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) end do twinContibution - call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) + if (prm%sum_N_tr > 0) 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) & @@ -641,13 +644,12 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) real(pReal) :: & mu, & nu, & - Gamma + Gamma_sf associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) - Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref) f_matrix = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,en)) & @@ -674,8 +676,11 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) dot_rho_dip_climb(i) = 0.0_pReal else ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 + Gamma_sf = prm%Gamma_sf(1) & + + prm%Gamma_sf(2) * (T-prm%T_ref) & + + prm%Gamma_sf(3) * (T-prm%T_ref)**2 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 - nu)/(2.0_pReal + nu) * Gamma / (mu*prm%b_sl(i)), & + b_d = merge(24.0_pReal*PI*(1.0_pReal - nu)/(2.0_pReal + nu) * Gamma_sf / (mu*prm%b_sl(i)), & 1.0_pReal, & prm%ExtendedDislocations) v_cl = 2.0_pReal*prm%omega*b_d**2*exp(-prm%Q_cl/(K_B*T)) & @@ -695,10 +700,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & - dot_rho_dip_climb - call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) + if (prm%sum_N_tw > 0) call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) 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) + if (prm%sum_N_tr > 0) call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) dot%f_tr(:,en) = f_matrix*dot_gamma_tr end associate @@ -718,7 +723,8 @@ module subroutine dislotwin_dependentState(T,ph,en) T real(pReal) :: & - sumf_tw,Gamma,sumf_tr + Gamma_sf, Delta_G, & + sumf_tw, sumf_tr real(pReal), dimension(param(ph)%sum_N_sl) :: & inv_lambda_sl real(pReal), dimension(param(ph)%sum_N_tw) :: & @@ -739,7 +745,12 @@ module subroutine dislotwin_dependentState(T,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)) - Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref) + Gamma_sf = prm%Gamma_sf(1) & + + prm%Gamma_sf(2) * (T-prm%T_ref) & + + prm%Gamma_sf(3) * (T-prm%T_ref)**2 + Delta_G = prm%Delta_G(1) & + + prm%Delta_G(2) * (T-prm%T_ref) & + + prm%Delta_G(3) * (T-prm%T_ref)**2 !* rescaled volume fraction for topology f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw ! this is per system ... @@ -763,11 +774,11 @@ module subroutine dislotwin_dependentState(T,ph,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) & + dst%tau_hat_tw(:,en) = Gamma_sf/(3.0_pReal*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_sf/(3.0_pReal*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*Delta_G/(3.0_pReal*prm%b_tr) end associate @@ -936,7 +947,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& tau, tau_r, & dot_N_0, & x0, V, & - Gamma, & + Gamma_sf, & mu, nu, & P_ncs, dP_ncs_dtau, & P, dP_dtau @@ -950,11 +961,13 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& isFCC: if (prm%fccTwinTransNucleation) then mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) - Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref) + Gamma_sf = prm%Gamma_sf(1) & + + prm%Gamma_sf(2) * (T-prm%T_ref) & + + prm%Gamma_sf(3) * (T-prm%T_ref)**2 do i = 1, prm%sum_N_tw tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) - x0 = mu*prm%b_tw(i)**2*(2.0_pReal+nu)/(Gamma*8.0_pReal*PI*(1.0_pReal-nu)) ! ToDo: In the paper, the Burgers vector for slip is used + x0 = mu*prm%b_tw(i)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu)) ! ToDo: In the paper, the Burgers vector for slip is used tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(PI/3.0_pReal)/x0) ! ToDo: In the paper, the Burgers vector for slip is used if (tau > tol_math_check .and. tau < tau_r) then @@ -1026,7 +1039,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& tau, tau_r, & dot_N_0, & x0, V, & - Gamma, & + Gamma_sf, & mu, nu, & P_ncs, dP_ncs_dtau, & P, dP_dtau @@ -1039,11 +1052,13 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) - Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref) + Gamma_sf = prm%Gamma_sf(1) & + + prm%Gamma_sf(2) * (T-prm%T_ref) & + + prm%Gamma_sf(3) * (T-prm%T_ref)**2 do i = 1, prm%sum_N_tr tau = math_tensordot(Mp,prm%P_tr(1:3,1:3,i)) - x0 = mu*prm%b_tr(i)**2*(2.0_pReal+nu)/(Gamma*8.0_pReal*PI*(1.0_pReal-nu)) ! ToDo: In the paper, the Burgers vector for slip is used + x0 = mu*prm%b_tr(i)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu)) ! ToDo: In the paper, the Burgers vector for slip is used tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(PI/3.0_pReal)/x0) ! ToDo: In the paper, the Burgers vector for slip is used if (tau > tol_math_check .and. tau < tau_r) then From 60e6e90874f720b56f307647ff861d28fb063188 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 26 Jan 2022 22:21:34 +0100 Subject: [PATCH 02/17] avoid global states --- src/phase_mechanical_plastic_dislotwin.f90 | 53 ++++++++-------------- 1 file changed, 18 insertions(+), 35 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index f200f6091..485b5522a 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -99,9 +99,7 @@ submodule(phase:plastic) dislotwin Lambda_sl, & !< mean free path between 2 obstacles seen by a moving dislocation Lambda_tw, & !< mean free path between 2 obstacles seen by a growing twin Lambda_tr, & !< mean free path between 2 obstacles seen by a growing martensite - tau_pass, & !< threshold stress for slip - tau_hat_tw, & !< threshold stress for twinning - tau_hat_tr !< threshold stress for transformation + tau_pass !< threshold stress for slip end type tDislotwinDependentState !-------------------------------------------------------------------------------------------------- @@ -429,14 +427,10 @@ module function plastic_dislotwin_init() result(myPlasticity) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tr',defaultVal=1.0e-6_pReal) if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tr' - allocate(dst%Lambda_sl (prm%sum_N_sl,Nmembers),source=0.0_pReal) allocate(dst%tau_pass (prm%sum_N_sl,Nmembers),source=0.0_pReal) - + allocate(dst%Lambda_sl (prm%sum_N_sl,Nmembers),source=0.0_pReal) allocate(dst%Lambda_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal) - allocate(dst%tau_hat_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal) - allocate(dst%Lambda_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal) - allocate(dst%tau_hat_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal) end associate @@ -723,7 +717,6 @@ module subroutine dislotwin_dependentState(T,ph,en) T real(pReal) :: & - Gamma_sf, Delta_G, & sumf_tw, sumf_tr real(pReal), dimension(param(ph)%sum_N_sl) :: & inv_lambda_sl @@ -740,18 +733,9 @@ module subroutine dislotwin_dependentState(T,ph,en) 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)) - Gamma_sf = prm%Gamma_sf(1) & - + prm%Gamma_sf(2) * (T-prm%T_ref) & - + prm%Gamma_sf(3) * (T-prm%T_ref)**2 - Delta_G = prm%Delta_G(1) & - + prm%Delta_G(2) * (T-prm%T_ref) & - + prm%Delta_G(3) * (T-prm%T_ref)**2 - !* rescaled volume fraction for topology f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw ! this is per system ... f_over_t_tr = sumf_tr/prm%t_tr ! but this not @@ -773,13 +757,6 @@ module subroutine dislotwin_dependentState(T,ph,en) !* threshold stress for dislocation motion 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_sf/(3.0_pReal*prm%b_tw) & - + 3.0_pReal*prm%b_tw*mu/(prm%L_tw*prm%b_tw) - dst%tau_hat_tr(:,en) = Gamma_sf/(3.0_pReal*prm%b_tr) & - + 3.0_pReal*prm%b_tr*mu/(prm%L_tr*prm%b_tr) & - + prm%h*Delta_G/(3.0_pReal*prm%b_tr) - end associate end subroutine dislotwin_dependentState @@ -824,9 +801,6 @@ module subroutine plastic_dislotwin_results(ph,group) case('Lambda_tw') call results_writeDataset(dst%Lambda_tw,group,trim(prm%output(ou)), & 'mean free path for twinning','m',prm%systems_tw) - case('tau_hat_tw') - call results_writeDataset(dst%tau_hat_tw,group,trim(prm%output(ou)), & - 'threshold stress for twinning','Pa',prm%systems_tw) case('f_tr') if (prm%sum_N_tr>0) call results_writeDataset(stt%f_tr,group,trim(prm%output(ou)), & @@ -944,7 +918,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& ddot_gamma_dtau_tw real :: & - tau, tau_r, & + tau, tau_r, tau_hat, & dot_N_0, & x0, V, & Gamma_sf, & @@ -959,11 +933,14 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) isFCC: if (prm%fccTwinTransNucleation) then + mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) Gamma_sf = prm%Gamma_sf(1) & + prm%Gamma_sf(2) * (T-prm%T_ref) & + prm%Gamma_sf(3) * (T-prm%T_ref)**2 + tau_hat = 3.0_pReal*prm%b_tw(1)*mu/prm%L_tw & + + Gamma_sf/(3.0_pReal*prm%b_tw(1)) do i = 1, prm%sum_N_tw tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) @@ -971,8 +948,8 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(PI/3.0_pReal)/x0) ! ToDo: In the paper, the Burgers vector for slip is used if (tau > tol_math_check .and. tau < tau_r) then - P = exp(-(dst%tau_hat_tw(i,en)/tau)**prm%r(i)) - dP_dTau = prm%r(i) * (dst%tau_hat_tw(i,en)/tau)**prm%r(i)/tau * P + P = exp(-(tau_hat/tau)**prm%r(i)) + dP_dTau = prm%r(i) * (tau_hat/tau)**prm%r(i)/tau * P s = prm%fcc_twinNucleationSlipPair(1:2,i) dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) & @@ -1036,10 +1013,10 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& ddot_gamma_dtau_tr real :: & - tau, tau_r, & + tau, tau_r, tau_hat, & dot_N_0, & x0, V, & - Gamma_sf, & + Gamma_sf, Delta_G, & mu, nu, & P_ncs, dP_ncs_dtau, & P, dP_dtau @@ -1055,6 +1032,12 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& Gamma_sf = prm%Gamma_sf(1) & + prm%Gamma_sf(2) * (T-prm%T_ref) & + prm%Gamma_sf(3) * (T-prm%T_ref)**2 + Delta_G = prm%Delta_G(1) & + + prm%Delta_G(2) * (T-prm%T_ref) & + + prm%Delta_G(3) * (T-prm%T_ref)**2 + tau_hat = 3.0_pReal*prm%b_tr(1)*mu/prm%L_tr & + + Gamma_sf/(3.0_pReal*prm%b_tr(1)) & + + prm%h*Delta_G/(3.0_pReal*prm%b_tr(1)) do i = 1, prm%sum_N_tr tau = math_tensordot(Mp,prm%P_tr(1:3,1:3,i)) @@ -1062,8 +1045,8 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(PI/3.0_pReal)/x0) ! ToDo: In the paper, the Burgers vector for slip is used if (tau > tol_math_check .and. tau < tau_r) then - P = exp(-(dst%tau_hat_tr(i,en)/tau)**prm%s(i)) - dP_dTau = prm%s(i) * (dst%tau_hat_tr(i,en)/tau)**prm%s(i)/tau * P + P = exp(-(tau_hat/tau)**prm%s(i)) + dP_dTau = prm%s(i) * (tau_hat/tau)**prm%s(i)/tau * P s = prm%fcc_twinNucleationSlipPair(1:2,i) dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) & From d9f96bc0ecba829032ace1ccaf12817d94453528 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 26 Jan 2022 22:28:24 +0100 Subject: [PATCH 03/17] following paper, eq. (36) and eq. (37) --- src/phase_mechanical_plastic_dislotwin.f90 | 27 +++++++++++----------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 485b5522a..36d2899b9 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -17,15 +17,14 @@ submodule(phase:plastic) dislotwin p_sb = 1.0_pReal, & !< p-exponent in shear band velocity q_sb = 1.0_pReal, & !< q-exponent in shear band velocity i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning - L_tw = 1.0_pReal, & !< length of twin nuclei in Burgers vectors: TODO unit should be meters - L_tr = 1.0_pReal, & !< length of trans nuclei in Burgers vectors: TODO unit should be meters - x_c_tw = 1.0_pReal, & !< critical distance for formation of twin nucleus - x_c_tr = 1.0_pReal, & !< critical distance for formation of trans nucleus + i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation + L_tw = 1.0_pReal, & !< length of twin nuclei + L_tr = 1.0_pReal, & !< length of trans nuclei + x_c = 1.0_pReal, & !< critical distance for formation of twin/trans nucleus V_cs = 1.0_pReal, & !< cross slip volume xi_sb = 1.0_pReal, & !< value for shearband resistance v_sb = 1.0_pReal, & !< value for shearband velocity_0 E_sb = 1.0_pReal, & !< activation energy for shear bands - i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation h = 1.0_pReal, & !< stack height of hex nucleus T_ref = T_ROOM, & a_cI = 1.0_pReal, & @@ -260,7 +259,6 @@ module function plastic_dislotwin_init() result(myPlasticity) 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') @@ -277,7 +275,6 @@ module function plastic_dislotwin_init() result(myPlasticity) 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' if ( prm%L_tw < 0.0_pReal) extmsg = trim(extmsg)//' L_tw' if ( prm%i_tw < 0.0_pReal) extmsg = trim(extmsg)//' i_tw' if (any(prm%b_tw < 0.0_pReal)) extmsg = trim(extmsg)//' b_tw' @@ -304,7 +301,6 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%Delta_G(1) = pl%get_asFloat('Delta_G') prm%Delta_G(2) = pl%get_asFloat('Delta_G,T', defaultVal=0.0_pReal) prm%Delta_G(3) = pl%get_asFloat('Delta_G,T^2',defaultVal=0.0_pReal) - prm%x_c_tr = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: This is not optional! 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) @@ -325,7 +321,6 @@ module function plastic_dislotwin_init() result(myPlasticity) 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' if ( prm%L_tr < 0.0_pReal) extmsg = trim(extmsg)//' L_tr' if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr' if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr' @@ -356,8 +351,12 @@ module function plastic_dislotwin_init() result(myPlasticity) if (prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) & prm%D = pl%get_asFloat('D') - if (prm%sum_N_tw + prm%sum_N_tr > 0) & + if (prm%sum_N_tw + prm%sum_N_tr > 0) then + prm%x_c = pl%get_asFloat('x_c') prm%V_cs = pl%get_asFloat('V_cs') + if (prm%x_c < 0.0_pReal) extmsg = trim(extmsg)//' x_c' + if (prm%V_cs < 0.0_pReal) extmsg = trim(extmsg)//' V_cs' + end if if (prm%sum_N_tw + prm%sum_N_tr > 0 .or. prm%ExtendedDislocations) then prm%T_ref = pl%get_asFloat('T_ref') @@ -941,11 +940,11 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& + prm%Gamma_sf(3) * (T-prm%T_ref)**2 tau_hat = 3.0_pReal*prm%b_tw(1)*mu/prm%L_tw & + Gamma_sf/(3.0_pReal*prm%b_tw(1)) + x0 = mu*prm%b_sl(1)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu)) + tau_r = mu*prm%b_sl(1)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c)+cos(PI/3.0_pReal)/x0) do i = 1, prm%sum_N_tw tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) - x0 = mu*prm%b_tw(i)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu)) ! ToDo: In the paper, the Burgers vector for slip is used - tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(PI/3.0_pReal)/x0) ! ToDo: In the paper, the Burgers vector for slip is used if (tau > tol_math_check .and. tau < tau_r) then P = exp(-(tau_hat/tau)**prm%r(i)) @@ -1038,11 +1037,11 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& tau_hat = 3.0_pReal*prm%b_tr(1)*mu/prm%L_tr & + Gamma_sf/(3.0_pReal*prm%b_tr(1)) & + prm%h*Delta_G/(3.0_pReal*prm%b_tr(1)) + x0 = mu*prm%b_sl(1)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu)) + tau_r = mu*prm%b_sl(1)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c)+cos(PI/3.0_pReal)/x0) do i = 1, prm%sum_N_tr tau = math_tensordot(Mp,prm%P_tr(1:3,1:3,i)) - x0 = mu*prm%b_tr(i)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu)) ! ToDo: In the paper, the Burgers vector for slip is used - tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(PI/3.0_pReal)/x0) ! ToDo: In the paper, the Burgers vector for slip is used if (tau > tol_math_check .and. tau < tau_r) then P = exp(-(tau_hat/tau)**prm%s(i)) From 8ceab5326f73bb25d66d75137ce5e4aab702a521 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 26 Jan 2022 22:49:11 +0100 Subject: [PATCH 04/17] more systematic handling of bcc/hcp case --- src/lattice.f90 | 39 ++++++++++++++++++++------------------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index e18c71edf..97f207024 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -528,22 +528,22 @@ end function lattice_C66_twin !> @brief Rotated elasticity matrices for transformation in 6x6-matrix notation !-------------------------------------------------------------------------------------------------- function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & - cOverA_trans,a_bcc,a_fcc) + cOverA_trans,a_fcc,a_bcc) integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) real(pReal), dimension(6,6), intent(in) :: C_parent66 + real(pReal), optional, intent(in) :: a_bcc, a_fcc, cOverA_trans real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66 real(pReal), dimension(3,3,sum(Ntrans)) :: Q,S type(rotation) :: R - real(pReal) :: a_bcc, a_fcc, cOverA_trans integer :: i !-------------------------------------------------------------------------------------------------- ! elasticity matrix of the target phase in cube orientation - if (lattice_target == 'hP') then + if (lattice_target == 'hP' .and. present(cOverA_trans)) then ! https://doi.org/10.1063/1.1663858 eq. (16), eq. (18), eq. (19) ! https://doi.org/10.1016/j.actamat.2016.07.032 eq. (47), eq. (48) if (cOverA_trans < 1.0_pReal .or. cOverA_trans > 2.0_pReal) & @@ -562,7 +562,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & C_target_unrotated66(3,3) = C_bar66(3,3) C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2))) C_target_unrotated66 = lattice_symmetrize_C66(C_target_unrotated66,'hP') - elseif (lattice_target == 'cI') then + elseif (lattice_target == 'cI' .and. present(a_fcc) .and. present(a_bcc)) then if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) & call IO_error(134,ext_msg='lattice_C66_trans: '//trim(lattice_target)) C_target_unrotated66 = C_parent66 @@ -1484,26 +1484,27 @@ end function lattice_SchmidMatrix_twin !> @brief Schmid matrix for transformation !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix) +function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_fcc,a_bcc) result(SchmidMatrix) integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) - real(pReal), intent(in) :: cOverA !< c/a ratio + real(pReal), optional, intent(in) :: cOverA, a_bcc, a_fcc real(pReal), dimension(3,3,sum(Ntrans)) :: SchmidMatrix real(pReal), dimension(3,3,sum(Ntrans)) :: devNull - real(pReal) :: a_bcc, a_fcc - if (lattice_target /= 'cI' .and. lattice_target /= 'hP') & - call IO_error(137,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) - if (lattice_target == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & + if (lattice_target == 'hP' .and. present(cOverA)) then + if (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal) & + call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) + call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA=cOverA) + else if (lattice_target == 'cI' .and. present(a_fcc) .and. present(a_bcc)) then + if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) & + call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) + call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,a_fcc=a_fcc,a_bcc=a_bcc) + else call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) - - if (lattice_target == 'cI' .and. (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal)) & - call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) - - call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA,a_fcc,a_bcc) + end if end function lattice_SchmidMatrix_trans @@ -1966,12 +1967,12 @@ end function buildCoordinateSystem !-------------------------------------------------------------------------------------------------- subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) - integer, dimension(:), intent(in) :: & + integer, dimension(:), intent(in) :: & Ntrans real(pReal), dimension(3,3,sum(Ntrans)), intent(out) :: & Q, & !< Total rotation: Q = R*B S !< Eigendeformation tensor for phase transformation - real(pReal), intent(in) :: & + real(pReal), optional, intent(in) :: & cOverA, & !< c/a for target hex lattice a_bcc, & !< lattice parameter a for bcc target lattice a_fcc !< lattice parameter a for fcc parent lattice @@ -2050,7 +2051,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) 0.0, 0.0, 1.0, 45.0 & ],shape(FCCTOBCC_BAINROT)) - if (a_bcc > 0.0_pReal .and. a_fcc > 0.0_pReal .and. dEq0(cOverA)) then ! fcc -> bcc + if (present(a_bcc) .and. present(a_fcc)) then do i = 1,sum(Ntrans) call R%fromAxisAngle(FCCTOBCC_SYSTEMTRANS(:,i),degrees=.true.,P=1) call B%fromAxisAngle(FCCTOBCC_BAINROT(:,i), degrees=.true.,P=1) @@ -2064,7 +2065,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) Q(1:3,1:3,i) = matmul(R%asMatrix(),B%asMatrix()) S(1:3,1:3,i) = matmul(R%asMatrix(),U) - MATH_I3 enddo - elseif (cOverA > 0.0_pReal .and. dEq0(a_bcc)) then ! fcc -> hex + else if (present(cOverA)) then ss = MATH_I3 sd = MATH_I3 ss(1,3) = sqrt(2.0_pReal)/4.0_pReal From 95fff5d6f705fe7e264558939841d3bab46019d9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 26 Jan 2022 23:01:47 +0100 Subject: [PATCH 05/17] limit to TWINNING to fcc --- src/phase_mechanical_plastic_dislotwin.f90 | 94 +++++++++------------- 1 file changed, 36 insertions(+), 58 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 36d2899b9..38ed6aec9 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -247,7 +247,7 @@ module function plastic_dislotwin_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! twin related parameters - prm%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(prm%N_tw)) twinActive: if (prm%sum_N_tw > 0) then prm%systems_tw = lattice_labels_twin(prm%N_tw,phase_lattice(ph)) @@ -264,40 +264,33 @@ module function plastic_dislotwin_init() result(myPlasticity) 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,prm%N_tw) - endif - ! expand: family => system 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 (.not. prm%fccTwinTransNucleation) extmsg = trim(extmsg)//' TWIP for non-fcc' if ( prm%L_tw < 0.0_pReal) extmsg = trim(extmsg)//' L_tw' if ( prm%i_tw < 0.0_pReal) extmsg = trim(extmsg)//' i_tw' if (any(prm%b_tw < 0.0_pReal)) extmsg = trim(extmsg)//' b_tw' if (any(prm%t_tw < 0.0_pReal)) extmsg = trim(extmsg)//' t_tw' if (any(prm%r < 0.0_pReal)) extmsg = trim(extmsg)//' p_tw' - if (.not. prm%fccTwinTransNucleation) then - if (any(prm%dot_N_0_tw < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tw' - end if else twinActive - allocate(prm%gamma_char,prm%b_tw,prm%dot_N_0_tw,prm%t_tw,prm%r,source=emptyRealArray) + allocate(prm%gamma_char,prm%b_tw,prm%t_tw,prm%r,source=emptyRealArray) allocate(prm%h_tw_tw(0,0)) end if twinActive !-------------------------------------------------------------------------------------------------- ! transformation related parameters - prm%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(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,prm%N_tr) - prm%h = pl%get_asFloat('h', defaultVal=0.0_pReal) ! ToDo: This is not optional! - prm%i_tr = pl%get_asFloat('i_tr', defaultVal=0.0_pReal) ! ToDo: This is not optional! + prm%h = pl%get_asFloat('h') + prm%i_tr = pl%get_asFloat('i_tr') prm%Delta_G(1) = pl%get_asFloat('Delta_G') prm%Delta_G(2) = pl%get_asFloat('Delta_G,T', defaultVal=0.0_pReal) prm%Delta_G(3) = pl%get_asFloat('Delta_G,T^2',defaultVal=0.0_pReal) @@ -317,10 +310,11 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%t_tr = pl%get_as1dFloat('t_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') prm%s = math_expand(prm%s,prm%N_tr) ! sanity checks + if (.not. prm%fccTwinTransNucleation) extmsg = trim(extmsg)//' TRIP for non-fcc' if ( prm%L_tr < 0.0_pReal) extmsg = trim(extmsg)//' L_tr' if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr' if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr' @@ -931,55 +925,39 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - isFCC: if (prm%fccTwinTransNucleation) then + mu = elastic_mu(ph,en) + nu = elastic_nu(ph,en) + Gamma_sf = prm%Gamma_sf(1) & + + prm%Gamma_sf(2) * (T-prm%T_ref) & + + prm%Gamma_sf(3) * (T-prm%T_ref)**2 + tau_hat = 3.0_pReal*prm%b_tw(1)*mu/prm%L_tw & + + Gamma_sf/(3.0_pReal*prm%b_tw(1)) + x0 = mu*prm%b_sl(1)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu)) + tau_r = mu*prm%b_sl(1)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c)+cos(PI/3.0_pReal)/x0) - mu = elastic_mu(ph,en) - nu = elastic_nu(ph,en) - Gamma_sf = prm%Gamma_sf(1) & - + prm%Gamma_sf(2) * (T-prm%T_ref) & - + prm%Gamma_sf(3) * (T-prm%T_ref)**2 - tau_hat = 3.0_pReal*prm%b_tw(1)*mu/prm%L_tw & - + Gamma_sf/(3.0_pReal*prm%b_tw(1)) - x0 = mu*prm%b_sl(1)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu)) - tau_r = mu*prm%b_sl(1)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c)+cos(PI/3.0_pReal)/x0) + do i = 1, prm%sum_N_tw + tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) - do i = 1, prm%sum_N_tw - tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) + if (tau > tol_math_check .and. tau < tau_r) then + P = exp(-(tau_hat/tau)**prm%r(i)) + dP_dTau = prm%r(i) * (tau_hat/tau)**prm%r(i)/tau * P - if (tau > tol_math_check .and. tau < tau_r) then - P = exp(-(tau_hat/tau)**prm%r(i)) - dP_dTau = prm%r(i) * (tau_hat/tau)**prm%r(i)/tau * P + s = prm%fcc_twinNucleationSlipPair(1:2,i) + dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) & + / (prm%L_tw*prm%b_sl(i)) - s = prm%fcc_twinNucleationSlipPair(1:2,i) - dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) & - / (prm%L_tw*prm%b_sl(i)) + P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) + dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) - P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) - dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) - - V = PI/4.0_pReal*dst%Lambda_tw(i,en)**2*prm%t_tw(i) - dot_gamma_tw(i) = V*dot_N_0*P_ncs*P - if (present(ddot_gamma_dtau_tw)) & - ddot_gamma_dtau_tw(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau) - else - dot_gamma_tw(i) = 0.0_pReal - if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal - end if - end do - - else isFCC - do i = 1, prm%sum_N_tw - error stop 'not implemented' - tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) - if (tau > tol_math_check) then - dot_gamma_tw(i) = 0.0_pReal - if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal - else - dot_gamma_tw(i) = 0.0_pReal - if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal - end if - end do - end if isFCC + V = PI/4.0_pReal*dst%Lambda_tw(i,en)**2*prm%t_tw(i) + dot_gamma_tw(i) = V*dot_N_0*P_ncs*P + if (present(ddot_gamma_dtau_tw)) & + ddot_gamma_dtau_tw(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau) + else + dot_gamma_tw(i) = 0.0_pReal + if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal + end if + end do end associate From 9e17efd1f89804de30dcbfece4bba0fccab8dbfd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 26 Jan 2022 23:13:56 +0100 Subject: [PATCH 06/17] transformation is for fcc-> hcp --- src/lattice.f90 | 2 +- src/phase_mechanical_plastic_dislotwin.f90 | 30 ++++++++++------------ 2 files changed, 14 insertions(+), 18 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 97f207024..b085baeae 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -533,7 +533,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) real(pReal), dimension(6,6), intent(in) :: C_parent66 - real(pReal), optional, intent(in) :: a_bcc, a_fcc, cOverA_trans + real(pReal), optional, intent(in) :: cOverA_trans, a_fcc, a_bcc real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66 diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 38ed6aec9..b46c894a2 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -27,8 +27,8 @@ submodule(phase:plastic) dislotwin E_sb = 1.0_pReal, & !< activation energy for shear bands h = 1.0_pReal, & !< stack height of hex nucleus T_ref = T_ROOM, & - a_cI = 1.0_pReal, & - a_cF = 1.0_pReal + a_cF = 1.0_pReal, & + cOverA_hP = 1.0_pReal real(pReal), dimension(3) :: & Gamma_sf = 0.0_pReal, & !< stacking fault energy Delta_G = 0.0_pReal !< free energy difference between austensite and martensite @@ -126,6 +126,7 @@ module function plastic_dislotwin_init() result(myPlasticity) startIndex, endIndex integer, dimension(:), allocatable :: & N_sl + real(pReal) :: a_cF real(pReal), allocatable, dimension(:) :: & rho_mob_0, & !< initial unipolar dislocation density per slip system rho_dip_0 !< initial dipole dislocation density per slip system @@ -289,24 +290,19 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%b_tr = pl%get_as1dFloat('b_tr') prm%b_tr = math_expand(prm%b_tr,prm%N_tr) - prm%h = pl%get_asFloat('h') - prm%i_tr = pl%get_asFloat('i_tr') - prm%Delta_G(1) = pl%get_asFloat('Delta_G') - prm%Delta_G(2) = pl%get_asFloat('Delta_G,T', defaultVal=0.0_pReal) - prm%Delta_G(3) = pl%get_asFloat('Delta_G,T^2',defaultVal=0.0_pReal) - 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%lattice_tr = pl%get_asString('lattice_tr') + prm%h = pl%get_asFloat('h') + prm%i_tr = pl%get_asFloat('i_tr') + prm%Delta_G(1) = pl%get_asFloat('Delta_G') + prm%Delta_G(2) = pl%get_asFloat('Delta_G,T', defaultVal=0.0_pReal) + prm%Delta_G(3) = pl%get_asFloat('Delta_G,T^2',defaultVal=0.0_pReal) + prm%L_tr = pl%get_asFloat('L_tr') + a_cF = pl%get_asFloat('a_cF') + prm%cOverA_hP = pl%get_asFloat('c/a_hP') prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,pl%get_as1dFloat('h_tr-tr'),& phase_lattice(ph)) - prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr,prm%lattice_tr, & - 0.0_pReal, & - prm%a_cI, & - prm%a_cF) + prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr,'hP',prm%cOverA_hP) prm%t_tr = pl%get_as1dFloat('t_tr') prm%t_tr = math_expand(prm%t_tr,prm%N_tr) @@ -472,7 +468,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) end if twinActive transActive: if (prm%sum_N_tr > 0) then - C66_tr = lattice_C66_trans(prm%N_tr,C,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF) + C66_tr = lattice_C66_trans(prm%N_tr,C,'hP',prm%cOverA_hP) do i=1,prm%sum_N_tr homogenizedC = homogenizedC & + stt%f_tr(i,en)*C66_tr(1:6,1:6,i) From 7d0d13bfed2893b1530a7a8ce1c01873fba503f8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 26 Jan 2022 23:25:47 +0100 Subject: [PATCH 07/17] following paper eq (10), eq (14), eq (49), eq (50), and Fig. 1 --- src/phase_mechanical_plastic_dislotwin.f90 | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index b46c894a2..b09303301 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -28,7 +28,9 @@ submodule(phase:plastic) dislotwin h = 1.0_pReal, & !< stack height of hex nucleus T_ref = T_ROOM, & a_cF = 1.0_pReal, & - cOverA_hP = 1.0_pReal + cOverA_hP = 1.0_pReal, & + V_mol = 1.0_pReal, & + rho = 1.0_pReal real(pReal), dimension(3) :: & Gamma_sf = 0.0_pReal, & !< stacking fault energy Delta_G = 0.0_pReal !< free energy difference between austensite and martensite @@ -290,15 +292,15 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%b_tr = pl%get_as1dFloat('b_tr') prm%b_tr = math_expand(prm%b_tr,prm%N_tr) - prm%h = pl%get_asFloat('h') prm%i_tr = pl%get_asFloat('i_tr') prm%Delta_G(1) = pl%get_asFloat('Delta_G') prm%Delta_G(2) = pl%get_asFloat('Delta_G,T', defaultVal=0.0_pReal) prm%Delta_G(3) = pl%get_asFloat('Delta_G,T^2',defaultVal=0.0_pReal) prm%L_tr = pl%get_asFloat('L_tr') a_cF = pl%get_asFloat('a_cF') + prm%h = 5.0_pReal * a_cF/sqrt(3.0_pReal) prm%cOverA_hP = pl%get_asFloat('c/a_hP') - + prm%rho = 4.0_pReal/(sqrt(3.0_pReal)*a_cF**2)/N_A prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,pl%get_as1dFloat('h_tr-tr'),& phase_lattice(ph)) @@ -1009,8 +1011,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& + prm%Delta_G(2) * (T-prm%T_ref) & + prm%Delta_G(3) * (T-prm%T_ref)**2 tau_hat = 3.0_pReal*prm%b_tr(1)*mu/prm%L_tr & - + Gamma_sf/(3.0_pReal*prm%b_tr(1)) & - + prm%h*Delta_G/(3.0_pReal*prm%b_tr(1)) + + (Gamma_sf + (prm%h/prm%V_mol - 2.0_pReal*prm%rho)*Delta_G)/(3.0_pReal*prm%b_tr(1)) x0 = mu*prm%b_sl(1)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu)) tau_r = mu*prm%b_sl(1)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c)+cos(PI/3.0_pReal)/x0) From 3bb2ae45d33d1e2a9f3ef72914c457161cc6a4c9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 26 Jan 2022 23:51:49 +0100 Subject: [PATCH 08/17] bugfixes - L_tw/tr are not in lengths of Burgers vectors anymore - need to define characteristic shear for tr --- src/phase_mechanical_plastic_dislotwin.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index b09303301..0fb885f1b 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -49,7 +49,7 @@ submodule(phase:plastic) dislotwin r, & !< exponent in twin nucleation rate s, & !< exponent in trans nucleation rate tau_0, & !< strength due to elements in solid solution - gamma_char, & !< characteristic shear for twins + gamma_char_tw, & !< characteristic shear for twins B, & !< drag coefficient d_caron !< distance of spontaneous annhihilation real(pReal), allocatable, dimension(:,:) :: & @@ -265,7 +265,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%L_tw = pl%get_asFloat('L_tw') prm%i_tw = pl%get_asFloat('i_tw') - prm%gamma_char= lattice_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph)) + prm%gamma_char_tw = lattice_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph)) ! expand: family => system prm%b_tw = math_expand(prm%b_tw,prm%N_tw) @@ -280,7 +280,7 @@ module function plastic_dislotwin_init() result(myPlasticity) if (any(prm%t_tw < 0.0_pReal)) extmsg = trim(extmsg)//' t_tw' if (any(prm%r < 0.0_pReal)) extmsg = trim(extmsg)//' p_tw' else twinActive - allocate(prm%gamma_char,prm%b_tw,prm%t_tw,prm%r,source=emptyRealArray) + allocate(prm%gamma_char_tw,prm%b_tw,prm%t_tw,prm%r,source=emptyRealArray) allocate(prm%h_tw_tw(0,0)) end if twinActive @@ -569,7 +569,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) shearBandingContribution: if (dNeq0(prm%v_sb)) then E_kB_T = prm%E_sb/(K_B*T) - call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? + call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? do i = 1,6 P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& @@ -626,6 +626,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) dot_gamma_tw real(pReal), dimension(param(ph)%sum_N_tr) :: & dot_gamma_tr + real(pReal), parameter :: & + gamma_char_tr = 1.0_pReal !!! TODO real(pReal) :: & mu, & nu, & @@ -686,10 +688,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) - dot_rho_dip_climb if (prm%sum_N_tw > 0) call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) - dot%f_tw(:,en) = f_matrix*dot_gamma_tw/prm%gamma_char + dot%f_tw(:,en) = f_matrix*dot_gamma_tw/prm%gamma_char_tw if (prm%sum_N_tr > 0) call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) - dot%f_tr(:,en) = f_matrix*dot_gamma_tr + dot%f_tr(:,en) = f_matrix*dot_gamma_tr/gamma_char_tr end associate @@ -941,8 +943,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& dP_dTau = prm%r(i) * (tau_hat/tau)**prm%r(i)/tau * P s = prm%fcc_twinNucleationSlipPair(1:2,i) - dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) & - / (prm%L_tw*prm%b_sl(i)) + dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/prm%L_tw P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) @@ -1023,8 +1024,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& dP_dTau = prm%s(i) * (tau_hat/tau)**prm%s(i)/tau * P s = prm%fcc_twinNucleationSlipPair(1:2,i) - dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) & - / (prm%L_tr*prm%b_sl(i)) + dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/prm%L_tr P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) From e7ed32ba4ebaff128f49d6fee52ec0ddabee65a9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 27 Jan 2022 00:01:53 +0100 Subject: [PATCH 09/17] not needed --- src/phase_mechanical_plastic.f90 | 6 ++---- src/phase_mechanical_plastic_dislotwin.f90 | 4 +--- 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 0a24b0cbf..f72251bad 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -167,12 +167,10 @@ submodule(phase:mechanical) plastic el !< current element number end subroutine nonlocal_dotState - module subroutine dislotwin_dependentState(T,ph,en) + module subroutine dislotwin_dependentState(ph,en) integer, intent(in) :: & ph, & en - real(pReal), intent(in) :: & - T end subroutine dislotwin_dependentState module subroutine dislotungsten_dependentState(ph,en) @@ -363,7 +361,7 @@ module subroutine plastic_dependentState(co, ip, el) plasticType: select case (phase_plasticity(ph)) case (PLASTIC_DISLOTWIN_ID) plasticType - call dislotwin_dependentState(thermal_T(ph,en),ph,en) + call dislotwin_dependentState(ph,en) case (PLASTIC_DISLOTUNGSTEN_ID) plasticType call dislotungsten_dependentState(ph,en) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 0fb885f1b..d886322ac 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -701,13 +701,11 @@ end subroutine dislotwin_dotState !-------------------------------------------------------------------------------------------------- !> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- -module subroutine dislotwin_dependentState(T,ph,en) +module subroutine dislotwin_dependentState(ph,en) integer, intent(in) :: & ph, & en - real(pReal), intent(in) :: & - T real(pReal) :: & sumf_tw, sumf_tr From 1e748dc54ebd7bc5e39b65bd7826a635896d2806 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 27 Jan 2022 07:28:33 +0100 Subject: [PATCH 10/17] commenting --- src/phase_mechanical_plastic_dislotwin.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index d886322ac..8660ee1ad 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -627,7 +627,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) real(pReal), dimension(param(ph)%sum_N_tr) :: & dot_gamma_tr real(pReal), parameter :: & - gamma_char_tr = 1.0_pReal !!! TODO + gamma_char_tr = 1.0_pReal ! ToDo: get correct value real(pReal) :: & mu, & nu, & @@ -941,7 +941,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& dP_dTau = prm%r(i) * (tau_hat/tau)**prm%r(i)/tau * P s = prm%fcc_twinNucleationSlipPair(1:2,i) - dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/prm%L_tw + dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/prm%L_tw ! ToDo: factor 3 P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) @@ -1022,7 +1022,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& dP_dTau = prm%s(i) * (tau_hat/tau)**prm%s(i)/tau * P s = prm%fcc_twinNucleationSlipPair(1:2,i) - dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/prm%L_tr + dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/prm%L_tr ! ToDo: factor 3 P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) From ea33b321aa8db723ba814da938f42f6389fa42f3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 27 Jan 2022 08:08:42 +0100 Subject: [PATCH 11/17] correct calculation of shear need to convert volume to shear (and back for output) --- src/phase_mechanical_plastic_dislotwin.f90 | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 8660ee1ad..046d8f2ca 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -27,6 +27,7 @@ submodule(phase:plastic) dislotwin E_sb = 1.0_pReal, & !< activation energy for shear bands h = 1.0_pReal, & !< stack height of hex nucleus T_ref = T_ROOM, & + gamma_char_tr = sqrt(0.125_pReal), & !< Characteristic shear for transformation a_cF = 1.0_pReal, & cOverA_hP = 1.0_pReal, & V_mol = 1.0_pReal, & @@ -626,8 +627,6 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) dot_gamma_tw real(pReal), dimension(param(ph)%sum_N_tr) :: & dot_gamma_tr - real(pReal), parameter :: & - gamma_char_tr = 1.0_pReal ! ToDo: get correct value real(pReal) :: & mu, & nu, & @@ -691,7 +690,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) dot%f_tw(:,en) = f_matrix*dot_gamma_tw/prm%gamma_char_tw if (prm%sum_N_tr > 0) call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) - dot%f_tr(:,en) = f_matrix*dot_gamma_tr/gamma_char_tr + dot%f_tr(:,en) = f_matrix*dot_gamma_tr/prm%gamma_char_tr end associate @@ -947,9 +946,9 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) V = PI/4.0_pReal*dst%Lambda_tw(i,en)**2*prm%t_tw(i) - dot_gamma_tw(i) = V*dot_N_0*P_ncs*P + dot_gamma_tw(i) = V*dot_N_0*P_ncs*P*prm%gamma_char_tw(i) if (present(ddot_gamma_dtau_tw)) & - ddot_gamma_dtau_tw(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau) + ddot_gamma_dtau_tw(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)*prm%gamma_char_tw(i) else dot_gamma_tw(i) = 0.0_pReal if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal @@ -1028,9 +1027,9 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) V = PI/4.0_pReal*dst%Lambda_tr(i,en)**2*prm%t_tr(i) - dot_gamma_tr(i) = V*dot_N_0*P_ncs*P + dot_gamma_tr(i) = V*dot_N_0*P_ncs*P*prm%gamma_char_tr if (present(ddot_gamma_dtau_tr)) & - ddot_gamma_dtau_tr(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau) + ddot_gamma_dtau_tr(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)*prm%gamma_char_tr else dot_gamma_tr(i) = 0.0_pReal if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr(i) = 0.0_pReal From 20b2079135c85e1c6e23975c313b0971519a561e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 6 Feb 2022 08:10:23 +0100 Subject: [PATCH 12/17] need to read V_mol, default of 1.0 gives strange results --- src/phase_mechanical_plastic_dislotwin.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 046d8f2ca..016e9baa8 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -302,6 +302,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%h = 5.0_pReal * a_cF/sqrt(3.0_pReal) prm%cOverA_hP = pl%get_asFloat('c/a_hP') prm%rho = 4.0_pReal/(sqrt(3.0_pReal)*a_cF**2)/N_A + prm%V_mol = pl%get_asFloat('V_mol') prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,pl%get_as1dFloat('h_tr-tr'),& phase_lattice(ph)) @@ -315,6 +316,7 @@ module function plastic_dislotwin_init() result(myPlasticity) ! sanity checks if (.not. prm%fccTwinTransNucleation) extmsg = trim(extmsg)//' TRIP for non-fcc' if ( prm%L_tr < 0.0_pReal) extmsg = trim(extmsg)//' L_tr' + if ( prm%V_mol < 0.0_pReal) extmsg = trim(extmsg)//' V_mol' if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr' if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr' if (any(prm%s < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr' @@ -1025,7 +1027,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) - + V = PI/4.0_pReal*dst%Lambda_tr(i,en)**2*prm%t_tr(i) dot_gamma_tr(i) = V*dot_N_0*P_ncs*P*prm%gamma_char_tr if (present(ddot_gamma_dtau_tr)) & From 1883919b3cd9d8f18ecfccfd7af7f1b2e0490931 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 6 Feb 2022 08:11:05 +0100 Subject: [PATCH 13/17] working configuration simulations are not particular stable, but trends are ok and values seem reasonable --- .../Phase_Dislotwin_TWIP-Steel-FeMnC.yaml | 73 +++++++++++++------ 1 file changed, 52 insertions(+), 21 deletions(-) diff --git a/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml b/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml index 6c19b9f93..e0f9e4cc9 100644 --- a/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml +++ b/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml @@ -1,31 +1,62 @@ type: dislotwin -output: [rho_mob, rho_dip, gamma_sl, Lambda_sl, tau_pass, f_tw, Lambda_tw, tau_hat_tw, f_tr] -D: 2.0e-5 +output: [rho_mob, rho_dip, gamma_sl, Lambda_sl, tau_pass, f_tw, Lambda_tw, f_tr] + +# Slip N_sl: [12] b_sl: [2.56e-10] +Q_sl: [3.5e-19] +p_sl: [0.325] +q_sl: [1.55] +B: [0.001] +i_sl: [30.0] + rho_mob_0: [1.0e+12] rho_dip_0: [1.0] + v_0: [1.0e+4] -Q_sl: [3.7e-19] -p_sl: [1.0] -q_sl: [1.0] -tau_0: [1.5e+8] -i_sl: [10.0] # Adj. parameter controlling dislocation mean free path -D_0: 4.0e-5 # Vacancy diffusion prefactor / m^2/s +tau_0: [3.5e+8] D_a: 1.0 # minimum dipole distance / b Q_cl: 4.5e-19 # Activation energy for climb / J + h_sl-sl: [0.122, 0.122, 0.625, 0.07, 0.137, 0.137, 0.122] # Interaction coefficients (Kubin et al. 2008) -# twinning parameters +h_tw-sl: [0.0] # ToDo: values not known + +# Twin N_tw: [12] -b_tw: [1.47e-10] # Burgers vector length of twin system / b -t_tw: [5.0e-8] # Twin stack mean thickness / m -L_tw: 442.0 # Length of twin nuclei / b -x_c_tw: 1.0e-9 # critical distance for formation of twin nucleus / m -V_cs: 1.67e-29 # cross slip volume / m^3 -p_tw: [10.0] # r-exponent in twin formation probability -i_tw: 1.0 # Adj. parameter controlling twin mean free path -h_sl-tw: [0.0, 1.0, 1.0] # dislocation-twin interaction coefficients -h_tw-tw: [0.0, 1.0] # twin-twin interaction coefficients -T_ref: 0.0 -Gamma_sf: -0.0396 # stacking fault energy / J/m^2 at zero K; TWIP steel: -0.0526; Cu: -0.0396 -Gamma_sf,T: 0.0002 # temperature dependence / J/(m^2 K) of stacking fault energy +b_tw: [1.40e-10] +L_tw: 2.e-7 # 800 b_sl +i_tw: 10.0 +t_tw: [5.0e-8] +p_tw: [2.5] # A + +h_tw-tw: [0.0, 1.0] # ToDo: values not known +h_sl-tw: [0.0, 1.0, 1.0] # ToDo: values not known + +# Transformation +N_tr: [12] +b_tr: [1.40e-10] +L_tr: 2.e-7 # 800 b_sl +i_tr: 3.0 +t_tr: [1.0e-7] +p_tr: [2.5] # B +V_mol: 7.09e-6 +Delta_G: 1.33343932e+02 +Delta_G,T: 2.56640412 +Delta_G,T^2: 1.49524179e-03 + +h_tr-tr: [1.0, 1.0] # guessing +h_sl-tr: [0.0, 1.0, 1.0] # guessing + +# Twin & Transformation +T_ref: 298.15 +Gamma_sf: 2.89394017e-02 +Gamma_sf,T: 1.22823814e-04 +Gamma_sf,T^2: 1.47322968e-07 +x_c: 1.0e-9 +V_cs: 1.67e-29 + +a_cF: 3.62e-10 +c/a_hP: 1.633 + +# Slip & Twin & Transformation +D: 2.0e-5 From 59fdd9b58625f1f149b1d7c3addf2545acf49bb8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 6 Feb 2022 10:07:15 +0100 Subject: [PATCH 14/17] in agreement with paper this is the sytem-wise formulation of D. Steinmetz' model --- src/phase_mechanical_plastic_dislotwin.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 016e9baa8..e4afd305a 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -942,7 +942,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& dP_dTau = prm%r(i) * (tau_hat/tau)**prm%r(i)/tau * P s = prm%fcc_twinNucleationSlipPair(1:2,i) - dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/prm%L_tw ! ToDo: factor 3 + dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/(prm%L_tw*3.0_pReal) P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) @@ -1023,7 +1023,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& dP_dTau = prm%s(i) * (tau_hat/tau)**prm%s(i)/tau * P s = prm%fcc_twinNucleationSlipPair(1:2,i) - dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/prm%L_tr ! ToDo: factor 3 + dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/(prm%L_tr*3.0_pReal) P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) From 8f0a16e6549230ab7ecc8ec28f9f16cafeb9e21e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 6 Feb 2022 11:27:53 +0100 Subject: [PATCH 15/17] fine tuning parameters more stable, still physically reasonable --- examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml b/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml index e0f9e4cc9..864a8d2bd 100644 --- a/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml +++ b/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml @@ -24,10 +24,10 @@ h_tw-sl: [0.0] # ToDo: values not known # Twin N_tw: [12] b_tw: [1.40e-10] -L_tw: 2.e-7 # 800 b_sl +L_tw: 1.8e-7 i_tw: 10.0 t_tw: [5.0e-8] -p_tw: [2.5] # A +p_tw: [3.0] # A h_tw-tw: [0.0, 1.0] # ToDo: values not known h_sl-tw: [0.0, 1.0, 1.0] # ToDo: values not known @@ -35,10 +35,10 @@ h_sl-tw: [0.0, 1.0, 1.0] # ToDo: values not known # Transformation N_tr: [12] b_tr: [1.40e-10] -L_tr: 2.e-7 # 800 b_sl +L_tr: 1.8e-7 i_tr: 3.0 t_tr: [1.0e-7] -p_tr: [2.5] # B +p_tr: [3.0] # B V_mol: 7.09e-6 Delta_G: 1.33343932e+02 Delta_G,T: 2.56640412 From 7c18f3f276f24ff4ecfc79b445287c75d66738fa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 11 Feb 2022 08:30:24 +0100 Subject: [PATCH 16/17] simplified --- src/phase_mechanical_elastic.f90 | 2 +- src/phase_mechanical_plastic_dislotwin.f90 | 59 +++++++++------------- 2 files changed, 24 insertions(+), 37 deletions(-) diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 7f30c8165..84ad7a20e 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -48,7 +48,7 @@ module subroutine elastic_init(phases) prm%C_11 = polynomial(elastic%asDict(),'C_11','T') prm%C_12 = polynomial(elastic%asDict(),'C_12','T') prm%C_44 = polynomial(elastic%asDict(),'C_44','T') - + if (any(phase_lattice(ph) == ['hP','tI'])) then prm%C_13 = polynomial(elastic%asDict(),'C_13','T') prm%C_33 = polynomial(elastic%asDict(),'C_33','T') diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index e4afd305a..a81860be8 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -26,15 +26,14 @@ submodule(phase:plastic) dislotwin v_sb = 1.0_pReal, & !< value for shearband velocity_0 E_sb = 1.0_pReal, & !< activation energy for shear bands h = 1.0_pReal, & !< stack height of hex nucleus - T_ref = T_ROOM, & gamma_char_tr = sqrt(0.125_pReal), & !< Characteristic shear for transformation a_cF = 1.0_pReal, & cOverA_hP = 1.0_pReal, & V_mol = 1.0_pReal, & rho = 1.0_pReal - real(pReal), dimension(3) :: & - Gamma_sf = 0.0_pReal, & !< stacking fault energy - Delta_G = 0.0_pReal !< free energy difference between austensite and martensite + type(tPolynomial) :: & + Gamma_sf, & !< stacking fault energy + Delta_G !< free energy difference between austensite and martensite real(pReal), allocatable, dimension(:) :: & b_sl, & !< absolute length of Burgers vector [m] for each slip system b_tw, & !< absolute length of Burgers vector [m] for each twin system @@ -79,7 +78,7 @@ submodule(phase:plastic) dislotwin character(len=pStringLen), allocatable, dimension(:) :: & output logical :: & - ExtendedDislocations, & !< consider split into partials for climb calculation + extendedDislocations, & !< consider split into partials for climb calculation fccTwinTransNucleation, & !< twinning and transformation models are for fcc omitDipoles !< flag controlling consideration of dipole formation character(len=:), allocatable, dimension(:) :: & @@ -210,7 +209,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%Q_cl = pl%get_asFloat('Q_cl') - prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.) + prm%extendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.) prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.) ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) @@ -294,9 +293,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%b_tr = math_expand(prm%b_tr,prm%N_tr) prm%i_tr = pl%get_asFloat('i_tr') - prm%Delta_G(1) = pl%get_asFloat('Delta_G') - prm%Delta_G(2) = pl%get_asFloat('Delta_G,T', defaultVal=0.0_pReal) - prm%Delta_G(3) = pl%get_asFloat('Delta_G,T^2',defaultVal=0.0_pReal) + prm%Delta_G = polynomial(pl%asDict(),'Delta_G','T') prm%L_tr = pl%get_asFloat('L_tr') a_cF = pl%get_asFloat('a_cF') prm%h = 5.0_pReal * a_cF/sqrt(3.0_pReal) @@ -353,12 +350,8 @@ module function plastic_dislotwin_init() result(myPlasticity) if (prm%V_cs < 0.0_pReal) extmsg = trim(extmsg)//' V_cs' end if - if (prm%sum_N_tw + prm%sum_N_tr > 0 .or. prm%ExtendedDislocations) then - prm%T_ref = pl%get_asFloat('T_ref') - prm%Gamma_sf(1) = pl%get_asFloat('Gamma_sf') - prm%Gamma_sf(2) = pl%get_asFloat('Gamma_sf,T',defaultVal=0.0_pReal) - prm%Gamma_sf(3) = pl%get_asFloat('Gamma_sf,T^2',defaultVal=0.0_pReal) - end if + if (prm%sum_N_tw + prm%sum_N_tr > 0 .or. prm%extendedDislocations) & + prm%Gamma_sf = polynomial(pl%asDict(),'Gamma_sf','T') slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dFloat('h_sl-tw'), & @@ -631,8 +624,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) dot_gamma_tr real(pReal) :: & mu, & - nu, & - Gamma_sf + nu + associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) @@ -664,13 +657,12 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) dot_rho_dip_climb(i) = 0.0_pReal else ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 - Gamma_sf = prm%Gamma_sf(1) & - + prm%Gamma_sf(2) * (T-prm%T_ref) & - + prm%Gamma_sf(3) * (T-prm%T_ref)**2 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 - nu)/(2.0_pReal + nu) * Gamma_sf / (mu*prm%b_sl(i)), & - 1.0_pReal, & - prm%ExtendedDislocations) + if (prm%extendedDislocations) then + b_d = 24.0_pReal*PI*(1.0_pReal - nu)/(2.0_pReal + nu) * prm%Gamma_sf%at(T) / (mu*prm%b_sl(i)) + else + b_d = 1.0_pReal + end if v_cl = 2.0_pReal*prm%omega*b_d**2*exp(-prm%Q_cl/(K_B*T)) & * (exp(abs(sigma_cl)*prm%b_sl(i)**3/(K_B*T)) - 1.0_pReal) @@ -719,8 +711,8 @@ module subroutine dislotwin_dependentState(ph,en) inv_lambda_tr_tr, & !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite f_over_t_tr real(pReal) :: & - mu, & - nu + mu + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) @@ -926,9 +918,8 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) - Gamma_sf = prm%Gamma_sf(1) & - + prm%Gamma_sf(2) * (T-prm%T_ref) & - + prm%Gamma_sf(3) * (T-prm%T_ref)**2 + Gamma_sf = prm%Gamma_sf%at(T) + tau_hat = 3.0_pReal*prm%b_tw(1)*mu/prm%L_tw & + Gamma_sf/(3.0_pReal*prm%b_tw(1)) x0 = mu*prm%b_sl(1)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu)) @@ -991,7 +982,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& tau, tau_r, tau_hat, & dot_N_0, & x0, V, & - Gamma_sf, Delta_G, & + Gamma_sf, & mu, nu, & P_ncs, dP_ncs_dtau, & P, dP_dtau @@ -1004,14 +995,10 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) - Gamma_sf = prm%Gamma_sf(1) & - + prm%Gamma_sf(2) * (T-prm%T_ref) & - + prm%Gamma_sf(3) * (T-prm%T_ref)**2 - Delta_G = prm%Delta_G(1) & - + prm%Delta_G(2) * (T-prm%T_ref) & - + prm%Delta_G(3) * (T-prm%T_ref)**2 + Gamma_sf = prm%Gamma_sf%at(T) + tau_hat = 3.0_pReal*prm%b_tr(1)*mu/prm%L_tr & - + (Gamma_sf + (prm%h/prm%V_mol - 2.0_pReal*prm%rho)*Delta_G)/(3.0_pReal*prm%b_tr(1)) + + (Gamma_sf + (prm%h/prm%V_mol - 2.0_pReal*prm%rho)*prm%Delta_G%at(T))/(3.0_pReal*prm%b_tr(1)) x0 = mu*prm%b_sl(1)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu)) tau_r = mu*prm%b_sl(1)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c)+cos(PI/3.0_pReal)/x0) From 9d1cc93a81e64156498f1d4cdaacb81c20c910cc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 11 Feb 2022 14:23:41 +0100 Subject: [PATCH 17/17] no need to store dotstate --- src/phase_mechanical_plastic.f90 | 11 +-- src/phase_mechanical_plastic_dislotwin.f90 | 106 +++++++++++---------- 2 files changed, 63 insertions(+), 54 deletions(-) diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 92530496f..22333560c 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -140,15 +140,15 @@ submodule(phase:mechanical) plastic dotState end function plastic_kinehardening_dotState - module subroutine dislotwin_dotState(Mp,T,ph,en) + module function dislotwin_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - real(pReal), intent(in) :: & - T integer, intent(in) :: & ph, & en - end subroutine dislotwin_dotState + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState + end function dislotwin_dotState module function dislotungsten_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & @@ -331,8 +331,7 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) dotState = plastic_kinehardening_dotState(Mp,ph,en) case (PLASTIC_DISLOTWIN_ID) plasticType - call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en) - dotState = plasticState(ph)%dotState(:,en) + dotState = dislotwin_dotState(Mp,ph,en) case (PLASTIC_DISLOTUNGSTEN_ID) plasticType dotState = dislotungsten_dotState(Mp,ph,en) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index a81860be8..37c86c137 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -84,7 +84,16 @@ submodule(phase:plastic) dislotwin character(len=:), allocatable, dimension(:) :: & systems_sl, & systems_tw - end type !< container type for internal constitutive parameters + end type tParameters !< container type for internal constitutive parameters + + type :: tIndexDotState + integer, dimension(2) :: & + rho_mob, & + rho_dip, & + gamma_sl, & + f_tw, & + f_tr + end type tIndexDotState type :: tDislotwinState real(pReal), dimension(:,:), pointer :: & @@ -106,9 +115,8 @@ submodule(phase:plastic) dislotwin !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param - type(tDislotwinState), allocatable, dimension(:) :: & - dotState, & - state + type(tIndexDotState), allocatable, dimension(:) :: indexDotState + type(tDislotwinState), allocatable, dimension(:) :: state type(tDislotwinDependentState), allocatable, dimension(:) :: dependentState contains @@ -159,15 +167,15 @@ module function plastic_dislotwin_init() result(myPlasticity) phases => config_material%get('phase') allocate(param(phases%length)) + allocate(indexDotState(phases%length)) allocate(state(phases%length)) - allocate(dotState(phases%length)) allocate(dependentState(phases%length)) - do ph = 1, phases%length if (.not. myPlasticity(ph)) cycle - associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), & + idx_dot => indexDotState(ph)) phase => phases%get(ph) mech => phase%get('mechanical') @@ -179,7 +187,6 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray) #endif - !-------------------------------------------------------------------------------------------------- ! slip related parameters N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) @@ -373,51 +380,51 @@ module function plastic_dislotwin_init() result(myPlasticity) + size(['f_tr']) * prm%sum_N_tr sizeState = sizeDotState - call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) + deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely !-------------------------------------------------------------------------------------------------- -! locally defined state aliases and initialization of state0 and atol +! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl + idx_dot%rho_mob = [startIndex,endIndex] stt%rho_mob=>plasticState(ph)%state(startIndex:endIndex,:) stt%rho_mob= spread(rho_mob_0,2,Nmembers) - dot%rho_mob=>plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl + idx_dot%rho_dip = [startIndex,endIndex] stt%rho_dip=>plasticState(ph)%state(startIndex:endIndex,:) stt%rho_dip= spread(rho_dip_0,2,Nmembers) - dot%rho_dip=>plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl + idx_dot%gamma_sl = [startIndex,endIndex] stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:) - dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw + idx_dot%f_tw = [startIndex,endIndex] stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:) - dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tw',defaultVal=1.0e-6_pReal) if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tw' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tr + idx_dot%f_tr = [startIndex,endIndex] stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:) - dot%f_tr=>plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tr',defaultVal=1.0e-6_pReal) if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tr' - allocate(dst%tau_pass (prm%sum_N_sl,Nmembers),source=0.0_pReal) - allocate(dst%Lambda_sl (prm%sum_N_sl,Nmembers),source=0.0_pReal) - allocate(dst%Lambda_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal) - allocate(dst%Lambda_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal) + allocate(dst%tau_pass (prm%sum_N_sl,Nmembers),source=0.0_pReal) + allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers),source=0.0_pReal) + allocate(dst%Lambda_tw(prm%sum_N_tw,Nmembers),source=0.0_pReal) + allocate(dst%Lambda_tr(prm%sum_N_tr,Nmembers),source=0.0_pReal) end associate @@ -596,15 +603,15 @@ end subroutine dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine dislotwin_dotState(Mp,T,ph,en) +module function dislotwin_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in):: & Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature at integration point integer, intent(in) :: & ph, & en + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState integer :: i real(pReal) :: & @@ -623,21 +630,27 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) real(pReal), dimension(param(ph)%sum_N_tr) :: & dot_gamma_tr real(pReal) :: & - mu, & - nu + mu, nu, & + T - associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), & + dot_rho_mob => dotState(indexDotState(ph)%rho_mob(1):indexDotState(ph)%rho_mob(2)), & + dot_rho_dip => dotState(indexDotState(ph)%rho_dip(1):indexDotState(ph)%rho_dip(2)), & + abs_dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), & + dot_f_tw => dotState(indexDotState(ph)%f_tw(1):indexDotState(ph)%f_tw(2)), & + dot_f_tr => dotState(indexDotState(ph)%f_tr(1):indexDotState(ph)%f_tr(2))) mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) + T = thermal_T(ph,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) + abs_dot_gamma_sl = abs(dot_gamma_sl) slipState: do i = 1, prm%sum_N_sl tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) @@ -651,7 +664,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) d_hat = math_clip(d_hat, left = prm%d_caron(i)) dot_rho_dip_formation(i) = 2.0_pReal*(d_hat-prm%d_caron(i))/prm%b_sl(i) & - * stt%rho_mob(i,en)*abs(dot_gamma_sl(i)) + * stt%rho_mob(i,en)*abs_dot_gamma_sl(i) if (dEq(d_hat,prm%d_caron(i))) then dot_rho_dip_climb(i) = 0.0_pReal @@ -672,23 +685,23 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) end if significantSlipStress end do slipState - dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) & - - dot_rho_dip_formation & - - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl) + dot_rho_mob = abs_dot_gamma_sl/(prm%b_sl*dst%Lambda_sl(:,en)) & + - dot_rho_dip_formation & + - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_mob(:,en)*abs_dot_gamma_sl - dot%rho_dip(:,en) = dot_rho_dip_formation & - - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & - - dot_rho_dip_climb + dot_rho_dip = dot_rho_dip_formation & + - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_dip(:,en)*abs_dot_gamma_sl & + - dot_rho_dip_climb - if (prm%sum_N_tw > 0) call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) - dot%f_tw(:,en) = f_matrix*dot_gamma_tw/prm%gamma_char_tw + if (prm%sum_N_tw > 0) call kinetics_tw(Mp,T,abs_dot_gamma_sl,ph,en,dot_gamma_tw) + dot_f_tw = f_matrix*dot_gamma_tw/prm%gamma_char_tw - if (prm%sum_N_tr > 0) call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) - dot%f_tr(:,en) = f_matrix*dot_gamma_tr/prm%gamma_char_tr + if (prm%sum_N_tr > 0) call kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,dot_gamma_tr) + dot_f_tr = f_matrix*dot_gamma_tr/prm%gamma_char_tr end associate -end subroutine dislotwin_dotState +end function dislotwin_dotState !-------------------------------------------------------------------------------------------------- @@ -816,15 +829,14 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & integer, intent(in) :: & ph, & en - real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: & dot_gamma_sl real(pReal), dimension(param(ph)%sum_N_sl), optional, intent(out) :: & ddot_gamma_dtau_sl, & tau_sl + real(pReal), dimension(param(ph)%sum_N_sl) :: & ddot_gamma_dtau - real(pReal), dimension(param(ph)%sum_N_sl) :: & tau, & stressRatio, & @@ -883,7 +895,7 @@ end subroutine kinetics_sl ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& +pure subroutine kinetics_tw(Mp,T,abs_dot_gamma_sl,ph,en,& dot_gamma_tw,ddot_gamma_dtau_tw) real(pReal), dimension(3,3), intent(in) :: & @@ -894,8 +906,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& ph, & en real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & - dot_gamma_sl - + abs_dot_gamma_sl real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & dot_gamma_tw real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: & @@ -933,7 +944,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& dP_dTau = prm%r(i) * (tau_hat/tau)**prm%r(i)/tau * P s = prm%fcc_twinNucleationSlipPair(1:2,i) - dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/(prm%L_tw*3.0_pReal) + dot_N_0 = sum(abs_dot_gamma_sl(s(2:1:-1))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/(prm%L_tw*3.0_pReal) P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) @@ -960,7 +971,7 @@ end subroutine kinetics_tw ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& +pure subroutine kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,& dot_gamma_tr,ddot_gamma_dtau_tr) real(pReal), dimension(3,3), intent(in) :: & @@ -971,8 +982,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& ph, & en real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & - dot_gamma_sl - + abs_dot_gamma_sl real(pReal), dimension(param(ph)%sum_N_tr), intent(out) :: & dot_gamma_tr real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: & @@ -1010,7 +1020,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& dP_dTau = prm%s(i) * (tau_hat/tau)**prm%s(i)/tau * P s = prm%fcc_twinNucleationSlipPair(1:2,i) - dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/(prm%L_tr*3.0_pReal) + dot_N_0 = sum(abs_dot_gamma_sl(s(2:1:-1))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/(prm%L_tr*3.0_pReal) P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)