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

standard notation

This commit is contained in:
Martin Diehl 2021-07-21 16:34:33 +02:00
parent e3e9dfdc60
commit 76572e4ca9
4 changed files with 98 additions and 98 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -321,7 +321,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
real(pReal) :: & real(pReal) :: &
VacancyDiffusion VacancyDiffusion
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos, gdot_neg,& dot_gamma_pos, dot_gamma_neg,&
tau_pos,& tau_pos,&
tau_neg, & tau_neg, &
v_cl, & v_cl, &
@ -333,10 +333,10 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
dot => dotState(ph), dst => dependentState(ph)) dot => dotState(ph), dst => dependentState(ph))
call kinetics(Mp,T,ph,en,& call kinetics(Mp,T,ph,en,&
gdot_pos,gdot_neg, & dot_gamma_pos,dot_gamma_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg) tau_pos_out = tau_pos,tau_neg_out = tau_neg)
dot%gamma_sl(:,en) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs dot%gamma_sl(:,en) = (dot_gamma_pos+dot_gamma_neg) ! ToDo: needs to be abs
VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T))
where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg

60
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -240,22 +240,22 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos,gdot_neg, & dot_gamma_pos,dot_gamma_neg, &
dgdot_dtau_pos,dgdot_dtau_neg ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
associate(prm => param(ph)) associate(prm => param(ph))
call kinetics(Mp,ph,en,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
do i = 1, prm%sum_N_sl do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i) Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P(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) &
+ dgdot_dtau_pos(i) * prm%P(k,l,i) * prm%P_nS_pos(m,n,i) & + ddot_gamma_dtau_pos(i) * prm%P(k,l,i) * prm%P_nS_pos(m,n,i) &
+ dgdot_dtau_neg(i) * prm%P(k,l,i) * prm%P_nS_neg(m,n,i) + ddot_gamma_dtau_neg(i) * prm%P(k,l,i) * prm%P_nS_neg(m,n,i)
enddo enddo
end associate end associate
@ -277,14 +277,14 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en)
real(pReal) :: & real(pReal) :: &
sumGamma sumGamma
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos,gdot_neg dot_gamma_pos,dot_gamma_neg
associate(prm => param(ph), stt => state(ph),& associate(prm => param(ph), stt => state(ph),&
dot => dotState(ph)) dot => dotState(ph))
call kinetics(Mp,ph,en,gdot_pos,gdot_neg) call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg)
dot%accshear(:,en) = abs(gdot_pos+gdot_neg) dot%accshear(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
sumGamma = sum(stt%accshear(:,en)) sumGamma = sum(stt%accshear(:,en))
@ -318,15 +318,15 @@ module subroutine plastic_kinehardening_deltaState(Mp,ph,en)
en en
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos,gdot_neg, & dot_gamma_pos,dot_gamma_neg, &
sense sense
associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph)) associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph))
call kinetics(Mp,ph,en,gdot_pos,gdot_neg) call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg)
sense = merge(state(ph)%sense(:,en), & ! keep existing... sense = merge(state(ph)%sense(:,en), & ! keep existing...
sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined sign(1.0_pReal,dot_gamma_pos+dot_gamma_neg), & ! ...or have a defined
dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction dEq0(dot_gamma_pos+dot_gamma_neg,1e-10_pReal)) ! current sense of shear direction
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -392,7 +392,7 @@ end subroutine plastic_kinehardening_results
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics(Mp,ph,en, & pure subroutine kinetics(Mp,ph,en, &
gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -401,11 +401,11 @@ pure subroutine kinetics(Mp,ph,en, &
en en
real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: &
gdot_pos, & dot_gamma_pos, &
gdot_neg dot_gamma_neg
real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: & real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: &
dgdot_dtau_pos, & ddot_gamma_dtau_pos, &
dgdot_dtau_neg ddot_gamma_dtau_neg
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau_pos, & tau_pos, &
@ -421,31 +421,31 @@ pure subroutine kinetics(Mp,ph,en, &
enddo enddo
where(dNeq0(tau_pos)) where(dNeq0(tau_pos))
gdot_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active dot_gamma_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_pos/stt%crss(:,en))**prm%n, tau_pos) * sign(abs(tau_pos/stt%crss(:,en))**prm%n, tau_pos)
else where else where
gdot_pos = 0.0_pReal dot_gamma_pos = 0.0_pReal
end where end where
where(dNeq0(tau_neg)) where(dNeq0(tau_neg))
gdot_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 dot_gamma_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_neg/stt%crss(:,en))**prm%n, tau_neg) * sign(abs(tau_neg/stt%crss(:,en))**prm%n, tau_neg)
else where else where
gdot_neg = 0.0_pReal dot_gamma_neg = 0.0_pReal
end where end where
if (present(dgdot_dtau_pos)) then if (present(ddot_gamma_dtau_pos)) then
where(dNeq0(gdot_pos)) where(dNeq0(dot_gamma_pos))
dgdot_dtau_pos = gdot_pos*prm%n/tau_pos ddot_gamma_dtau_pos = dot_gamma_pos*prm%n/tau_pos
else where else where
dgdot_dtau_pos = 0.0_pReal ddot_gamma_dtau_pos = 0.0_pReal
end where end where
endif endif
if (present(dgdot_dtau_neg)) then if (present(ddot_gamma_dtau_neg)) then
where(dNeq0(gdot_neg)) where(dNeq0(dot_gamma_neg))
dgdot_dtau_neg = gdot_neg*prm%n/tau_neg ddot_gamma_dtau_neg = dot_gamma_neg*prm%n/tau_neg
else where else where
dgdot_dtau_neg = 0.0_pReal ddot_gamma_dtau_neg = 0.0_pReal
end where end where
endif endif
end associate end associate

52
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -773,7 +773,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
dv_dtauNS !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau, & !< resolved shear stress including backstress terms tau, & !< resolved shear stress including backstress terms
gdotTotal !< shear rate dot_gamma !< shear rate
associate(prm => param(ph),dst=>microstructure(ph),& associate(prm => param(ph),dst=>microstructure(ph),&
stt=>state(ph)) stt=>state(ph))
@ -823,12 +823,12 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
gdotTotal = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
do s = 1,ns do s = 1,ns
Lp = Lp + gdotTotal(s) * prm%Schmid(1:3,1:3,s) Lp = Lp + dot_gamma(s) * prm%Schmid(1:3,1:3,s)
forall (i=1:3,j=1:3,k=1:3,l=1:3) & forall (i=1:3,j=1:3,k=1:3,l=1:3) &
dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
+ prm%Schmid(i,j,s) * prm%Schmid(k,l,s) & + prm%Schmid(i,j,s) * prm%Schmid(k,l,s) &
@ -972,7 +972,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
real(pReal), dimension(param(ph)%sum_N_sl,4) :: & real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
v0, & v0, &
gdot !< shear rates dot_gamma !< shear rates
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau, & !< current resolved shear stress tau, & !< current resolved shear stress
v_climb !< climb velocity of edge dipoles v_climb !< climb velocity of edge dipoles
@ -995,7 +995,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
ns = prm%sum_N_sl ns = prm%sum_N_sl
tau = 0.0_pReal tau = 0.0_pReal
gdot = 0.0_pReal dot_gamma = 0.0_pReal
rho = getRho(ph,en) rho = getRho(ph,en)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
@ -1004,7 +1004,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
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) dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
@ -1030,17 +1030,17 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
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)
rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication rhoDotMultiplication(s,1:2) = sum(abs(dot_gamma(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation
rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication rhoDotMultiplication(s,3:4) = sum(abs(dot_gamma(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation
endforall endforall
else isBCC else isBCC
rhoDotMultiplication(:,1:4) = spread( & rhoDotMultiplication(:,1:4) = spread( &
(sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) & (sum(abs(dot_gamma(:,1:2)),2) * prm%f_ed_mult + sum(abs(dot_gamma(:,3:4)),2)) &
* sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26 * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26
endif isBCC endif isBCC
@ -1053,20 +1053,20 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
! 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(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c) = -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(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile + abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c))) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile * rhoSgl(:,2*c+3) * abs(dot_gamma(:,2*c)) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile * rhoSgl(:,2*c+4) * abs(dot_gamma(:,2*c-1)) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) &
+ abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) &
@ -1079,9 +1079,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
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 &
* ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single * ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1))) & ! was single hitting single
+ 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single + 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single
+ rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent + rhoDip(:,c) * (abs(dot_gamma(:,2*c-1)) + abs(dot_gamma(:,2*c)))) ! single knocks dipole constituent
! annihilated screw dipoles leave edge jogs behind on the colinear system ! annihilated screw dipoles leave edge jogs behind on the colinear system
if (phase_lattice(ph) == 'cF') & if (phase_lattice(ph) == 'cF') &
@ -1118,7 +1118,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
else else
dot%rho(:,en) = pack(rhoDot,.true.) dot%rho(:,en) = pack(rhoDot,.true.)
dot%gamma(:,en) = sum(gdot,2) dot%gamma(:,en) = sum(dot_gamma,2)
endif endif
end associate end associate
@ -1168,7 +1168,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
v0, & v0, &
neighbor_v0, & !< dislocation glide velocity of enighboring ip neighbor_v0, & !< dislocation glide velocity of enighboring ip
gdot !< shear rates dot_gamma !< shear rates
real(pReal), dimension(3,param(ph)%sum_N_sl,4) :: & real(pReal), dimension(3,param(ph)%sum_N_sl,4) :: &
m !< direction of dislocation motion m !< direction of dislocation motion
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
@ -1194,7 +1194,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
stt => state(ph)) stt => state(ph))
ns = prm%sum_N_sl ns = prm%sum_N_sl
gdot = 0.0_pReal dot_gamma = 0.0_pReal
rho = getRho(ph,en) rho = getRho(ph,en)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
@ -1202,7 +1202,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
@ -1212,14 +1212,14 @@ function rhoDotFlux(timestep,ph,en,ip,el)
if (plasticState(ph)%nonlocal) then if (plasticState(ph)%nonlocal) then
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ...
.and. prm%C_CFL * abs(v0) * timestep & .and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG #ifdef DEBUG
if (debugConstitutive%extensive) then if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', & print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', &
maxval(abs(v0), abs(gdot) > 0.0_pReal & maxval(abs(v0), abs(dot_gamma) > 0.0_pReal &
.and. prm%C_CFL * abs(v0) * timestep & .and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el))), & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
' at a timestep of ',timestep ' at a timestep of ',timestep

78
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -287,31 +287,31 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_sl_pos,gdot_sl_neg, & dot_gamma_sl_pos,dot_gamma_sl_neg, &
dgdot_dtauslip_pos,dgdot_dtauslip_neg ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
gdot_tw,dgdot_dtautwin dot_gamma_tw,ddot_gamma_dtautwin
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
associate(prm => param(ph)) associate(prm => param(ph))
call kinetics_sl(Mp,ph,en,gdot_sl_pos,gdot_sl_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg,ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg)
slipSystems: do i = 1, prm%sum_N_sl slipSystems: do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_sl_pos(i)+gdot_sl_neg(i))*prm%P_sl(1:3,1:3,i) Lp = Lp + (dot_gamma_sl_pos(i)+dot_gamma_sl_neg(i))*prm%P_sl(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) &
+ dgdot_dtauslip_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) & + ddot_gamma_dtau_sl_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) &
+ dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i) + ddot_gamma_dtau_sl_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i)
enddo slipSystems enddo slipSystems
call kinetics_tw(Mp,ph,en,gdot_tw,dgdot_dtautwin) call kinetics_tw(Mp,ph,en,dot_gamma_tw,ddot_gamma_dtautwin)
twinSystems: do i = 1, prm%sum_N_tw twinSystems: do i = 1, prm%sum_N_tw
Lp = Lp + gdot_tw(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) &
+ dgdot_dtautwin(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + ddot_gamma_dtautwin(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinSystems enddo twinSystems
end associate end associate
@ -334,14 +334,14 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en)
xi_sl_sat_offset,& xi_sl_sat_offset,&
sumF sumF
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_sl_pos,gdot_sl_neg, & dot_gamma_sl_pos,dot_gamma_sl_neg, &
right_SlipSlip right_SlipSlip
associate(prm => param(ph), stt => state(ph), dot => dotState(ph)) associate(prm => param(ph), stt => state(ph), dot => dotState(ph))
call kinetics_sl(Mp,ph,en,gdot_sl_pos,gdot_sl_neg) call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg)
dot%gamma_sl(:,en) = abs(gdot_sl_pos+gdot_sl_neg) dot%gamma_sl(:,en) = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en)) call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en))
@ -406,7 +406,7 @@ end subroutine plastic_phenopowerlaw_results
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_sl(Mp,ph,en, & pure subroutine kinetics_sl(Mp,ph,en, &
gdot_sl_pos,gdot_sl_neg,dgdot_dtau_sl_pos,dgdot_dtau_sl_neg) dot_gamma_sl_pos,dot_gamma_sl_neg,ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -415,11 +415,11 @@ pure subroutine kinetics_sl(Mp,ph,en, &
en en
real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: &
gdot_sl_pos, & dot_gamma_sl_pos, &
gdot_sl_neg dot_gamma_sl_neg
real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: & real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: &
dgdot_dtau_sl_pos, & ddot_gamma_dtau_sl_pos, &
dgdot_dtau_sl_neg ddot_gamma_dtau_sl_neg
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau_sl_pos, & tau_sl_pos, &
@ -435,31 +435,31 @@ pure subroutine kinetics_sl(Mp,ph,en, &
enddo enddo
where(dNeq0(tau_sl_pos)) where(dNeq0(tau_sl_pos))
gdot_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active dot_gamma_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_sl_pos/stt%xi_sl(:,en))**prm%n_sl, tau_sl_pos) * sign(abs(tau_sl_pos/stt%xi_sl(:,en))**prm%n_sl, tau_sl_pos)
else where else where
gdot_sl_pos = 0.0_pReal dot_gamma_sl_pos = 0.0_pReal
end where end where
where(dNeq0(tau_sl_neg)) where(dNeq0(tau_sl_neg))
gdot_sl_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2 dot_gamma_sl_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_sl_neg/stt%xi_sl(:,en))**prm%n_sl, tau_sl_neg) * sign(abs(tau_sl_neg/stt%xi_sl(:,en))**prm%n_sl, tau_sl_neg)
else where else where
gdot_sl_neg = 0.0_pReal dot_gamma_sl_neg = 0.0_pReal
end where end where
if (present(dgdot_dtau_sl_pos)) then if (present(ddot_gamma_dtau_sl_pos)) then
where(dNeq0(gdot_sl_pos)) where(dNeq0(dot_gamma_sl_pos))
dgdot_dtau_sl_pos = gdot_sl_pos*prm%n_sl/tau_sl_pos ddot_gamma_dtau_sl_pos = dot_gamma_sl_pos*prm%n_sl/tau_sl_pos
else where else where
dgdot_dtau_sl_pos = 0.0_pReal ddot_gamma_dtau_sl_pos = 0.0_pReal
end where end where
endif endif
if (present(dgdot_dtau_sl_neg)) then if (present(ddot_gamma_dtau_sl_neg)) then
where(dNeq0(gdot_sl_neg)) where(dNeq0(dot_gamma_sl_neg))
dgdot_dtau_sl_neg = gdot_sl_neg*prm%n_sl/tau_sl_neg ddot_gamma_dtau_sl_neg = dot_gamma_sl_neg*prm%n_sl/tau_sl_neg
else where else where
dgdot_dtau_sl_neg = 0.0_pReal ddot_gamma_dtau_sl_neg = 0.0_pReal
end where end where
endif endif
end associate end associate
@ -475,7 +475,7 @@ end subroutine kinetics_sl
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_tw(Mp,ph,en,& pure subroutine kinetics_tw(Mp,ph,en,&
gdot_tw,dgdot_dtau_tw) 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
@ -484,9 +484,9 @@ pure subroutine kinetics_tw(Mp,ph,en,&
en en
real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
gdot_tw dot_gamma_tw
real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: & real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: &
dgdot_dtau_tw ddot_gamma_dtau_tw
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
tau_tw tau_tw
@ -499,17 +499,17 @@ pure subroutine kinetics_tw(Mp,ph,en,&
enddo enddo
where(tau_tw > 0.0_pReal) where(tau_tw > 0.0_pReal)
gdot_tw = (1.0_pReal-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction dot_gamma_tw = (1.0_pReal-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction
* prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw * prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw
else where else where
gdot_tw = 0.0_pReal dot_gamma_tw = 0.0_pReal
end where end where
if (present(dgdot_dtau_tw)) then if (present(ddot_gamma_dtau_tw)) then
where(dNeq0(gdot_tw)) where(dNeq0(dot_gamma_tw))
dgdot_dtau_tw = gdot_tw*prm%n_tw/tau_tw ddot_gamma_dtau_tw = dot_gamma_tw*prm%n_tw/tau_tw
else where else where
dgdot_dtau_tw = 0.0_pReal ddot_gamma_dtau_tw = 0.0_pReal
end where end where
endif endif