all parameters as slip family-dependent

This commit is contained in:
Philip Eisenlohr 2023-08-04 12:43:56 -04:00
parent a2c9cfbbea
commit ffdf478b62
13 changed files with 131 additions and 114 deletions

@ -1 +1 @@
Subproject commit 8c05965ef4437598898f467a213ffa88938e860a
Subproject commit 4e195733964bfa36de70abb768a65a7216180d0f

View File

@ -10,10 +10,10 @@ output: [xi_sl, gamma_sl]
N_sl: [12]
n_sl: 20
a_sl: 3.7
h_0_sl-sl: 1.02e+9
dot_gamma_0_sl: [0.001]
n_sl: [20]
a_sl: [3.7]
xi_0_sl: [76.e+6]
xi_inf_sl: [266.e+6]
h_0_sl-sl: [1.02e+9]
h_sl-sl: [1, 1, 5.123, 0.574, 1.123, 1.123, 1]
dot_gamma_0_sl: 0.001

View File

@ -10,10 +10,10 @@ output: [xi_sl, gamma_sl]
N_sl: [12]
n_sl: 20
a_sl: 5.4
h_0_sl-sl: 281.5e+6
dot_gamma_0_sl: [7.5e-5]
n_sl: [20]
a_sl: [5.4]
xi_0_sl: [2.69e+6]
xi_inf_sl: [67.5e+6]
h_0_sl-sl: [0.2815e+9]
h_sl-sl: [1, 1, 5.123, 0.574, 1.123, 1.123, 1]
dot_gamma_0_sl: 7.5e-5

View File

@ -15,10 +15,10 @@ output: [xi_sl, gamma_sl]
N_sl: [12]
n_sl: 83.3
a_sl: 1.0
h_0_sl-sl: 75.0e+6
dot_gamma_0_sl: [0.001]
n_sl: [83.3]
a_sl: [1.0]
xi_0_sl: [26.25e+6]
xi_inf_sl: [53.0e+6]
h_0_sl-sl: [75.0e+6]
h_sl-sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4]
dot_gamma_0_sl: 0.001

View File

@ -10,10 +10,10 @@ output: [xi_sl, gamma_sl]
N_sl: [12]
n_sl: 20
a_sl: 0.6
h_0_sl-sl: 3.5e+8
dot_gamma_0_sl: [3.e-3]
n_sl: [20]
a_sl: [0.6]
xi_0_sl: [1.6e+6]
xi_inf_sl: [96.4e+6]
h_0_sl-sl: [0.35e+9]
h_sl-sl: [1, 1, 5.123, 0.574, 1.123, 1.123, 1]
dot_gamma_0_sl: 3.e-3

View File

@ -12,10 +12,10 @@ output: [xi_sl, gamma_sl]
N_sl: [12, 12]
n_sl: 20
a_sl: 2.25
h_0_sl-sl: 1.0e+9
dot_gamma_0_sl: [0.001, 0.001]
n_sl: [20, 20]
a_sl: [2.25, 2.25]
xi_0_sl: [95.e+6, 96.e+6]
xi_inf_sl: [222.e+6, 412.e+6]
h_0_sl-sl: [1.0e+9, 1.0e+9]
h_sl-sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4]
dot_gamma_0_sl: 0.001

View File

@ -14,16 +14,16 @@ xi_0_sl: [10.e+6, 55.e+6, 60.e+6, 0., 60.e+6]
xi_inf_sl: [40.e+6, 135.e+6, 150.e+6, 0., 150.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
a_sl: [2.25, 2.25, 2.25, 1, 2.25]
dot_gamma_0_sl: [0.001, 0.001, 0.001, 0, 0.001]
dot_gamma_0_tw: [0.001, 0, 0.001]
n_sl: [20, 20, 20, 1, 20]
n_tw: [20, 1, 20]
f_sat_sl-tw: [10.0, 10.0, 10.0, 0, 10.0]
h_0_sl-sl: 500.0e+6
h_0_tw-tw: 50.0e+6
h_0_tw-sl: 150.0e+6
h_0_sl-sl: [0.5e+9, 0.5e+9, 0.5e+9, 0, 0.5e+9]
h_0_tw-tw: [50.0e+6, 0, 50.0e+6]
h_0_tw-sl: [0.15e+9, 0, 0.15e+9]
h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,

View File

@ -10,10 +10,10 @@ output: [xi_sl, gamma_sl]
N_sl: [12]
n_sl: 20
a_sl: 0.9
h_0_sl-sl: 781.2e+6
xi_0_sl: [114.e+6]
xi_inf_sl: [207.e+6]
dot_gamma_0_sl: [0.001]
n_sl: [20]
a_sl: [0.9]
xi_0_sl: [0.114e+9]
xi_inf_sl: [0.207e+9]
h_0_sl-sl: [0.7812e+9]
h_sl-sl: [1, 1, 5.123, 0.574, 1.123, 1.123, 1]
dot_gamma_0_sl: 0.001

View File

@ -9,9 +9,9 @@ output: [xi_sl, gamma_sl]
N_sl: [2, 2, 2, 4, 2, 4, 2, 2, 4, 0, 0, 8]
n_sl: 6.0
a_sl: 2.0
h_0_sl-sl: 20.0e+6
n_sl: [6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 1, 1, 6.0]
a_sl: [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 1.0, 1.0, 2.0]
h_0_sl-sl: [20e+6, 20e+6, 20e+6, 20e+6, 20e+6, 20e+6, 20e+6, 20e+6, 20e+6, 0.0, 0.0, 20e+6]
xi_0_sl: [8.5e+6, 4.3e+6, 10.4e+6, 4.5e+6, 5.6e+6, 5.1e+6, 7.4e+6, 15.0e+6, 6.6e+6, 0.0, 0.0, 12.0e+6]
xi_inf_sl: [11.0e+6, 9.0e+6, 11.0e+6, 9.0e+6, 10.0e+6, 10.0e+6, 10.0e+6, 10.0e+6, 9.0e+6, 0.0, 0.0, 13.0e+6]
h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
@ -30,4 +30,4 @@ h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
-1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, # 150
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0] # unused entries are indicated by -1.0
dot_gamma_0_sl: 2.6e-8
dot_gamma_0_sl: [2.6e-8, 2.6e-8, 2.6e-8, 2.6e-8, 2.6e-8, 2.6e-8, 2.6e-8, 2.6e-8, 2.6e-8, 1.0, 1.0, 2.6e-8]

View File

@ -12,17 +12,17 @@ output: [gamma_sl]
N_sl: [3, 3, 0, 12] # basal, prism, -, 1. pyr<c+a>
n_sl: 20
a_sl: 2.0
dot_gamma_0_sl: 0.001
h_0_sl-sl: 200.e+6
dot_gamma_0_sl: [0.001, 0.001, 0.0, 0.001]
n_sl: [20, 20, 1, 20]
a_sl: [2.0, 2.0, 1.0, 2.0]
# C. Zambaldi et al.:
xi_0_sl: [349.e+6, 150.e+6, 0.0, 1107.e+6]
xi_inf_sl: [568.e+6, 150.e+7, 0.0, 3420.e+6]
xi_0_sl: [0.349e+9, 0.15e+9, 0.0, 1.107e+9]
xi_inf_sl: [0.568e+9, 1.50e+9, 0.0, 3.420e+9]
# L. Wang et al. :
# xi_0_sl: [127.e+6, 96.e+6, 0.0, 240.e+6]
h_0_sl-sl: [0.2e+9, 0.2e+9, 0.0, 0.2e+9]
h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,

View File

@ -1,12 +1,14 @@
# Tasan et.al. 2015 Acta Materalia
# Tasan et.al. 2015 International Journal of Plasticity
# Diehl et.al. 2015 Meccanica
N_sl: [12, 12]
a_sl: 2.0
dot_gamma_0_sl: 0.001
h_0_sl-sl: 563.0e+9
h_sl-sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
type: phenopowerlaw
N_sl: [12, 12]
dot_gamma_0_sl: [0.001, 0.001]
n_sl: [20, 20]
a_sl: [2.0, 2.0]
xi_0_sl: [405.8e+6, 456.7e+6]
xi_inf_sl: [872.9e+6, 971.2e+6]
h_0_sl-sl: [563.0e+9, 563.0e+9]
h_sl-sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4]

View File

@ -13,12 +13,12 @@ phase:
plastic:
type: phenopowerlaw
N_sl: [12]
a_sl: 2.25
a_sl: [2.25]
atol_xi: 1.0
dot_gamma_0_sl: 0.001
h_0_sl-sl: 75.e+6
dot_gamma_0_sl: [0.001]
h_0_sl-sl: [75.e+6]
h_sl-sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
n_sl: [20]
output: [xi_sl]
xi_0_sl: [31.e+6]
xi_inf_sl: [63.e+6]

View File

@ -7,23 +7,21 @@
submodule(phase:plastic) phenopowerlaw
type :: tParameters
real(pREAL) :: &
dot_gamma_0_sl = 1.0_pREAL, & !< reference shear strain rate for slip
dot_gamma_0_tw = 1.0_pREAL, & !< reference shear strain rate for twin
n_sl = 1.0_pREAL, & !< stress exponent for slip
n_tw = 1.0_pREAL, & !< stress exponent for twin
f_sat_sl_tw = 1.0_pREAL, & !< push-up factor for slip saturation due to twinning
c_1 = 1.0_pREAL, &
c_2 = 1.0_pREAL, &
c_3 = 1.0_pREAL, &
c_4 = 1.0_pREAL, &
h_0_sl_sl = 1.0_pREAL, & !< reference hardening slip - slip
h_0_tw_sl = 1.0_pREAL, & !< reference hardening twin - slip
h_0_tw_tw = 1.0_pREAL, & !< reference hardening twin - twin
a_sl = 1.0_pREAL
real(pREAL), allocatable, dimension(:) :: &
dot_gamma_0_sl, & !< reference shear strain rate for slip
dot_gamma_0_tw, & !< reference shear strain rate for twin
a_sl, &
n_sl, & !< stress exponent for slip
n_tw, & !< stress exponent for twin
xi_inf_sl, & !< maximum critical shear stress for slip
h_int, & !< per family hardening activity (optional)
f_sat_sl_tw, & !< push-up factor for slip saturation due to twinning
c_1, &
c_2, &
c_3, &
c_4, &
h_0_sl_sl, & !< reference hardening slip - slip
h_0_tw_sl, & !< reference hardening twin - slip
h_0_tw_tw, & !< reference hardening twin - twin
gamma_char !< characteristic shear for twins
real(pREAL), allocatable, dimension(:,:) :: &
h_sl_sl, & !< slip resistance from slip activity
@ -81,6 +79,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
ph, i, &
N_fam, &
Nmembers, &
sizeState, sizeDotState, &
startIndex, endIndex
@ -90,7 +89,8 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
real(pREAL), dimension(:), allocatable :: &
xi_0_sl, & !< initial critical shear stress for slip
xi_0_tw, & !< initial critical shear stress for twin
a !< non-Schmid coefficients
a, & !< non-Schmid coefficients
ones
character(len=:), allocatable :: &
refs, &
extmsg
@ -134,20 +134,26 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
prm%output = pl%get_as1dStr('output',defaultVal=emptyStrArray)
#endif
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
N_tw = pl%get_as1dInt('N_tw',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl))
prm%sum_N_tw = sum(abs(N_tw))
!--------------------------------------------------------------------------------------------------
! slip related parameters
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then
prm%dot_gamma_0_sl = pl%get_asReal('dot_gamma_0_sl')
prm%n_sl = pl%get_asReal('n_sl')
prm%a_sl = pl%get_asReal('a_sl')
prm%h_0_sl_sl = pl%get_asReal('h_0_sl-sl')
xi_0_sl = math_expand(pl%get_as1dReal('xi_0_sl', requiredSize=size(N_sl)),N_sl)
prm%xi_inf_sl = math_expand(pl%get_as1dReal('xi_inf_sl', requiredSize=size(N_sl)),N_sl)
prm%h_int = math_expand(pl%get_as1dReal('h_int', requiredSize=size(N_sl), &
defaultVal=[(0.0_pREAL,i=1,size(N_sl))]),N_sl)
N_fam = size(N_sl)
ones = [(1.0_pREAL,i=1,N_fam)]
prm%dot_gamma_0_sl = math_expand(pl%get_as1dReal('dot_gamma_0_sl',requiredSize=N_fam,defaultVal= ones), N_sl)
prm%n_sl = math_expand(pl%get_as1dReal('n_sl', requiredSize=N_fam,defaultVal= ones), N_sl)
prm%a_sl = math_expand(pl%get_as1dReal('a_sl', requiredSize=N_fam,defaultVal= ones), N_sl)
prm%h_0_sl_sl = math_expand(pl%get_as1dReal('h_0_sl-sl', requiredSize=N_fam,defaultVal=0.0_pREAL*ones), N_sl)
prm%c_1 = math_expand(pl%get_as1dReal('c_1', requiredSize=N_fam,defaultVal=0.0_pREAL*ones), N_sl)
prm%c_2 = math_expand(pl%get_as1dReal('c_2', requiredSize=N_fam,defaultVal= ones), N_sl)
xi_0_sl = math_expand(pl%get_as1dReal('xi_0_sl', requiredSize=N_fam), N_sl)
prm%xi_inf_sl = math_expand(pl%get_as1dReal('xi_inf_sl', requiredSize=N_fam), N_sl)
prm%f_sat_sl_tw = math_expand(pl%get_as1dReal('f_sat_sl-tw',requiredSize=N_fam,defaultVal=0.0_pREAL*ones), N_sl)
prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph))
@ -166,35 +172,37 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
prm%systems_sl = crystal_labels_slip(N_sl,phase_lattice(ph))
! sanity checks
if ( prm%dot_gamma_0_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_sl'
if ( prm%a_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' a_sl'
if ( prm%n_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' n_sl'
if (any(prm%dot_gamma_0_sl <= 0.0_pREAL)) extmsg = trim(extmsg)//' dot_gamma_0_sl'
if (any(prm%a_sl <= 0.0_pREAL)) extmsg = trim(extmsg)//' a_sl'
if (any(prm%n_sl <= 0.0_pREAL)) extmsg = trim(extmsg)//' n_sl'
if (any(xi_0_sl <= 0.0_pREAL)) extmsg = trim(extmsg)//' xi_0_sl'
if (any(prm%xi_inf_sl <= 0.0_pREAL)) extmsg = trim(extmsg)//' xi_inf_sl'
else slipActive
xi_0_sl = emptyRealArray
allocate(prm%xi_inf_sl, &
prm%h_int, &
allocate(prm%dot_gamma_0_sl, & !< reference shear strain rate for slip
prm%a_sl, &
prm%n_sl, & !< stress exponent for slip
prm%xi_inf_sl, & !< maximum critical shear stress for slip
prm%f_sat_sl_tw, & !< push-up factor for slip saturation due to twinning
prm%c_1, &
prm%c_2, &
prm%h_0_sl_sl, & !< reference hardening slip - slip
source=emptyRealArray)
allocate(prm%h_sl_sl(0,0))
end if slipActive
!--------------------------------------------------------------------------------------------------
! twin related parameters
N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(abs(N_tw))
twinActive: if (prm%sum_N_tw > 0) then
prm%c_1 = pl%get_asReal('c_1',defaultVal=0.0_pREAL)
prm%c_2 = pl%get_asReal('c_2',defaultVal=1.0_pREAL)
prm%c_3 = pl%get_asReal('c_3',defaultVal=0.0_pREAL)
prm%c_4 = pl%get_asReal('c_4',defaultVal=0.0_pREAL)
prm%dot_gamma_0_tw = pl%get_asReal('dot_gamma_0_tw')
prm%n_tw = pl%get_asReal('n_tw')
prm%f_sat_sl_tw = pl%get_asReal('f_sat_sl-tw')
prm%h_0_tw_tw = pl%get_asReal('h_0_tw-tw')
xi_0_tw = math_expand(pl%get_as1dReal('xi_0_tw',requiredSize=size(N_tw)),N_tw)
N_fam = size(N_tw)
ones = [(1.0_pREAL,i=1,N_fam)]
prm%dot_gamma_0_tw = math_expand(pl%get_as1dReal('dot_gamma_0_tw', requiredSize=N_fam,defaultVal=ones), N_tw)
prm%n_tw = math_expand(pl%get_as1dReal('n_tw', requiredSize=N_fam,defaultVal= ones), N_tw)
prm%c_3 = math_expand(pl%get_as1dReal('c_3', requiredSize=N_fam,defaultVal=0.0_pREAL*ones), N_tw)
prm%c_4 = math_expand(pl%get_as1dReal('c_4', requiredSize=N_fam,defaultVal=0.0_pREAL*ones), N_tw)
prm%h_0_tw_tw = math_expand(pl%get_as1dReal('h_0_tw-tw', requiredSize=N_fam,defaultVal=0.0_pReal*ones), N_tw)
xi_0_tw = math_expand(pl%get_as1dReal('xi_0_tw', requiredSize=N_fam), N_tw)
prm%gamma_char = crystal_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%h_tw_tw = crystal_interaction_TwinByTwin(N_tw,pl%get_as1dReal('h_tw-tw'),phase_lattice(ph))
@ -203,27 +211,33 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
prm%systems_tw = crystal_labels_twin(N_tw,phase_lattice(ph))
! sanity checks
if (prm%dot_gamma_0_tw <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_tw'
if (prm%n_tw <= 0.0_pREAL) extmsg = trim(extmsg)//' n_tw'
if (any(prm%dot_gamma_0_tw <= 0.0_pREAL)) extmsg = trim(extmsg)//' dot_gamma_0_tw'
if (any(prm%n_tw <= 0.0_pREAL)) extmsg = trim(extmsg)//' n_tw'
if (any(xi_0_tw <= 0.0_pREAL)) extmsg = trim(extmsg)//' xi_0_tw'
else twinActive
xi_0_tw = emptyRealArray
allocate(prm%gamma_char,source=emptyRealArray)
allocate(prm%dot_gamma_0_tw, & !< reference shear strain rate for twin
prm%n_tw, & !< stress exponent for twin
prm%c_3, &
prm%c_4, &
prm%gamma_char, & !< characteristic shear for twins
prm%h_0_tw_sl, & !< reference hardening twin - slip
prm%h_0_tw_tw, & !< reference hardening twin - twin
source=emptyRealArray)
allocate(prm%h_tw_tw(0,0))
end if twinActive
!--------------------------------------------------------------------------------------------------
! slip-twin related parameters
slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
prm%h_0_tw_sl = pl%get_asReal('h_0_tw-sl')
prm%h_sl_tw = crystal_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dReal('h_sl-tw'), &
phase_lattice(ph))
prm%h_tw_sl = crystal_interaction_TwinBySlip(N_tw,N_sl,pl%get_as1dReal('h_tw-sl'), &
phase_lattice(ph))
prm%h_0_tw_sl = math_expand(pl%get_as1dReal('h_0_tw-sl',requiredSize=N_fam,defaultVal=0.0_pReal*ones), N_tw)
prm%h_sl_tw = crystal_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dReal('h_sl-tw'),phase_lattice(ph))
prm%h_tw_sl = crystal_interaction_TwinBySlip(N_tw,N_sl,pl%get_as1dReal('h_tw-sl'),phase_lattice(ph))
else slipAndTwinActive
allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0
allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0
prm%h_0_tw_sl = 0.0_pREAL
prm%h_0_tw_sl = [(0.0_pREAL,i=1,size(N_tw))]
end if slipAndTwinActive
!--------------------------------------------------------------------------------------------------
@ -344,10 +358,10 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
dotState
real(pREAL) :: &
xi_sl_sat_offset,&
sumF
real(pREAL), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_sl_pos,dot_gamma_sl_neg, &
xi_sl_sat_offset, &
left_SlipSlip
associate(prm => param(ph), stt => state(ph), &
@ -362,10 +376,11 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
left_SlipSlip = sign(abs(1.0_pREAL-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, &
1.0_pREAL-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))
dot_xi_sl = prm%h_0_sl_sl * (1.0_pREAL + prm%c_1 * sumF**prm%c_2) * (1.0_pREAL + prm%h_int) &
dot_xi_sl = prm%h_0_sl_sl * (1.0_pREAL + prm%c_1 * sumF**prm%c_2) &
* left_SlipSlip * matmul(prm%h_sl_sl,dot_gamma_sl) &
+ matmul(prm%h_sl_tw,dot_gamma_tw)