example files for TWIP/TRIP

I'm not able to reproduce the results from https://doi.org/10.1016/j.actamat.2016.07.032,
but these parameters give results that are qualitatively ok

examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml

@@ -13,7 +13,7 @@ output: [rho_dip, rho_mob]
 N_sl: [12, 12]
 
 b_sl: [2.49e-10, 2.49e-10]
-rho_mob_0: [2.81e12, 2.8e+12]
+rho_mob_0: [2.81e+12, 2.8e+12]
 rho_dip_0: [1.0, 1.0] # not given
 v_0: [1.4e+3, 1.4e+3]
 Q_sl: [1.57e-19, 1.57e-19] # Delta_F

examples/config/phase/mechanical/plastic/dislotwin_TWIP-TRIP.yaml

@@ -0,0 +1,72 @@
+type: dislotwin
+
+references:
+  - S.L. Wong et al.,
+    Acta Materialia 118:140-151, 2016,
+    https://doi.org/10.1016/j.actamat.2016.07.032
+  - K. Sedighiani et al.,
+    Mechanics of Materials, 164:104117, 2022,
+    https://doi.org/10.1016/j.mechmat.2021.104117
+  - L. Kubin et al.,
+    Acta Materialia 56:6040-6049,
+    https://doi.org/10.1016/j.actamat.2008.08.012
+
+output: [rho_mob, rho_dip, gamma_sl, Lambda_sl, tau_pass, f_tw, Lambda_tw, f_tr]
+
+# Glide
+N_sl: [12]
+b_sl: [2.56e-10] # a/sqrt(2)
+Q_sl: [3.5e-19]
+p_sl: [0.325]
+q_sl: [1.55]
+B: [0.001]
+i_sl: [30.0]
+v_0: [1.4e+3]
+tau_0: [5.5e+8] # adjusted
+D_a: 2.0
+Q_cl: 3.0e-19
+
+rho_mob_0: [5.0e+10]
+rho_dip_0: [5.0e+10]
+
+h_sl-sl: [0.122, 0.122, 0.625, 0.07, 0.137, 0.137, 0.122]
+
+# Twin
+N_tw: [12]
+b_tw: [1.47e-10] # a/sqrt(6)
+L_tw: 1.91e-7 # 1300 *b_tw
+i_tw: 10.0
+t_tw: [5.0e-8]
+p_tw: [7] # A, adjusted
+
+h_tw-tw: [1.0, 1.0]
+h_sl-tw: [1.0, 1.0, 1.0]
+
+# Transformation
+N_tr: [12]
+b_tr: [1.47e-10] # a/sqrt(6)
+L_tr: 2.21e-7 # 1500 *b_tr
+i_tr: 10.0 # adjusted
+t_tr: [1.0e-7]
+p_tr: [4] # B, adjusted
+V_mol: 7.09e-6
+c/a_hP: 1.633
+
+Delta_G: 1.2055e+2
+Delta_G,T: 2.5515
+Delta_G,T^2: 1.4952e-3
+
+h _tr-tr: [1.0, 1.0]
+h_sl-tr: [1.5, 1.5, 1.5]
+
+# Twin & Transformation
+T_ref: 293.15
+Gamma_sf: 2.833e-2
+Gamma_sf,T: 1.214e-4
+Gamma_sf,T^2: 1.473e-7
+
+x_c: 1.0e-9
+V_cs: 1.67e-29
+
+# Slip & Twin & Transformation
+D: 5.0e-5

src/phase_mechanical_plastic_dislotwin.f90

@@ -27,7 +27,6 @@ submodule(phase:plastic) dislotwin
       gamma_0_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
-      a_cF       = 1.0_pReal, &
       cOverA_hP  = 1.0_pReal, &
       V_mol      = 1.0_pReal, &
       rho        = 1.0_pReal
@@ -55,8 +54,8 @@ submodule(phase:plastic) dislotwin
     real(pReal),               allocatable, dimension(:,:) :: &
       h_sl_sl, &                                                                                    !< components of slip-slip interaction matrix
       h_sl_tw, &                                                                                    !< components of slip-twin interaction matrix
-      h_tw_tw, &                                                                                    !< components of twin-twin interaction matrix
       h_sl_tr, &                                                                                    !< components of slip-trans interaction matrix
+      h_tw_tw, &                                                                                    !< components of twin-twin interaction matrix
       h_tr_tr, &                                                                                    !< components of trans-trans interaction matrix
       n0_sl, &                                                                                      !< slip system normal
       forestProjection
@@ -302,7 +301,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
       prm%i_tr       = pl%get_asFloat('i_tr')
       prm%Delta_G    = polynomial(pl%asDict(),'Delta_G','T')
       prm%L_tr       = pl%get_asFloat('L_tr')
-      a_cF           = pl%get_asFloat('a_cF')
+      a_cF           = prm%b_tr(1)*sqrt(6.0_pReal)                                                  ! b_tr is Shockley partial
       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