From 6936ecb0910066fbfa59fcd19331a93fae4caa50 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 21:36:26 +0200 Subject: [PATCH 01/10] sl is the subscript for slip --- .../Phase_Dislotwin_TWIP-Steel-FeMnC.yaml | 2 +- .../plastic/dislotwin_IF-steel.yaml | 2 +- src/phase_mechanical_plastic_dislotwin.f90 | 48 +++++++++---------- 3 files changed, 25 insertions(+), 27 deletions(-) diff --git a/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml b/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml index fe51585be..6c19b9f93 100644 --- a/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml +++ b/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml @@ -6,7 +6,7 @@ b_sl: [2.56e-10] rho_mob_0: [1.0e+12] rho_dip_0: [1.0] v_0: [1.0e+4] -Q_s: [3.7e-19] +Q_sl: [3.7e-19] p_sl: [1.0] q_sl: [1.0] tau_0: [1.5e+8] diff --git a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml index 5478aa869..e48d5437f 100644 --- a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml +++ b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml @@ -11,7 +11,7 @@ b_sl: [2.49e-10, 2.49e-10] rho_mob_0: [2.81e12, 2.8e+12] rho_dip_0: [1.0, 1.0] # not given v_0: [1.4e+3, 1.4e+3] -Q_s: [1.57e-19, 1.57e-19] # Delta_F +Q_sl: [1.57e-19, 1.57e-19] # Delta_F tau_0: [454.e+6, 454.e+6] p_sl: [0.325, 0.325] q_sl: [1.55, 1.55] diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 6ecf9ee7a..4aa62bcb8 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -22,7 +22,6 @@ submodule(phase:plastic) dislotwin D = 1.0_pReal, & !< grain size p_sb = 1.0_pReal, & !< p-exponent in shear band velocity q_sb = 1.0_pReal, & !< q-exponent in shear band velocity - D_a = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors @@ -42,7 +41,7 @@ submodule(phase:plastic) dislotwin b_sl, & !< absolute length of Burgers vector [m] for each slip system b_tw, & !< absolute length of Burgers vector [m] for each twin system b_tr, & !< absolute length of Burgers vector [m] for each transformation system - Q_s,& !< activation energy for glide [J] for each slip system + Q_sl,& !< activation energy for glide [J] for each slip system v_0, & !< dislocation velocity prefactor [m/s] for each slip system dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system @@ -55,7 +54,8 @@ submodule(phase:plastic) dislotwin s, & !< s-exponent in trans nucleation rate tau_0, & !< strength due to elements in solid solution gamma_char, & !< characteristic shear for twins - B !< drag coefficient + B, & !< drag coefficient + d_caron !< distance of spontaneous annhihilation real(pReal), allocatable, dimension(:,:) :: & h_sl_sl, & !< components of slip-slip interaction matrix h_sl_tw, & !< components of slip-twin interaction matrix @@ -206,7 +206,7 @@ module function plastic_dislotwin_init() result(myPlasticity) rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl)) prm%v_0 = pl%get_as1dFloat('v_0', requiredSize=size(N_sl)) prm%b_sl = pl%get_as1dFloat('b_sl', requiredSize=size(N_sl)) - prm%Q_s = pl%get_as1dFloat('Q_s', requiredSize=size(N_sl)) + prm%Q_sl = pl%get_as1dFloat('Q_sl', requiredSize=size(N_sl)) prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(N_sl)) prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl)) prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl)) @@ -214,9 +214,8 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%B = pl%get_as1dFloat('B', requiredSize=size(N_sl), & defaultVal=[(0.0_pReal, i=1,size(N_sl))]) - prm%D_a = pl%get_asFloat('D_a') - prm%D_0 = pl%get_asFloat('D_0') - prm%Q_cl = pl%get_asFloat('Q_cl') + prm%D_0 = pl%get_asFloat('D_0') + prm%Q_cl = pl%get_asFloat('Q_cl') prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.) prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.) @@ -231,12 +230,13 @@ module function plastic_dislotwin_init() result(myPlasticity) rho_dip_0 = math_expand(rho_dip_0, N_sl) prm%v_0 = math_expand(prm%v_0, N_sl) prm%b_sl = math_expand(prm%b_sl, N_sl) - prm%Q_s = math_expand(prm%Q_s, N_sl) + prm%Q_sl = math_expand(prm%Q_sl, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl) prm%p = math_expand(prm%p, N_sl) prm%q = math_expand(prm%q, N_sl) prm%tau_0 = math_expand(prm%tau_0, N_sl) prm%B = math_expand(prm%B, N_sl) + prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl ! sanity checks if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' @@ -245,14 +245,15 @@ module function plastic_dislotwin_init() result(myPlasticity) if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0' if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' - if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' + if (any(prm%Q_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_sl' if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl' if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' + if (any(prm%d_caron < 0.0_pReal)) extmsg = trim(extmsg)//' d_caron(D_a,b_sl)' if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl' if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl' else slipActive rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray - allocate(prm%b_sl,prm%Q_s,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray) + allocate(prm%b_sl,prm%Q_sl,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray) allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0)) endif slipActive @@ -631,7 +632,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) integer :: i real(pReal) :: & f_unrotated, & - rho_dip_distance, & + d_hat, & v_cl, & !< climb velocity tau, & sigma_cl, & !< climb stress @@ -639,15 +640,14 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_rho_dip_formation, & dot_rho_dip_climb, & - rho_dip_distance_min, & dot_gamma_sl real(pReal), dimension(param(ph)%sum_N_tw) :: & dot_gamma_tw real(pReal), dimension(param(ph)%sum_N_tr) :: & dot_gamma_tr - associate(prm => param(ph), stt => state(ph), & - dot => dotState(ph), dst => dependentState(ph)) + + associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,en)) & @@ -656,8 +656,6 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) call kinetics_sl(Mp,T,ph,en,dot_gamma_sl) dot%gamma_sl(:,en) = abs(dot_gamma_sl) - rho_dip_distance_min = prm%D_a*prm%b_sl - slipState: do i = 1, prm%sum_N_sl tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) @@ -665,14 +663,14 @@ 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 - rho_dip_distance = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau)) - rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,en)) - rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i)) + d_hat = 3.0_pReal*prm%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)) - dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(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)) - if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then + if (dEq(d_hat,prm%d_caron(i))) then dot_rho_dip_climb(i) = 0.0_pReal else ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 @@ -685,17 +683,17 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & - / (rho_dip_distance-rho_dip_distance_min(i)) + / (d_hat-prm%d_caron(i)) endif endif significantSlipStress enddo slipState dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) & - dot_rho_dip_formation & - - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl) + - 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*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & + - 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) @@ -885,7 +883,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & significantStress: where(tau_eff > tol_math_check) stressRatio = tau_eff/prm%tau_0 StressRatio_p = stressRatio** prm%p - BoltzmannRatio = prm%Q_s/(kB*T) + BoltzmannRatio = prm%Q_sl/(kB*T) v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) v_run_inverse = prm%B/(tau_eff*prm%b_sl) From f96b0371ac27ee38331613cbbf3a4d18b80b3405 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 21:48:40 +0200 Subject: [PATCH 02/10] proper output name/format --- PRIVATE | 2 +- src/phase_mechanical_plastic_kinehardening.f90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/PRIVATE b/PRIVATE index 174ecac2d..347932db3 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 174ecac2d3ab7596bdb60184d6bb9e1a52cb7378 +Subproject commit 347932db368ec5c0d9e9092ef106c225605d830a diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 880b128a9..f1c9b0332 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -359,11 +359,11 @@ module subroutine plastic_kinehardening_results(ph,group) case ('xi') if(prm%sum_N_sl>0) call results_writeDataset(stt%xi,group,trim(prm%output(o)), & 'resistance against plastic slip','Pa') - case ('tau_b') !ToDo: chi + case ('chi') if(prm%sum_N_sl>0) call results_writeDataset(stt%chi,group,trim(prm%output(o)), & 'back stress','Pa') case ('sgn(gamma)') - if(prm%sum_N_sl>0) call results_writeDataset(stt%sgn_gamma,group,trim(prm%output(o)), & ! ToDo: could be int + if(prm%sum_N_sl>0) call results_writeDataset(int(stt%sgn_gamma),group,trim(prm%output(o)), & 'sense of shear','1') case ('chi_0') if(prm%sum_N_sl>0) call results_writeDataset(stt%chi_0,group,trim(prm%output(o)), & From 2d1d99542c4ab6e341c3b66bf8370e7224088f89 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 21:51:51 +0200 Subject: [PATCH 03/10] takeover from original paper --- .../config/Phase_DisloUCLA_Tungsten.config | 42 ------------------- examples/config/Phase_Dislotungsten_W.yaml | 26 ++++++++++++ 2 files changed, 26 insertions(+), 42 deletions(-) delete mode 100644 examples/config/Phase_DisloUCLA_Tungsten.config create mode 100644 examples/config/Phase_Dislotungsten_W.yaml diff --git a/examples/config/Phase_DisloUCLA_Tungsten.config b/examples/config/Phase_DisloUCLA_Tungsten.config deleted file mode 100644 index 1eef32577..000000000 --- a/examples/config/Phase_DisloUCLA_Tungsten.config +++ /dev/null @@ -1,42 +0,0 @@ -[Tungsten] -elasticity hooke -plasticity disloucla - -(output) edge_density -(output) dipole_density -(output) shear_rate_slip -(output) accumulated_shear_slip -(output) resolved_stress_slip -(output) threshold_stress_slip - -grainsize 2.7e-5 # Average grain size [m] 2.0e-5 -SolidSolutionStrength 0.0 # Strength due to elements in solid solution - -### Dislocation glide parameters ### -#per family -Nslip 12 -slipburgers 2.72e-10 # Burgers vector of slip system [m] -rhoedge0 1.0e12 # Initial edge dislocation density [m/m**3] -rhoedgedip0 1.0 # Initial edged dipole dislocation density [m/m**3] -Qedge 2.61154e-19 # Activation energy for dislocation glide [J], 1.63 eV -v0 1 # Initial glide velocity [m/s] -p_slip 0.86 # p-exponent in glide velocity -q_slip 1.69 # q-exponent in glide velocity -tau_peierls 2.03e9 # peierls stress [Pa] - -#mobility law -kink_height 2.567e-10 # kink height sqrt(6)/3*lattice_parameter [m] -omega 9.1e11 # attemp frequency (from kMC paper) [s^(-1)] -kink_width 29.95e-10 # kink pair width ~ 11 b (kMC paper) [m] -dislolength 78e-10 # dislocation length (ideally lambda) [m] initial value 11b -friction_coeff 8.3e-5 # friction coeff. B [Pa*s] - -#hardening -dipoleformationfactor 0 # to have hardening due to dipole formation off -CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path -D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s] -Qsd 4.5e-19 # Activation energy for climb [J] -Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b] -Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b] -interaction_slipslip 0.009 0.72 0.009 0.05 0.05 0.06 0.09 -nonschmid_coefficients 0.938 0.71 4.43 0.0 0.0 0.0 diff --git a/examples/config/Phase_Dislotungsten_W.yaml b/examples/config/Phase_Dislotungsten_W.yaml new file mode 100644 index 000000000..38cf2d4e4 --- /dev/null +++ b/examples/config/Phase_Dislotungsten_W.yaml @@ -0,0 +1,26 @@ +type: dislotungsten + +N_sl: [12] + +rho_mob_0: [1.0e+9] +rho_dip_0: [1.0] + +nu_a: [9.1e+11] +b_sl: [2.72e-10] +Delta_H_kp,0: [2.61154e-19] # 1.63 eV, Delta_H0 + +tau_Peierls: [2.03e+9] +p_sl: [0.86] +q_sl: [1.69] +h: [2.566e-10] +w: [2.992e-09] +B: [8.3e-5] +D_a: 1.0 # d_edge + +# climb (disabled) +D_0: 0.0 +Q_cl: 0.0 +V_cl: [0.0] + +h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09] +a_nonSchmid: [0.938, 0.71, 4.43] From c46813657f38e0a898bba78eb2d78fc71dfc1d60 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 21:56:34 +0200 Subject: [PATCH 04/10] better to understand --- ...phase_mechanical_plastic_dislotungsten.f90 | 33 +++++++++---------- 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 6ca8e440e..0fc0076fc 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -56,7 +56,7 @@ submodule(phase:plastic) dislotungsten type :: tDisloTungstendependentState real(pReal), dimension(:,:), allocatable :: & Lambda_sl, & - threshold_stress + tau_pass end type tDisloTungstendependentState !-------------------------------------------------------------------------------------------------- @@ -246,8 +246,8 @@ module function plastic_dislotungsten_init() result(myPlasticity) 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' - allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal) - allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers), source=0.0_pReal) + 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) end associate @@ -380,8 +380,8 @@ module subroutine dislotungsten_dependentState(ph,en) associate(prm => param(ph), stt => state(ph),dst => dependentState(ph)) dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) - dst%threshold_stress(:,en) = prm%mu*prm%b_sl & - * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) + dst%tau_pass(:,en) = prm%mu*prm%b_sl & + * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) @@ -416,7 +416,7 @@ module subroutine plastic_dislotungsten_results(ph,group) if(prm%sum_N_sl>0) call results_writeDataset(dst%Lambda_sl,group,trim(prm%output(o)), & 'mean free path for slip','m') case('tau_pass') - if(prm%sum_N_sl>0) call results_writeDataset(dst%threshold_stress,group,trim(prm%output(o)), & + if(prm%sum_N_sl>0) call results_writeDataset(dst%tau_pass,group,trim(prm%output(o)), & 'threshold stress for slip','Pa') end select enddo outputsLoop @@ -456,8 +456,7 @@ pure subroutine kinetics(Mp,T,ph,en, & StressRatio_p,StressRatio_pminus1, & dvel, vel, & tau_pos,tau_neg, & - t_n, t_k, dtk,dtn, & - needsGoodName ! ToDo: @Karo: any idea? + t_n, t_k, dtk,dtn integer :: j associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) @@ -475,13 +474,12 @@ pure subroutine kinetics(Mp,T,ph,en, & dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, & effectiveLength => dst%Lambda_sl(:,en) - prm%w) - significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check) - StressRatio = (abs(tau_pos)-dst%threshold_stress(:,en))/prm%tau_Peierls + significantPositiveTau: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) + StressRatio = (abs(tau_pos)-dst%tau_pass(:,en))/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) + t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) vel = prm%h/(t_n + t_k) @@ -492,7 +490,7 @@ pure subroutine kinetics(Mp,T,ph,en, & end where significantPositiveTau if (present(ddot_gamma_dtau_pos)) then - significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check) + significantPositiveTau2: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_pos @@ -505,13 +503,12 @@ pure subroutine kinetics(Mp,T,ph,en, & end where significantPositiveTau2 endif - significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check) - StressRatio = (abs(tau_neg)-dst%threshold_stress(:,en))/prm%tau_Peierls + significantNegativeTau: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) + StressRatio = (abs(tau_neg)-dst%tau_pass(:,en))/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) + t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) vel = prm%h/(t_n + t_k) @@ -522,7 +519,7 @@ pure subroutine kinetics(Mp,T,ph,en, & end where significantNegativeTau if (present(ddot_gamma_dtau_neg)) then - significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check) + significantNegativeTau2: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_neg From 2a78174547b473f1d96ff5916ad5e805c2f53c88 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 22:04:39 +0200 Subject: [PATCH 05/10] standard names in source code --- examples/config/Phase_Dislotungsten_W.yaml | 2 +- ...phase_mechanical_plastic_dislotungsten.f90 | 36 +++++++++---------- 2 files changed, 17 insertions(+), 21 deletions(-) diff --git a/examples/config/Phase_Dislotungsten_W.yaml b/examples/config/Phase_Dislotungsten_W.yaml index 38cf2d4e4..bf8796cfa 100644 --- a/examples/config/Phase_Dislotungsten_W.yaml +++ b/examples/config/Phase_Dislotungsten_W.yaml @@ -16,7 +16,7 @@ h: [2.566e-10] w: [2.992e-09] B: [8.3e-5] D_a: 1.0 # d_edge - + # climb (disabled) D_0: 0.0 Q_cl: 0.0 diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 0fc0076fc..e7cf32080 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -18,7 +18,7 @@ submodule(phase:plastic) dislotungsten Q_cl = 1.0_pReal !< activation energy for dislocation climb real(pReal), allocatable, dimension(:) :: & b_sl, & !< magnitude of Burgers vector [m] - D_a, & + d_caron, & !< distance of spontaneous annhihilation i_sl, & !< Adj. parameter for distance between 2 forest dislocations f_at, & !< factor to calculate atomic volume tau_Peierls, & !< Peierls stress @@ -172,7 +172,6 @@ module function plastic_dislotungsten_init() result(myPlasticity) prm%D_0 = pl%get_asFloat('D_0') prm%Q_cl = pl%get_asFloat('Q_cl') prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal - prm%D_a = pl%get_asFloat('D_a') * prm%b_sl prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.) @@ -191,7 +190,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) prm%B = math_expand(prm%B, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl) prm%f_at = math_expand(prm%f_at, N_sl) - prm%D_a = math_expand(prm%D_a, N_sl) + prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl ! sanity checks if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' @@ -202,12 +201,13 @@ module function plastic_dislotungsten_init() result(myPlasticity) if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' if (any(prm%tau_Peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls' - if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' D_a or b_sl' + if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' + if (any(prm%d_caron < 0.0_pReal)) extmsg = trim(extmsg)//' d_caron(D_a,b_sl)' if (any(prm%f_at <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl' else slipActive rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray - allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%f_at,prm%tau_Peierls, & + allocate(prm%b_sl,prm%d_caron,prm%i_sl,prm%f_at,prm%tau_Peierls, & prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, & source = emptyRealArray) allocate(prm%forestProjection(0,0)) @@ -316,8 +316,6 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) ph, & en - real(pReal) :: & - VacancyDiffusion real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_pos, dot_gamma_neg,& tau_pos,& @@ -325,38 +323,36 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) v_cl, & dot_rho_dip_formation, & dot_rho_dip_climb, & - dip_distance + d_hat - associate(prm => param(ph), stt => state(ph),& - dot => dotState(ph), dst => dependentState(ph)) + associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) call kinetics(Mp,T,ph,en,& dot_gamma_pos,dot_gamma_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) dot%gamma_sl(:,en) = (dot_gamma_pos+dot_gamma_neg) ! ToDo: needs to be abs - VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg dot_rho_dip_formation = 0.0_pReal dot_rho_dip_climb = 0.0_pReal else where - dip_distance = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & - prm%D_a, & ! lower limit - dst%Lambda_sl(:,en)) ! upper limit - dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation + d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & + prm%d_caron, & ! lower limit + dst%Lambda_sl(:,en)) ! upper limit + dot_rho_dip_formation = merge(2.0_pReal*d_hat* stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation 0.0_pReal, & prm%dipoleformation) - v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) & - * (1.0_pReal/(dip_distance+prm%D_a)) - dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? + 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)) & + * (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 dot%rho_mob(:,en) = abs(dot%gamma_sl(:,en))/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication - dot_rho_dip_formation & - - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en)) ! Spontaneous annihilation of 2 single edge dislocations + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en)) ! Spontaneous annihilation of 2 edges dot%rho_dip(:,en) = dot_rho_dip_formation & - - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,en)*abs(dot%gamma_sl(:,en)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*abs(dot%gamma_sl(:,en)) & ! Spontaneous annihilation of an edge with a dipole - dot_rho_dip_climb end associate From 527fd306e27b4f7731474ce9d945644f9a3ca2f2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 22:25:06 +0200 Subject: [PATCH 06/10] change of behavior - gamma (state) increases monotoneously - region of spontaneous annihilation is ignored in dipole formation --- src/phase_mechanical_plastic_dislotungsten.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index e7cf32080..4f86b23b1 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -331,16 +331,16 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot_gamma_pos,dot_gamma_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) - dot%gamma_sl(:,en) = (dot_gamma_pos+dot_gamma_neg) ! ToDo: needs to be abs + dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg) - where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg + where(dEq0(tau_pos)) ! ToDo: use avg of +/- 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)), & + d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & ! ToDo: use avg of +/- prm%d_caron, & ! lower limit dst%Lambda_sl(:,en)) ! upper limit - dot_rho_dip_formation = merge(2.0_pReal*d_hat* stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation + 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)) & @@ -348,11 +348,11 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) 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 - dot%rho_mob(:,en) = abs(dot%gamma_sl(:,en))/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication + dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication - dot_rho_dip_formation & - - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en)) ! Spontaneous annihilation of 2 edges + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges 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(:,en)) & ! Spontaneous annihilation of an edge with a dipole + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole - dot_rho_dip_climb end associate @@ -373,7 +373,7 @@ module subroutine dislotungsten_dependentState(ph,en) dislocationSpacing - associate(prm => param(ph), stt => state(ph),dst => dependentState(ph)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) dst%tau_pass(:,en) = prm%mu*prm%b_sl & From d4ffc778c2d491f7930c488197bda23f56c2d071 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 24 Jul 2021 07:09:35 +0200 Subject: [PATCH 07/10] easier to read removed comment regarding use of dot_state in kinetics_t(w/r). Data stored in dotState is not reliable, FPI integrator for writes to it and Runge-Kutta calls the dot state function at different time steps --- src/phase_mechanical_plastic_dislotwin.f90 | 243 +++++++++++---------- 1 file changed, 123 insertions(+), 120 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 4aa62bcb8..1d6fadd13 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -517,7 +517,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) integer :: i,k,l,m,n real(pReal) :: & f_unrotated,StressRatio_p,& - BoltzmannRatio, & + E_kB_T, & ddot_gamma_dtau, & tau real(pReal), dimension(param(ph)%sum_N_sl) :: & @@ -587,7 +587,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) shearBandingContribution: if(dNeq0(prm%v_sb)) then - BoltzmannRatio = prm%E_sb/(kB*T) + E_kB_T = prm%E_sb/(kB*T) call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? do i = 1,6 @@ -597,8 +597,8 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) significantShearBandStress: if (abs(tau) > tol_math_check) then StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb - dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) - ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb & + dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau) + ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb & * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) @@ -649,58 +649,58 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_unrotated = 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) + call kinetics_sl(Mp,T,ph,en,dot_gamma_sl) + dot%gamma_sl(:,en) = abs(dot_gamma_sl) - slipState: do i = 1, prm%sum_N_sl - tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) + slipState: do i = 1, prm%sum_N_sl + tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) - significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then - 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 = math_clip(d_hat, right = dst%Lambda_sl(i,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)) - - if (dEq(d_hat,prm%d_caron(i))) then + significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then + dot_rho_dip_formation(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal - 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)), & - 1.0_pReal, & - prm%ExtendedDislocations) - v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & - * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) + else significantSlipStress + d_hat = 3.0_pReal*prm%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)) - dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & - / (d_hat-prm%d_caron(i)) - endif - endif significantSlipStress - enddo slipState + 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)) - 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) + if (dEq(d_hat,prm%d_caron(i))) then + dot_rho_dip_climb(i) = 0.0_pReal + 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)), & + 1.0_pReal, & + prm%ExtendedDislocations) + v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & + * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) - 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_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & + / (d_hat-prm%d_caron(i)) + endif + endif significantSlipStress + enddo slipState - call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) - dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char + dot%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) - call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) - dot%f_tr(:,en) = f_unrotated*dot_gamma_tr + 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 + + call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) + dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char + + call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) + dot%f_tr(:,en) = f_unrotated*dot_gamma_tr end associate @@ -865,7 +865,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & tau, & stressRatio, & StressRatio_p, & - BoltzmannRatio, & + Q_kB_T, & v_wait_inverse, & !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned) v_run_inverse, & !< inverse of the velocity of a free moving dislocation (unsigned) dV_wait_inverse_dTau, & @@ -874,33 +874,34 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & tau_eff !< effective resolved stress integer :: i + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)] + tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)] - tau_eff = abs(tau)-dst%tau_pass(:,en) + tau_eff = abs(tau)-dst%tau_pass(:,en) - significantStress: where(tau_eff > tol_math_check) - stressRatio = tau_eff/prm%tau_0 - StressRatio_p = stressRatio** prm%p - BoltzmannRatio = prm%Q_sl/(kB*T) - v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) - v_run_inverse = prm%B/(tau_eff*prm%b_sl) + significantStress: where(tau_eff > tol_math_check) + stressRatio = tau_eff/prm%tau_0 + StressRatio_p = stressRatio** prm%p + Q_kB_T = prm%Q_sl/(kB*T) + v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q) + v_run_inverse = prm%B/(tau_eff*prm%b_sl) - dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) + dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) - dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio & - * (stressRatio**(prm%p-1.0_pReal)) & - * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & - / prm%tau_0 - dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff - dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) & - / (v_wait_inverse+v_run_inverse)**2.0_pReal - ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl - else where significantStress - dot_gamma_sl = 0.0_pReal - ddot_gamma_dtau = 0.0_pReal - end where significantStress + dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * Q_kB_T & + * (stressRatio**(prm%p-1.0_pReal)) & + * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & + / prm%tau_0 + dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff + dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) & + / (v_wait_inverse+v_run_inverse)**2.0_pReal + ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl + else where significantStress + dot_gamma_sl = 0.0_pReal + ddot_gamma_dtau = 0.0_pReal + end where significantStress end associate @@ -943,34 +944,35 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& integer :: i,s1,s2 + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - do i = 1, prm%sum_N_tw - tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) - isFCC: if (prm%fccTwinTransNucleation) then - s1=prm%fcc_twinNucleationSlipPair(1,i) - s2=prm%fcc_twinNucleationSlipPair(2,i) - if (tau(i) < dst%tau_r_tw(i,en)) then ! ToDo: correct? - Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& - abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& ! ToDo: MD: it would be more consistent to use shearrates from state - (prm%L_tw*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i)))) ! P_ncs - else - Ndot0=0.0_pReal - end if - else isFCC - Ndot0=prm%dot_N_0_tw(i) - endif isFCC - enddo + do i = 1, prm%sum_N_tw + tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) + isFCC: if (prm%fccTwinTransNucleation) then + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau(i) < dst%tau_r_tw(i,en)) then ! ToDo: correct? + Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& + abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& + (prm%L_tw*prm%b_sl(i))*& + (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i)))) + else + Ndot0=0.0_pReal + end if + else isFCC + Ndot0=prm%dot_N_0_tw(i) + endif isFCC + enddo - significantStress: where(tau > tol_math_check) - StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r - dot_gamma_tw = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r) - ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r - else where significantStress - dot_gamma_tw = 0.0_pReal - ddot_gamma_dtau = 0.0_pReal - end where significantStress + significantStress: where(tau > tol_math_check) + StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r + dot_gamma_tw = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r) + ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r + else where significantStress + dot_gamma_tw = 0.0_pReal + ddot_gamma_dtau = 0.0_pReal + end where significantStress end associate @@ -1009,36 +1011,37 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& Ndot0, & stressRatio_s, & ddot_gamma_dtau - integer :: i,s1,s2 + + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - do i = 1, prm%sum_N_tr - tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i)) - isFCC: if (prm%fccTwinTransNucleation) then - s1=prm%fcc_twinNucleationSlipPair(1,i) - s2=prm%fcc_twinNucleationSlipPair(2,i) - if (tau(i) < dst%tau_r_tr(i,en)) then ! ToDo: correct? - Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& - abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& ! ToDo: MD: it would be more consistent to use shearrates from state - (prm%L_tr*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i)))) ! P_ncs - else - Ndot0=0.0_pReal - end if - else isFCC - Ndot0=prm%dot_N_0_tr(i) - endif isFCC - enddo + do i = 1, prm%sum_N_tr + tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i)) + isFCC: if (prm%fccTwinTransNucleation) then + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau(i) < dst%tau_r_tr(i,en)) then ! ToDo: correct? + Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& + abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& + (prm%L_tr*prm%b_sl(i))*& + (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i)))) + else + Ndot0=0.0_pReal + end if + else isFCC + Ndot0=prm%dot_N_0_tr(i) + endif isFCC + enddo - significantStress: where(tau > tol_math_check) - StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s - dot_gamma_tr = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s) - ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s - else where significantStress - dot_gamma_tr = 0.0_pReal - ddot_gamma_dtau = 0.0_pReal - end where significantStress + significantStress: where(tau > tol_math_check) + StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s + dot_gamma_tr = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s) + ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s + else where significantStress + dot_gamma_tr = 0.0_pReal + ddot_gamma_dtau = 0.0_pReal + end where significantStress end associate From 6f906647eb1198d6563851089f274d01b73a2daa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 24 Jul 2021 09:10:16 +0200 Subject: [PATCH 08/10] documenting hex --- PRIVATE | 2 +- .../mechanical/plastic/phenopowerlaw_Mg.yaml} | 33 ++++++++++++------- .../mechanical/plastic/phenopowerlaw_Ti.yaml | 5 ++- 3 files changed, 26 insertions(+), 14 deletions(-) rename examples/config/{Phase_Phenopowerlaw_Magnesium.yaml => phase/mechanical/plastic/phenopowerlaw_Mg.yaml} (65%) diff --git a/PRIVATE b/PRIVATE index 347932db3..020cf0621 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 347932db368ec5c0d9e9092ef106c225605d830a +Subproject commit 020cf06218fc4f14cb55ae0527972d12a57343bc diff --git a/examples/config/Phase_Phenopowerlaw_Magnesium.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml similarity index 65% rename from examples/config/Phase_Phenopowerlaw_Magnesium.yaml rename to examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml index 4d4b9eb8e..b29075173 100644 --- a/examples/config/Phase_Phenopowerlaw_Magnesium.yaml +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml @@ -1,20 +1,29 @@ -N_sl: [3, 3, 0, 6, 0, 6] -N_tw: [6, 0, 0, 6] -h_0_tw-tw: 50.0e+6 -h_0_sl-sl: 500.0e+6 -h_0_tw-sl: 150.0e+6 -h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] -h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] -h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] -h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] -output: [xi_sl, xi_tw] type: phenopowerlaw -xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6] +references: + - F. Wang et al., + Acta Materialia 80:77-93, 2014, + https://doi.org/10.1016/j.actamat.2014.07.048 + +output: [xi_sl, xi_tw] + +N_sl: [3, 3, 0, 6, 0, 6] # basal, 1. prism, -, 1. pyr, -, 2. pyr +N_tw: [6, 0, 6] # tension, -, compression + +xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6] xi_inf_sl: [40.e+6, 135.e+6, 0., 150.e+6, 0., 150.e+6] -xi_0_tw: [40.e+6, 0., 0., 60.e+6] +xi_0_tw: [40.e+6, 0., 60.e+6] + a_sl: 2.25 dot_gamma_0_sl: 0.001 dot_gamma_0_tw: 0.001 n_sl: 20 n_tw: 20 f_sat_sl-tw: 10.0 + +h_0_sl-sl: 500.0e+6 +h_0_tw-tw: 50.0e+6 +h_0_tw-sl: 150.0e+6 +h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] +h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] +h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] +h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml index 5e05206df..804cf4541 100644 --- a/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml @@ -7,14 +7,17 @@ references: Acta Materialia 132:598-610, 2017, https://doi.org/10.1016/j.actamat.2017.05.015 output: [gamma_sl] -N_sl: [3, 3, 0, 0, 12] + +N_sl: [3, 3, 0, 0, 12] # basal, 1. prism, -, -, 2. pyr n_sl: 20 a_sl: 2.0 dot_gamma_0_sl: 0.001 h_0_sl-sl: 200.e+6 + # C. Zambaldi et al.: xi_0_sl: [349.e+6, 150.e+6, 0.0, 0.0, 1107.e+6] xi_inf_sl: [568.e+6, 150.e+7, 0.0, 0.0, 3420.e+6] # L. Wang et al. : # xi_0_sl: [127.e+6, 96.e+6, 0.0, 0.0, 240.e+6] + h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] From 931dc9955752ec99759c2aa36f9c137ee5ca0b89 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 24 Jul 2021 11:21:07 +0200 Subject: [PATCH 09/10] not used climb formulation was updated a while ago --- PRIVATE | 2 +- .../config/phase/mechanical/plastic/dislotwin_IF-steel.yaml | 1 - src/phase_mechanical_plastic_dislotwin.f90 | 5 +---- 3 files changed, 2 insertions(+), 6 deletions(-) diff --git a/PRIVATE b/PRIVATE index 020cf0621..c656b6f08 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 020cf06218fc4f14cb55ae0527972d12a57343bc +Subproject commit c656b6f08489756c9ee6a6e1a62858c8b7836f10 diff --git a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml index e48d5437f..362797d0b 100644 --- a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml +++ b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml @@ -19,6 +19,5 @@ i_sl: [23.3, 23.3] D_a: 7.4 # C_anni B: [0.001, 0.001] h_sl-sl: [0.1, 0.72, 0.1, 0.053, 0.053, 0.073, 0.137, 0.72, 0.72, 0.053, 0.053, 0.053, 0.053, 0.073, 0.073, 0.073, 0.073, 0.073, 0.073, 0.137, 0.073, 0.073, 0.137, 0.073] -D_0: 4.0e-05 Q_cl: 5.4e-19 # no recovery! D: 40.e-6 # estimated diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 1d6fadd13..2cc6f4af8 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -16,7 +16,6 @@ submodule(phase:plastic) dislotwin real(pReal) :: & mu = 1.0_pReal, & !< equivalent shear modulus nu = 1.0_pReal, & !< equivalent shear Poisson's ratio - D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient 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 @@ -214,7 +213,6 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%B = pl%get_as1dFloat('B', requiredSize=size(N_sl), & defaultVal=[(0.0_pReal, i=1,size(N_sl))]) - prm%D_0 = pl%get_asFloat('D_0') prm%Q_cl = pl%get_asFloat('Q_cl') prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.) @@ -230,7 +228,7 @@ module function plastic_dislotwin_init() result(myPlasticity) rho_dip_0 = math_expand(rho_dip_0, N_sl) prm%v_0 = math_expand(prm%v_0, N_sl) prm%b_sl = math_expand(prm%b_sl, N_sl) - prm%Q_sl = math_expand(prm%Q_sl, N_sl) + prm%Q_sl = math_expand(prm%Q_sl, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl) prm%p = math_expand(prm%p, N_sl) prm%q = math_expand(prm%q, N_sl) @@ -239,7 +237,6 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl ! sanity checks - if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' From b3f5e1223235043c9c736130058b1b564d19f194 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 24 Jul 2021 11:43:14 +0200 Subject: [PATCH 10/10] using burgers vector of tw/tr system seems to make more sense here otherwise twinning will not work in many cases. Matching number is only required for nucleation of tw/tr --- PRIVATE | 2 +- src/phase_mechanical_plastic_dislotwin.f90 | 16 +++++++--------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/PRIVATE b/PRIVATE index c656b6f08..5049664c5 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit c656b6f08489756c9ee6a6e1a62858c8b7836f10 +Subproject commit 5049664c5cac3cc1571c7b61f3345f1ba8d627f6 diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 2cc6f4af8..4242d0e7a 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -205,7 +205,7 @@ module function plastic_dislotwin_init() result(myPlasticity) rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl)) prm%v_0 = pl%get_as1dFloat('v_0', requiredSize=size(N_sl)) prm%b_sl = pl%get_as1dFloat('b_sl', requiredSize=size(N_sl)) - prm%Q_sl = pl%get_as1dFloat('Q_sl', requiredSize=size(N_sl)) + prm%Q_sl = pl%get_as1dFloat('Q_sl', requiredSize=size(N_sl)) prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(N_sl)) prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl)) prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl)) @@ -758,19 +758,17 @@ module subroutine dislotwin_dependentState(T,ph,en) dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) !* threshold stress for growing twin/martensite - if(prm%sum_N_tw == prm%sum_N_sl) & - 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_sl) ! slip Burgers here correct? - if(prm%sum_N_tr == prm%sum_N_sl) & - 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_sl) & ! slip Burgers here correct? - + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr) + 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) + 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) & + + 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 and is the same for twin and trans + 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 = 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