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

Merge branch 'dislotwin-correction' into 'development' 

dislotwin-correction

See merge request damask/DAMASK!347
This commit is contained in:
Franz Roters 2021-02-28 19:45:34 +00:00
parent 1be6950c0b d92393ab70
commit 7b89cb41f9
10 changed files with 121 additions and 179 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 0289c1bbfec1a1aef77a8cbaeed134035549e738 Subproject commit 48dd9972d9023caa8b04226112dcdd57fa0be6af

4
examples/ConfigFiles/Microstructure_DP_Steel.config
View File

@ -1,4 +0,0 @@
[DP_Steel]
crystallite 1
(constituent) phase 1 texture 1 fraction 0.82
(constituent) phase 2 texture 2 fraction 0.18

64
examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config
View File

@ -1,64 +0,0 @@
[TWIP_Steel_FeMnC]
elasticity hooke
plasticity dislotwin
(output) rho_mob
(output) rho_dip
(output) gamma_sl
(output) lambda_sl
(output) tau_pass
(output) f_tw
(output) lambda_tw
(output) tau_hat_tw
(output) f_tr
### Material parameters ###
lattice_structure fcc
C11 175.0e9 # From Music et al. Applied Physics Letters 91, 191904 (2007)
C12 115.0e9
C44 135.0e9
grainsize 2.0e-5 # Average grain size [m]
SolidSolutionStrength 1.5e8 # Strength due to elements in solid solution
### Dislocation glide parameters ###
Nslip 12
slipburgers 2.56e-10 # Burgers vector of slip system [m]
rhoedgedip0 1.0 # Initial dislocation density [m/m**3]
rhoedge0 1.0e12 # Initial dislocation density [m/m**3]
v0 1.0e-4 # Initial glide velocity [m/s]
Qedge 3.7e-19 # Activation energy for dislocation glide [J]
p_slip 1.0 # p-exponent in glide velocity
q_slip 1.0 # q-exponent in glide velocity
# hardening of glide
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^3]
Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b]
interactionSlipSlip 0.122 0.122 0.625 0.07 0.137 0.122 # Interaction coefficients (Kubin et al. 2008)
### Shearband parameters ###
shearbandresistance 180e6
shearbandvelocity 0e-4 # set to zero to turn shear banding of
QedgePerSbSystem 3.7e-19 # Activation energy for shear banding [J]
p_shearband 1.0 # p-exponent in glide velocity
q_shearband 1.0 # q-exponent in glide velocity
### Twinning parameters ###
Ntwin 12
twinburgers 1.47e-10 # Burgers vector of twin system [m]
twinsize 5.0e-8 # Twin stack mean thickness [m]
L0_twin 442.0 # Length of twin nuclei in Burgers vectors
maxtwinfraction 1.0 # Maximum admissible twin volume fraction
xc_twin 1.0e-9 # critical distance for formation of twin nucleus
VcrossSlip 1.67e-29 # cross slip volume
r_twin 10.0 # r-exponent in twin formation probability
Cmfptwin 1.0 # Adj. parameter controlling twin mean free path
Cthresholdtwin 1.0 # Adj. parameter controlling twin threshold stress
interactionSlipTwin 0.0 1.0 1.0 # Dislocation-Twin interaction coefficients
interactionTwinTwin 0.0 1.0 # Twin-Twin interaction coefficients
SFE_0K -0.0396 # stacking fault energy at zero K; TWIP steel: -0.0526; Cu: -0.0396
dSFE_dT 0.0002 # temperature dependance of stacking fault energy

41
examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml Normal file
View File

@ -0,0 +1,41 @@
TWIP_Steel_FeMnC:
lattice: cF
mechanics:
elasticity: {type: hooke, C_11: 175.0e9, C_12: 115.0e9, C_44: 135.0e9}
plasticity:
type: dislotwin
output: [rho_mob, rho_dip, gamma_sl, Lambda_sl, tau_pass, f_tw, Lambda_tw, tau_hat_tw, f_tr]
D: 2.0e-5
N_sl: [12]
b_sl: [2.56e-10]
rho_mob_0: [1.0e12]
rho_dip_0: [1.0]
v_0: [1.0e4]
Q_s: [3.7e-19]
p_sl: [1.0]
q_sl: [1.0]
tau_0: [1.5e8]
i_sl: [10.0] # Adj. parameter controlling dislocation mean free path
D_0: 4.0e-5 # Vacancy diffusion prefactor / m^2/s
D_a: 1.0 # minimum dipole distance / b
Q_cl: 4.5e-19 # Activation energy for climb / J
h_sl_sl: [0.122, 0.122, 0.625, 0.07, 0.137, 0.122] # Interaction coefficients (Kubin et al. 2008)
# shear band parameters
xi_sb: 180.0e6
Q_sb: 3.7e-19
p_sb: 1.0
q_sb: 1.0
v_sb: 0.0 # set to 0, to turn it off
# twinning parameters
N_tw: [12]
b_tw: [1.47e-10] # Burgers vector length of twin system / b
t_tw: [5.0e-8] # Twin stack mean thickness / m
L_tw: 442.0 # Length of twin nuclei / b
x_c_tw: 1.0e-9 # critical distance for formation of twin nucleus / m
V_cs: 1.67e-29 # cross slip volume / m^3
p_tw: [10.0] # r-exponent in twin formation probability
i_tw: 1.0 # Adj. parameter controlling twin mean free path
h_sl_tw: [0.0, 1.0, 1.0] # dislocation-twin interaction coefficients
h_tw_tw: [0.0, 1.0] # twin-twin interaction coefficients
Gamma_sf_0K: -0.0396 # stacking fault energy / J/m^2 at zero K; TWIP steel: -0.0526; Cu: -0.0396
dGamma_sf_dT: 0.0002 # temperature dependence / J/(m^2 K) of stacking fault energy

36
examples/ConfigFiles/Phase_Dislotwin_Tungsten.config
View File

@ -1,36 +0,0 @@
[Tungsten]
elasticity hooke
plasticity dislotwin
### Material parameters ###
lattice_structure bcc
C11 523.0e9 # From Marinica et al. Journal of Physics: Condensed Matter(2013)
C12 202.0e9
C44 161.0e9
grainsize 2.0e-5 # Average grain size [m]
SolidSolutionStrength 1.5e8 # 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]
v0 1.0e-4 # Initial glide velocity [m/s]
Qedge 2.725e-19 # Activation energy for dislocation glide [J]
p_slip 0.78 # p-exponent in glide velocity
q_slip 1.58 # q-exponent in glide velocity
tau_peierls 2.03e9 # peierls stress (for bcc)
dipoleformationfactor 0 # to have hardening due to dipole formation off
#hardening
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 1 1 1.4 1.4 1.4 1.4

21
examples/ConfigFiles/Phase_Dislotwin_Tungsten.yaml Normal file
View File

@ -0,0 +1,21 @@
Tungsten:
lattice: cI
mechanics:
elasticity: {type: hooke, C_11: 523.0e9, C_12: 202.0e9, C_44: 161.0e9} # Marinica et al. Journal of Physics: Condensed Matter(2013)
plasticity:
type: dislotwin
D: 2.0e-5 # Average grain size / m
N_sl: [12]
b_sl: [2.72e-10] # Burgers vector length of slip families / m
rho_mob_0: [1.0e12]
rho_dip_0: [1.0]
v_0: [1.0e4] # Initial glide velocity / m/s
Q_s: [2.725e-19] # Activation energy for dislocation glide / J
p_sl: [0.78] # p-exponent in glide velocity
q_sl: [1.58] # q-exponent in glide velocity
tau_0: [1.5e8] # solid solution strength / Pa
i_sl: [10.0] # Adj. parameter controlling dislocation mean free path
D_0: 4.0e-5 # Vacancy diffusion prefactor / m^2/s
D_a: 1.0 # minimum dipole distance / b
Q_cl: 4.5e-19 # Activation energy for climb / J
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]

3
examples/ConfigFiles/Phase_Hydrogen.config
View File

@ -1,3 +0,0 @@
hydrogenflux_diffusion11 1.0
hydrogenflux_mobility11 1.0
hydrogenVolume 1e-28

4
python/damask/util.py
View File

@ -133,6 +133,8 @@ def execute(cmd,
stdout = stdout.decode('utf-8').replace('\x08','') stdout = stdout.decode('utf-8').replace('\x08','')
stderr = stderr.decode('utf-8').replace('\x08','') stderr = stderr.decode('utf-8').replace('\x08','')
if process.returncode != 0: if process.returncode != 0:
print(stdout)
print(stderr)
raise RuntimeError(f"'{cmd}' failed with returncode {process.returncode}") raise RuntimeError(f"'{cmd}' failed with returncode {process.returncode}")
return stdout, stderr return stdout, stderr
@ -437,7 +439,7 @@ class _ProgressBar:
bar = '█' * filled_length + '░' * (self.bar_length - filled_length) bar = '█' * filled_length + '░' * (self.bar_length - filled_length)
delta_time = datetime.datetime.now() - self.start_time delta_time = datetime.datetime.now() - self.start_time
remaining_time = (self.total - (iteration+1)) * delta_time / (iteration+1) remaining_time = (self.total - (iteration+1)) * delta_time / (iteration+1)
remaining_time -= datetime.timedelta(microseconds=remaining_time.microseconds) # remove μs remaining_time -= datetime.timedelta(microseconds=remaining_time.microseconds) # remove μs
sys.stderr.write(f'\r{self.prefix} {bar} {fraction:>4.0%} ETA {remaining_time}') sys.stderr.write(f'\r{self.prefix} {bar} {fraction:>4.0%} ETA {remaining_time}')
sys.stderr.flush() sys.stderr.flush()

8
src/math.f90
View File

@ -1133,6 +1133,7 @@ real(pReal) pure function math_areaTriangle(v1,v2,v3)
real(pReal), dimension (3), intent(in) :: v1,v2,v3 real(pReal), dimension (3), intent(in) :: v1,v2,v3
math_areaTriangle = 0.5_pReal * norm2(math_cross(v1-v2,v1-v3)) math_areaTriangle = 0.5_pReal * norm2(math_cross(v1-v2,v1-v3))
end function math_areaTriangle end function math_areaTriangle
@ -1147,11 +1148,13 @@ real(pReal) pure elemental function math_clip(a, left, right)
real(pReal), intent(in) :: a real(pReal), intent(in) :: a
real(pReal), intent(in), optional :: left, right real(pReal), intent(in), optional :: left, right
math_clip = a math_clip = a
if (present(left)) math_clip = max(left,math_clip) if (present(left)) math_clip = max(left,math_clip)
if (present(right)) math_clip = min(right,math_clip) if (present(right)) math_clip = min(right,math_clip)
if (present(left) .and. present(right)) & if (present(left) .and. present(right)) then
math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right) if(left>right) error stop 'left > right'
endif
end function math_clip end function math_clip
@ -1182,6 +1185,7 @@ subroutine selfTest
integer :: d integer :: d
logical :: e logical :: e
if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - & if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - &
math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) & math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) &
error stop 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]' error stop 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]'

117
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -24,7 +24,6 @@ submodule(phase:plastic) dislotwin
q_sb = 1.0_pReal, & !< q-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 D_a = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance
i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning
tau_0 = 1.0_pReal, & !< strength due to elements in solid solution
L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors
L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors
x_c_tw = 1.0_pReal, & !< critical distance for formation of twin nucleus x_c_tw = 1.0_pReal, & !< critical distance for formation of twin nucleus
@ -53,6 +52,7 @@ submodule(phase:plastic) dislotwin
q, & !< q-exponent in glide velocity q, & !< q-exponent in glide velocity
r, & !< r-exponent in twin nucleation rate r, & !< r-exponent in twin nucleation rate
s, & !< s-exponent in trans nucleation rate s, & !< s-exponent in trans nucleation rate
tau_0, & !< strength due to elements in solid solution
gamma_char, & !< characteristic shear for twins gamma_char, & !< characteristic shear for twins
B !< drag coefficient B !< drag coefficient
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
@ -81,7 +81,7 @@ submodule(phase:plastic) dislotwin
logical :: & logical :: &
ExtendedDislocations, & !< consider split into partials for climb calculation ExtendedDislocations, & !< consider split into partials for climb calculation
fccTwinTransNucleation, & !< twinning and transformation models are for fcc fccTwinTransNucleation, & !< twinning and transformation models are for fcc
dipoleFormation !< flag indicating consideration of dipole formation omitDipoles !< flag controlling consideration of dipole formation
end type !< container type for internal constitutive parameters end type !< container type for internal constitutive parameters
type :: tDislotwinState type :: tDislotwinState
@ -213,10 +213,10 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl)) prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl))
prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl)) prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl))
prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl)) prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl))
prm%tau_0 = pl%get_asFloats('tau_0', requiredSize=size(N_sl))
prm%B = pl%get_asFloats('B', requiredSize=size(N_sl), & prm%B = pl%get_asFloats('B', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal, i=1,size(N_sl))]) defaultVal=[(0.0_pReal, i=1,size(N_sl))])
prm%tau_0 = pl%get_asFloat('tau_0')
prm%D_a = pl%get_asFloat('D_a') prm%D_a = pl%get_asFloat('D_a')
prm%D_0 = pl%get_asFloat('D_0') prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl') prm%Q_cl = pl%get_asFloat('Q_cl')
@ -226,7 +226,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%dGamma_sf_dT = pl%get_asFloat('dGamma_sf_dT') prm%dGamma_sf_dT = pl%get_asFloat('dGamma_sf_dT')
endif endif
prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation',defaultVal = .false.) prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.)
! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex)
! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 ! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
@ -242,6 +242,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%i_sl = math_expand(prm%i_sl, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%p = math_expand(prm%p, N_sl) prm%p = math_expand(prm%p, N_sl)
prm%q = math_expand(prm%q, 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%B = math_expand(prm%B, N_sl)
! sanity checks ! sanity checks
@ -443,15 +444,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:) stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:)
dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:) dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_twin',defaultVal=1.0e-7_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tw',defaultVal=1.0e-7_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_twin' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tw'
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tr endIndex = endIndex + prm%sum_N_tr
stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:) stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:)
dot%f_tr=>plasticState(ph)%dotState(startIndex:endIndex,:) dot%f_tr=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tr',defaultVal=1.0e-6_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tr'
allocate(dst%Lambda_sl (prm%sum_N_sl,Nconstituents),source=0.0_pReal) allocate(dst%Lambda_sl (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
allocate(dst%tau_pass (prm%sum_N_sl,Nconstituents),source=0.0_pReal) allocate(dst%tau_pass (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
@ -535,9 +536,9 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_sl,ddot_gamma_dtau_slip dot_gamma_sl,ddot_gamma_dtau_slip
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
dot_gamma_twin,ddot_gamma_dtau_twin dot_gamma_tw,ddot_gamma_dtau_tw
real(pReal), dimension(param(ph)%sum_N_tr) :: & real(pReal), dimension(param(ph)%sum_N_tr) :: &
dot_gamma_tr,ddot_gamma_dtau_trans dot_gamma_tr,ddot_gamma_dtau_tr
real(pReal):: dot_gamma_sb real(pReal):: dot_gamma_sb
real(pReal), dimension(3,3) :: eigVectors, P_sb real(pReal), dimension(3,3) :: eigVectors, P_sb
real(pReal), dimension(3) :: eigValues real(pReal), dimension(3) :: eigValues
@ -578,20 +579,20 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
+ ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
enddo slipContribution enddo slipContribution
call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin,ddot_gamma_dtau_twin) call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tw,ddot_gamma_dtau_tw)
twinContibution: do i = 1, prm%sum_N_tw twinContibution: do i = 1, prm%sum_N_tw
Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinContibution enddo twinContibution
call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_trans) call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_tr)
transContibution: do i = 1, prm%sum_N_tr transContibution: do i = 1, prm%sum_N_tr
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
enddo transContibution enddo transContibution
Lp = Lp * f_unrotated Lp = Lp * f_unrotated
@ -646,7 +647,6 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
f_unrotated, & f_unrotated, &
rho_dip_distance, & rho_dip_distance, &
v_cl, & !< climb velocity v_cl, & !< climb velocity
Gamma, & !< stacking fault energy
tau, & tau, &
sigma_cl, & !< climb stress sigma_cl, & !< climb stress
b_d !< ratio of Burgers vector to stacking fault width b_d !< ratio of Burgers vector to stacking fault width
@ -656,7 +656,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
rho_dip_distance_min, & rho_dip_distance_min, &
dot_gamma_sl dot_gamma_sl
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
dot_gamma_twin dot_gamma_tw
real(pReal), dimension(param(ph)%sum_N_tr) :: & real(pReal), dimension(param(ph)%sum_N_tr) :: &
dot_gamma_tr dot_gamma_tr
@ -675,7 +675,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
slipState: do i = 1, prm%sum_N_sl slipState: do i = 1, prm%sum_N_sl
tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
significantSlipStress: if (dEq0(tau)) then significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then
dot_rho_dip_formation(i) = 0.0_pReal dot_rho_dip_formation(i) = 0.0_pReal
dot_rho_dip_climb(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal
else significantSlipStress else significantSlipStress
@ -683,24 +683,18 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,me)) rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,me))
rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i)) rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i))
if (prm%dipoleFormation) then 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*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) & * stt%rho_mob(i,me)*abs(dot_gamma_sl(i))
* stt%rho_mob(i,me)*abs(dot_gamma_sl(i))
else
dot_rho_dip_formation(i) = 0.0_pReal
endif
if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then
dot_rho_dip_climb(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal
else else
!@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 ! 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))) sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
if (prm%ExtendedDislocations) then b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) &
Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T * (prm%Gamma_sf_0K + prm%dGamma_sf_dT * T) / (prm%mu*prm%b_sl(i)), &
b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i)) 1.0_pReal, &
else prm%ExtendedDislocations)
b_d = 1.0_pReal
endif
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & 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) * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
@ -718,8 +712,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) & - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) &
- dot_rho_dip_climb - dot_rho_dip_climb
call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin) call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tw)
dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char dot%f_tw(:,me) = f_unrotated*dot_gamma_tw/prm%gamma_char
call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr) call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr)
dot%f_tr(:,me) = f_unrotated*dot_gamma_tr dot%f_tr(:,me) = f_unrotated*dot_gamma_tr
@ -741,11 +735,9 @@ module subroutine dislotwin_dependentState(T,ph,me)
T T
real(pReal) :: & real(pReal) :: &
sumf_twin,Gamma,sumf_trans sumf_tw,Gamma,sumf_tr
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation inv_lambda_sl
inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
f_over_t_tw f_over_t_tw
@ -760,38 +752,27 @@ module subroutine dislotwin_dependentState(T,ph,me)
stt => state(ph),& stt => state(ph),&
dst => dependentState(ph)) dst => dependentState(ph))
sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,me)) sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,me))
sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,me)) sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,me))
Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T
!* rescaled volume fraction for topology !* rescaled volume fraction for topology
f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,me)/prm%t_tw ! this is per system ... f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,me)/prm%t_tw ! this is per system ...
f_over_t_tr = sumf_trans/prm%t_tr ! but this not f_over_t_tr = sumf_tr/prm%t_tr ! but this not
! ToDo ...Physically correct, but naming could be adjusted ! ToDo ...Physically correct, but naming could be adjusted
inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, & inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl
stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl
if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_twin)
if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl_tr = matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
dst%Lambda_sl(:,me) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl)
inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_trans)
if ((prm%sum_N_tw > 0) .or. (prm%sum_N_tr > 0)) then ! ToDo: better logic needed here
dst%Lambda_sl(:,me) = prm%D &
/ (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr))
else
dst%Lambda_sl(:,me) = prm%D &
/ (1.0_pReal+prm%D*inv_lambda_sl_sl) !!!!!! correct?
endif
inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
dst%Lambda_tw(:,me) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) dst%Lambda_tw(:,me) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
dst%Lambda_tr(:,me) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) dst%Lambda_tr(:,me) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
!* threshold stress for dislocation motion !* threshold stress for dislocation motion
@ -957,7 +938,7 @@ end subroutine kinetics_slip
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,& pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
dot_gamma_twin,ddot_gamma_dtau_twin) dot_gamma_tw,ddot_gamma_dtau_tw)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -970,9 +951,9 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
dot_gamma_sl dot_gamma_sl
real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
dot_gamma_twin dot_gamma_tw
real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: &
ddot_gamma_dtau_twin ddot_gamma_dtau_tw
real, dimension(param(ph)%sum_N_tw) :: & real, dimension(param(ph)%sum_N_tw) :: &
tau, & tau, &
@ -1004,16 +985,16 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
significantStress: where(tau > tol_math_check) significantStress: where(tau > tol_math_check)
StressRatio_r = (dst%tau_hat_tw(:,me)/tau)**prm%r StressRatio_r = (dst%tau_hat_tw(:,me)/tau)**prm%r
dot_gamma_twin = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r) dot_gamma_tw = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r)
ddot_gamma_dtau = (dot_gamma_twin*prm%r/tau)*StressRatio_r ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r
else where significantStress else where significantStress
dot_gamma_twin = 0.0_pReal dot_gamma_tw = 0.0_pReal
ddot_gamma_dtau = 0.0_pReal ddot_gamma_dtau = 0.0_pReal
end where significantStress end where significantStress
end associate end associate
if(present(ddot_gamma_dtau_twin)) ddot_gamma_dtau_twin = ddot_gamma_dtau if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau
end subroutine kinetics_twin end subroutine kinetics_twin
@ -1026,7 +1007,7 @@ end subroutine kinetics_twin
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,& pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
dot_gamma_tr,ddot_gamma_dtau_trans) dot_gamma_tr,ddot_gamma_dtau_tr)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -1041,7 +1022,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
real(pReal), dimension(param(ph)%sum_N_tr), intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tr), intent(out) :: &
dot_gamma_tr dot_gamma_tr
real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: &
ddot_gamma_dtau_trans ddot_gamma_dtau_tr
real, dimension(param(ph)%sum_N_tr) :: & real, dimension(param(ph)%sum_N_tr) :: &
tau, & tau, &
@ -1081,7 +1062,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
end associate end associate
if(present(ddot_gamma_dtau_trans)) ddot_gamma_dtau_trans = ddot_gamma_dtau if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau
end subroutine kinetics_trans end subroutine kinetics_trans