223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

works for no active slip systems

This commit is contained in:
Martin Diehl 2020-03-16 11:53:36 +01:00
parent 5d7ff888fc
commit f8049b85be
2 changed files with 101 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

111
src/constitutive_plastic_disloUCLA.f90
View File

@ -12,15 +12,12 @@ submodule(constitutive) plastic_disloUCLA
type :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: &
D, & !< grain size D = 1.0_pReal, & !< grain size
mu, & !< equivalent shear modulus mu = 1.0_pReal, & !< equivalent shear modulus
D_0, & !< prefactor for self-diffusion coefficient D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient
Q_cl !< activation energy for dislocation climb Q_cl = 1.0_pReal !< activation energy for dislocation climb
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
rho_mob_0, & !< initial dislocation density
rho_dip_0, & !< initial dipole density
b_sl, & !< magnitude of burgers vector [m] b_sl, & !< magnitude of burgers vector [m]
nonSchmidCoeff, &
D_a, & D_a, &
i_sl, & !< Adj. parameter for distance between 2 forest dislocations i_sl, & !< Adj. parameter for distance between 2 forest dislocations
atomicVolume, & atomicVolume, &
@ -36,7 +33,7 @@ submodule(constitutive) plastic_disloUCLA
omega !< attempt frequency for kink pair nucleation omega !< attempt frequency for kink pair nucleation
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
h_sl_sl, & !< slip resistance from slip activity h_sl_sl, & !< slip resistance from slip activity
forestProjectionEdge forestProjection
real(pReal), allocatable, dimension(:,:,:) :: & real(pReal), allocatable, dimension(:,:,:) :: &
Schmid, & Schmid, &
nonSchmid_pos, & nonSchmid_pos, &
@ -87,6 +84,10 @@ module subroutine plastic_disloUCLA_init
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
N_sl N_sl
real, dimension(:), allocatable :: &
rho_mob_0, & !< initial dislocation density
rho_dip_0, & !< initial dipole density
a !< non-Schmid coefficients
character(len=pStringLen) :: & character(len=pStringLen) :: &
extmsg = '' extmsg = ''
@ -104,7 +105,6 @@ module subroutine plastic_disloUCLA_init
allocate(dotState(Ninstance)) allocate(dotState(Ninstance))
allocate(dependentState(Ninstance)) allocate(dependentState(Ninstance))
do p = 1, size(phase_plasticity) do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_DISLOUCLA_ID) cycle if (phase_plasticity(p) /= PLASTICITY_DISLOUCLA_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), & associate(prm => param(phase_plasticityInstance(p)), &
@ -127,25 +127,22 @@ module subroutine plastic_disloUCLA_init
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& a = config%getFloats('nonschmid_coefficients',defaultVal = emptyRealArray)
defaultVal = emptyRealArray) prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,+1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,-1)
else else
allocate(prm%nonSchmidCoeff(0)) prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_pos = prm%Schmid prm%nonSchmid_neg = prm%Schmid
prm%nonSchmid_neg = prm%Schmid
endif endif
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl, & prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,config%getFloats('interaction_slipslip'), &
config%getFloats('interaction_slipslip'), & config%getString('lattice_structure'))
config%getString('lattice_structure')) prm%forestProjection = lattice_forestProjection_edge(N_sl,config%getString('lattice_structure'),&
prm%forestProjectionEdge = lattice_forestProjection_edge(N_sl,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal))
config%getFloat('c/a',defaultVal=0.0_pReal)) prm%forestProjection = transpose(prm%forestProjection)
prm%forestProjectionEdge = transpose(prm%forestProjectionEdge)
prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(N_sl)) rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(N_sl))
prm%rho_dip_0 = config%getFloats('rhoedgedip0', requiredSize=size(N_sl)) rho_dip_0 = config%getFloats('rhoedgedip0', requiredSize=size(N_sl))
prm%v0 = config%getFloats('v0', requiredSize=size(N_sl)) prm%v0 = config%getFloats('v0', requiredSize=size(N_sl))
prm%b_sl = config%getFloats('slipburgers', requiredSize=size(N_sl)) prm%b_sl = config%getFloats('slipburgers', requiredSize=size(N_sl))
prm%delta_F = config%getFloats('qedge', requiredSize=size(N_sl)) prm%delta_F = config%getFloats('qedge', requiredSize=size(N_sl))
@ -161,16 +158,16 @@ module subroutine plastic_disloUCLA_init
prm%omega = config%getFloats('omega', requiredSize=size(N_sl)) prm%omega = config%getFloats('omega', requiredSize=size(N_sl))
prm%B = config%getFloats('friction_coeff', requiredSize=size(N_sl)) prm%B = config%getFloats('friction_coeff', requiredSize=size(N_sl))
prm%D = config%getFloat('grainsize') prm%D = config%getFloat('grainsize')
prm%D_0 = config%getFloat('d0') prm%D_0 = config%getFloat('d0')
prm%Q_cl = config%getFloat('qsd') prm%Q_cl = config%getFloat('qsd')
prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal
prm%D_a = config%getFloat('cedgedipmindistance') * prm%b_sl prm%D_a = config%getFloat('cedgedipmindistance') * prm%b_sl
prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key
! expand: family => system ! expand: family => system
prm%rho_mob_0 = math_expand(prm%rho_mob_0, N_sl) rho_mob_0 = math_expand(rho_mob_0, N_sl)
prm%rho_dip_0 = math_expand(prm%rho_dip_0, N_sl) rho_dip_0 = math_expand(rho_dip_0, N_sl)
prm%q = math_expand(prm%q, N_sl) prm%q = math_expand(prm%q, N_sl)
prm%p = math_expand(prm%p, N_sl) prm%p = math_expand(prm%p, N_sl)
prm%delta_F = math_expand(prm%delta_F, N_sl) prm%delta_F = math_expand(prm%delta_F, N_sl)
@ -185,22 +182,25 @@ module subroutine plastic_disloUCLA_init
prm%atomicVolume = math_expand(prm%atomicVolume, N_sl) prm%atomicVolume = math_expand(prm%atomicVolume, N_sl)
prm%D_a = math_expand(prm%D_a, N_sl) prm%D_a = math_expand(prm%D_a, N_sl)
! sanity checks ! sanity checks
if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' 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 ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl'
if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0'
if (any(prm%rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0'
if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0'
if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' slipb_sl' if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge' if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge'
if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_0' if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_0'
if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipb_sl' if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or b_sl'
if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipb_sl' if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or b_sl'
else slipActive else slipActive
allocate(prm%rho_mob_0(0)) rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%rho_dip_0(0)) allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%atomicVolume,prm%tau_0, &
prm%delta_F,prm%v0,prm%p,prm%q,prm%B,prm%kink_height,prm%w,prm%omega, &
source = emptyRealArray)
allocate(prm%forestProjection(0,0))
allocate(prm%h_sl_sl (0,0))
endif slipActive endif slipActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -215,23 +215,23 @@ module subroutine plastic_disloUCLA_init
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) stt%rho_mob => plasticState(p)%state(startIndex:endIndex,:)
stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) stt%rho_mob = spread(rho_mob_0,2,NipcMyPhase)
dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) dot%rho_mob => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_rho',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) stt%rho_dip => plasticState(p)%state(startIndex:endIndex,:)
stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase) stt%rho_dip = spread(rho_dip_0,2,NipcMyPhase)
dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) dot%rho_dip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_rho',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_rho',defaultVal=1.0_pReal)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) stt%gamma_sl => plasticState(p)%state(startIndex:endIndex,:)
dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) dot%gamma_sl => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = 1.0e6_pReal ! ARRG plasticState(p)%atol(startIndex:endIndex) = 1.0e6_pReal ! ARRG
! global alias ! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
@ -256,7 +256,7 @@ end subroutine plastic_disloUCLA_init
!> @brief Calculate plastic velocity gradient and its tangent. !> @brief Calculate plastic velocity gradient and its tangent.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, &
Mp,T,instance,of) Mp,T,instance,of)
real(pReal), dimension(3,3), intent(out) :: & real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: & real(pReal), dimension(3,3,3,3), intent(out) :: &
@ -369,8 +369,7 @@ module subroutine plastic_disloUCLA_dependentState(instance,of)
associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) associate(prm => param(instance), stt => state(instance),dst => dependentState(instance))
dislocationSpacing = sqrt(matmul(prm%forestProjectionEdge, & dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,of)+stt%rho_dip(:,of)))
stt%rho_mob(:,of)+stt%rho_dip(:,of)))
dst%threshold_stress(:,of) = prm%mu*prm%b_sl & dst%threshold_stress(:,of) = prm%mu*prm%b_sl &
* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of))) * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of)))

93
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -8,24 +8,21 @@ submodule(constitutive) plastic_phenopowerlaw
type :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: &
gdot0_slip, & !< reference shear strain rate for slip gdot0_slip = 1.0_pReal, & !< reference shear strain rate for slip
gdot0_twin, & !< reference shear strain rate for twin gdot0_twin = 1.0_pReal, & !< reference shear strain rate for twin
n_slip, & !< stress exponent for slip n_slip = 1.0_pReal, & !< stress exponent for slip
n_twin, & !< stress exponent for twin n_twin = 1.0_pReal, & !< stress exponent for twin
spr, & !< push-up factor for slip saturation due to twinning spr = 1.0_pReal, & !< push-up factor for slip saturation due to twinning
c_1, & c_1 = 1.0_pReal, &
c_2, & c_2 = 1.0_pReal, &
c_3, & c_3 = 1.0_pReal, &
c_4, & c_4 = 1.0_pReal, &
h0_SlipSlip, & !< reference hardening slip - slip h0_SlipSlip = 1.0_pReal, & !< reference hardening slip - slip
h0_TwinSlip, & !< reference hardening twin - slip h0_TwinSlip = 1.0_pReal, & !< reference hardening twin - slip
h0_TwinTwin, & !< reference hardening twin - twin h0_TwinTwin = 1.0_pReal, & !< reference hardening twin - twin
a_slip a_slip = 1.0_pReal
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
xi_slip_0, & !< initial critical shear stress for slip
xi_twin_0, & !< initial critical shear stress for twin
xi_slip_sat, & !< maximum critical shear stress for slip xi_slip_sat, & !< maximum critical shear stress for slip
nonSchmidCoeff, &
H_int, & !< per family hardening activity (optional) H_int, & !< per family hardening activity (optional)
gamma_twin_char !< characteristic shear for twins gamma_twin_char !< characteristic shear for twins
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
@ -41,6 +38,8 @@ submodule(constitutive) plastic_phenopowerlaw
integer :: & integer :: &
sum_N_sl, & !< total number of active slip system sum_N_sl, & !< total number of active slip system
sum_N_tw !< total number of active twin systems sum_N_tw !< total number of active twin systems
logical :: &
nonSchmidActive = .false.
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
end type tParameters end type tParameters
@ -77,6 +76,10 @@ module subroutine plastic_phenopowerlaw_init
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
N_sl, N_tw N_sl, N_tw
real, dimension(:), allocatable :: &
xi_slip_0, & !< initial critical shear stress for slip
xi_twin_0, & !< initial critical shear stress for twin
a !< non-Schmid coefficients
character(len=pStringLen) :: & character(len=pStringLen) :: &
extmsg = '' extmsg = ''
@ -97,13 +100,6 @@ module subroutine plastic_phenopowerlaw_init
stt => state(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), &
config => config_phase(p)) config => config_phase(p))
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal)
prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal)
prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal)
prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters
N_sl = config%getInts('nslip',defaultVal=emptyIntArray) N_sl = config%getInts('nslip',defaultVal=emptyIntArray)
@ -113,12 +109,11 @@ module subroutine plastic_phenopowerlaw_init
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& a = config%getFloats('nonschmid_coefficients',defaultVal = emptyRealArray)
defaultVal = emptyRealArray) if(size(a) > 0) prm%nonSchmidActive = .true.
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,+1) prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,-1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
else else
allocate(prm%nonSchmidCoeff(0))
prm%nonSchmid_pos = prm%P_sl prm%nonSchmid_pos = prm%P_sl
prm%nonSchmid_neg = prm%P_sl prm%nonSchmid_neg = prm%P_sl
endif endif
@ -126,7 +121,7 @@ module subroutine plastic_phenopowerlaw_init
config%getFloats('interaction_slipslip'), & config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure')) config%getString('lattice_structure'))
prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(N_sl)) xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(N_sl))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(N_sl)) prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(N_sl))
prm%H_int = config%getFloats('h_int', requiredSize=size(N_sl), & prm%H_int = config%getFloats('h_int', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal,i=1,size(N_sl))]) defaultVal=[(0.0_pReal,i=1,size(N_sl))])
@ -137,7 +132,7 @@ module subroutine plastic_phenopowerlaw_init
prm%h0_SlipSlip = config%getFloat('h0_slipslip') prm%h0_SlipSlip = config%getFloat('h0_slipslip')
! expand: family => system ! expand: family => system
prm%xi_slip_0 = math_expand(prm%xi_slip_0, N_sl) xi_slip_0 = math_expand(xi_slip_0, N_sl)
prm%xi_slip_sat = math_expand(prm%xi_slip_sat,N_sl) prm%xi_slip_sat = math_expand(prm%xi_slip_sat,N_sl)
prm%H_int = math_expand(prm%H_int, N_sl) prm%H_int = math_expand(prm%H_int, N_sl)
@ -145,11 +140,13 @@ module subroutine plastic_phenopowerlaw_init
if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip' if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip'
if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip' if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip'
if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip' if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip'
if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_0' if (any(xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_0'
if (any(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_sat' if (any(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_sat'
else slipActive else slipActive
xi_slip_0 = emptyRealArray
allocate(prm%xi_slip_sat,prm%H_int,source=emptyRealArray)
allocate(prm%interaction_SlipSlip(0,0)) allocate(prm%interaction_SlipSlip(0,0))
allocate(prm%xi_slip_0(0))
endif slipActive endif slipActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -165,23 +162,28 @@ module subroutine plastic_phenopowerlaw_init
prm%gamma_twin_char = lattice_characteristicShear_twin(N_tw,config%getString('lattice_structure'),& prm%gamma_twin_char = lattice_characteristicShear_twin(N_tw,config%getString('lattice_structure'),&
config%getFloat('c/a')) config%getFloat('c/a'))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(N_tw)) xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(N_tw))
prm%gdot0_twin = config%getFloat('gdot0_twin') prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal)
prm%n_twin = config%getFloat('n_twin') prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal)
prm%spr = config%getFloat('s_pr') prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal)
prm%h0_TwinTwin = config%getFloat('h0_twintwin') prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal)
prm%gdot0_twin = config%getFloat('gdot0_twin')
prm%n_twin = config%getFloat('n_twin')
prm%spr = config%getFloat('s_pr')
prm%h0_TwinTwin = config%getFloat('h0_twintwin')
! expand: family => system ! expand: family => system
prm%xi_twin_0 = math_expand(prm%xi_twin_0, N_tw) xi_twin_0 = math_expand(xi_twin_0,N_tw)
! sanity checks ! sanity checks
if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin' if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin'
if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin' if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin'
else twinActive else twinActive
xi_twin_0 = emptyRealArray
allocate(prm%gamma_twin_char,source=emptyRealArray)
allocate(prm%interaction_TwinTwin(0,0)) allocate(prm%interaction_TwinTwin(0,0))
allocate(prm%xi_twin_0(0))
allocate(prm%gamma_twin_char(0))
endif twinActive endif twinActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -218,7 +220,7 @@ module subroutine plastic_phenopowerlaw_init
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase) stt%xi_slip = spread(xi_slip_0, 2, NipcMyPhase)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_xi',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
@ -226,7 +228,7 @@ module subroutine plastic_phenopowerlaw_init
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase) stt%xi_twin = spread(xi_twin_0, 2, NipcMyPhase)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_xi',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
@ -431,20 +433,17 @@ pure subroutine kinetics_slip(Mp,instance,of, &
tau_slip_pos, & tau_slip_pos, &
tau_slip_neg tau_slip_neg
integer :: i integer :: i
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance)) associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1, prm%sum_N_sl do i = 1, prm%sum_N_sl
tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))
tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), &
0.0_pReal, nonSchmidActive) 0.0_pReal, prm%nonSchmidActive)
enddo enddo
where(dNeq0(tau_slip_pos)) where(dNeq0(tau_slip_pos))
gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos)
else where else where
gdot_slip_pos = 0.0_pReal gdot_slip_pos = 0.0_pReal