This commit is contained in:
Martin Diehl 2021-06-24 15:25:43 +02:00
parent ae9a48c124
commit 12ca7428cf
1 changed files with 28 additions and 63 deletions

View File

@ -744,17 +744,6 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
enddo enddo
endif endif
#ifdef DEBUG
if (debugConstitutive%extensive &
.and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)&
.or. .not. debugConstitutive%selective)) then
print'(/,a,i8,1x,i2,1x,i1,/)', '<< CONST >> nonlocal_microstructure at el ip ',el,ip
print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,en)
print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,en)*1e-6
print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,en)*1e-6
endif
#endif
end associate end associate
end subroutine nonlocal_dependentState end subroutine nonlocal_dependentState
@ -1030,17 +1019,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
#ifdef DEBUG
if (debugConstitutive%basic &
.and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip) &
.or. .not. debugConstitutive%selective)) then
print'(a,/,10(12x,12(e12.5,1x),/))', '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip
print'(a,/,4(12x,12(e12.5,1x),/))', '<< CONST >> gdot / 1/s',gdot
endif
#endif
!****************************************************************************
!*** limits for stable dipole height ! limits for stable dipole height
do s = 1,ns do s = 1,ns
tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,en) tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,en)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
@ -1057,8 +1038,8 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
dUpper = max(dUpper,dLower) dUpper = max(dUpper,dLower)
!****************************************************************************
!*** dislocation multiplication ! dislocation multiplication
rhoDotMultiplication = 0.0_pReal rhoDotMultiplication = 0.0_pReal
isBCC: if (phase_lattice(ph) == 'cI') then isBCC: if (phase_lattice(ph) == 'cI') then
forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal)
@ -1080,9 +1061,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
!**************************************************************************** !****************************************************************************
!*** calculate dipole formation and annihilation ! dipole formation and annihilation
!*** formation by glide ! formation by glide
do c = 1,2 do c = 1,2
rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile
@ -1107,7 +1088,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
enddo enddo
!*** athermal annihilation ! athermal annihilation
rhoDotAthermalAnnihilation = 0.0_pReal rhoDotAthermalAnnihilation = 0.0_pReal
forall (c=1:2) & forall (c=1:2) &
rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl & rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl &
@ -1122,7 +1103,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
* 0.25_pReal * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed * 0.25_pReal * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed
!*** thermally activated annihilation of edge dipoles by climb ! thermally activated annihilation of edge dipoles by climb
rhoDotThermalAnnihilation = 0.0_pReal rhoDotThermalAnnihilation = 0.0_pReal
selfDiffusion = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) selfDiffusion = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature))
vClimb = prm%V_at * selfDiffusion * prm%mu & vClimb = prm%V_at * selfDiffusion * prm%mu &
@ -1650,13 +1631,13 @@ end subroutine stateInit
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates kinetics !> @brief calculates kinetics
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, ph) pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, ph)
integer, intent(in) :: & integer, intent(in) :: &
c, & !< dislocation character (1:edge, 2:screw) c, & !< dislocation character (1:edge, 2:screw)
ph ph
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
Temperature !< temperature T !< T
real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: &
tau, & !< resolved external shear stress (without non Schmid effects) tau, & !< resolved external shear stress (without non Schmid effects)
tauNS, & !< resolved external shear stress (including non Schmid effects) tauNS, & !< resolved external shear stress (including non Schmid effects)
@ -1677,16 +1658,11 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem
vViscous, & !< viscous glide velocity vViscous, & !< viscous glide velocity
dtPeierls_dtau, & !< derivative with respect to resolved shear stress dtPeierls_dtau, & !< derivative with respect to resolved shear stress
dtSolidSolution_dtau, & !< derivative with respect to resolved shear stress dtSolidSolution_dtau, & !< derivative with respect to resolved shear stress
meanfreepath_S, & !< mean free travel distance for dislocations between two solid solution obstacles lambda_S, & !< mean free distance between two solid solution obstacles
meanfreepath_P, & !< mean free travel distance for dislocations between two Peierls barriers lambda_P, & !< mean free distance between two Peierls barriers
jumpWidth_P, & !< depth of activated area
jumpWidth_S, & !< depth of activated area
activationLength_P, & !< length of activated dislocation line
activationLength_S, & !< length of activated dislocation line
activationVolume_P, & !< volume that needs to be activated to overcome barrier activationVolume_P, & !< volume that needs to be activated to overcome barrier
activationVolume_S, & !< volume that needs to be activated to overcome barrier activationVolume_S, & !< volume that needs to be activated to overcome barrier
activationEnergy_P, & !< energy that is needed to overcome barrier activationEnergy_P, & !< energy that is needed to overcome barrier
activationEnergy_S, & !< energy that is needed to overcome barrier
criticalStress_P, & !< maximum obstacle strength criticalStress_P, & !< maximum obstacle strength
criticalStress_S, & !< maximum obstacle strength criticalStress_S, & !< maximum obstacle strength
mobility !< dislocation mobility mobility !< dislocation mobility
@ -1700,40 +1676,32 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem
if (abs(tau(s)) > tauThreshold(s)) then if (abs(tau(s)) > tauThreshold(s)) then
!* Peierls contribution !* Peierls contribution
!* Effective stress includes non Schmid constributions tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s))
!* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity lambda_P = prm%b_sl(s)
tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive activationVolume_P = prm%w *prm%b_sl(s)**3.0_pReal
meanfreepath_P = prm%b_sl(s)
jumpWidth_P = prm%b_sl(s)
activationLength_P = prm%w *prm%b_sl(s)
activationVolume_P = activationLength_P * jumpWidth_P * prm%b_sl(s)
criticalStress_P = prm%peierlsStress(s,c) criticalStress_P = prm%peierlsStress(s,c)
activationEnergy_P = criticalStress_P * activationVolume_P activationEnergy_P = criticalStress_P * activationVolume_P
tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) ! ensure that the activation probability cannot become greater than one tauRel_P = min(1.0_pReal, tauEff / criticalStress_P)
tPeierls = 1.0_pReal / prm%nu_a & tPeierls = 1.0_pReal / prm%nu_a &
* exp(activationEnergy_P / (kB * Temperature) & * exp(activationEnergy_P / (kB * T) &
* (1.0_pReal - tauRel_P**prm%p)**prm%q) * (1.0_pReal - tauRel_P**prm%p)**prm%q)
if (tauEff < criticalStress_P) then if (tauEff < criticalStress_P) then
dtPeierls_dtau = tPeierls * prm%p * prm%q * activationVolume_P / (kB * Temperature) & dtPeierls_dtau = tPeierls * prm%p * prm%q * activationVolume_P / (kB * T) &
* (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal) * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal)
else else
dtPeierls_dtau = 0.0_pReal dtPeierls_dtau = 0.0_pReal
endif endif
!* Contribution from solid solution strengthening ! Contribution from solid solution strengthening
!* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity
tauEff = abs(tau(s)) - tauThreshold(s) tauEff = abs(tau(s)) - tauThreshold(s)
meanfreepath_S = prm%b_sl(s) / sqrt(prm%c_sol) lambda_S = prm%b_sl(s) / sqrt(prm%c_sol)
jumpWidth_S = prm%f_sol * prm%b_sl(s) activationVolume_S = prm%f_sol * prm%b_sl(s)**3.0_pReal / sqrt(prm%c_sol)
activationLength_S = prm%b_sl(s) / sqrt(prm%c_sol) criticalStress_S = prm%Q_sol / activationVolume_S
activationVolume_S = activationLength_S * jumpWidth_S * prm%b_sl(s) tauRel_S = min(1.0_pReal, tauEff / criticalStress_S)
activationEnergy_S = prm%Q_sol
criticalStress_S = activationEnergy_S / activationVolume_S
tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) ! ensure that the activation probability cannot become greater than one
tSolidSolution = 1.0_pReal / prm%nu_a & tSolidSolution = 1.0_pReal / prm%nu_a &
* exp(activationEnergy_S / (kB * Temperature)* (1.0_pReal - tauRel_S**prm%p)**prm%q) * exp(prm%Q_sol / (kB * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
if (tauEff < criticalStress_S) then if (tauEff < criticalStress_S) then
dtSolidSolution_dtau = tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * Temperature) & dtSolidSolution_dtau = tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * T) &
* (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal) * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal)
else else
dtSolidSolution_dtau = 0.0_pReal dtSolidSolution_dtau = 0.0_pReal
@ -1744,13 +1712,10 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem
mobility = prm%b_sl(s) / prm%eta mobility = prm%b_sl(s) / prm%eta
vViscous = mobility * tauEff vViscous = mobility * tauEff
!* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of
!* free flight at glide velocity in between.
!* adopt sign from resolved stress
v(s) = sign(1.0_pReal,tau(s)) & v(s) = sign(1.0_pReal,tau(s)) &
/ (tPeierls / meanfreepath_P + tSolidSolution / meanfreepath_S + 1.0_pReal / vViscous) / (tPeierls / lambda_P + tSolidSolution / lambda_S + 1.0_pReal / vViscous)
dv_dtau(s) = v(s)**2.0_pReal * (dtSolidSolution_dtau / meanfreepath_S + mobility /vViscous**2.0_pReal) dv_dtau(s) = v(s)**2.0_pReal * (dtSolidSolution_dtau / lambda_S + mobility /vViscous**2.0_pReal)
dv_dtauNS(s) = v(s)**2.0_pReal * dtPeierls_dtau / meanfreepath_P dv_dtauNS(s) = v(s)**2.0_pReal * dtPeierls_dtau / lambda_P
endif endif
enddo enddo