bug fixes and polishing

bugs: wrong bracket in calculation of threshold_stress_twin
      3**b_twin instead of 3*b_twin in calculation of threshold_stress_twin
      ph instead of instance in calculation of mfp_twin
This commit is contained in:
Martin Diehl 2018-07-16 19:53:26 +02:00
parent 84c29bb255
commit fc0499a9bb
1 changed files with 38 additions and 54 deletions

View File

@ -347,6 +347,7 @@ subroutine plastic_dislotwin_init(fileUnit)
prm%twinsize= math_expand(prm%twinsize,prm%Ntwin) prm%twinsize= math_expand(prm%twinsize,prm%Ntwin)
prm%r = config_phase(p)%getFloats('r_twin') prm%r = config_phase(p)%getFloats('r_twin')
prm%r = math_expand(prm%r,prm%Ntwin)
prm%L0_twin = config_phase(p)%getFloat('l0_twin') prm%L0_twin = config_phase(p)%getFloat('l0_twin')
@ -1042,8 +1043,7 @@ subroutine plastic_dislotwin_microstructure(temperature,ipc,ip,el)
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
x0 x0
real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%Ntwin) :: fOverStacksize real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%Ntwin) :: fOverStacksize
real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%Ntrans) :: & real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%Ntrans) :: ftransOverLamellarSize
ftransOverLamellarSize
type(tParameters), pointer :: prm type(tParameters), pointer :: prm
@ -1061,13 +1061,11 @@ subroutine plastic_dislotwin_microstructure(temperature,ipc,ip,el)
sum(state(instance)%strainTransFraction(1_pInt:prm%totalNtrans,of)) sum(state(instance)%strainTransFraction(1_pInt:prm%totalNtrans,of))
!* Stacking fault energy !* Stacking fault energy
sfe = param(instance)%SFE_0K + & sfe = param(instance)%SFE_0K + param(instance)%dSFE_dT * Temperature
param(instance)%dSFE_dT * Temperature
!* rescaled twin volume fraction for topology !* rescaled twin volume fraction for topology
fOverStacksize = & fOverStacksize = state(instance)%twinFraction(1_pInt:prm%totalNtwin,of)/prm%twinsize
state(instance)%twinFraction(1_pInt:prm%totalNtwin,of)/prm%twinsize
!* rescaled trans volume fraction for topology !* rescaled trans volume fraction for topology
ftransOverLamellarSize = & ftransOverLamellarSize = &
@ -1111,77 +1109,62 @@ subroutine plastic_dislotwin_microstructure(temperature,ipc,ip,el)
do s = 1_pInt,prm%totalNslip do s = 1_pInt,prm%totalNslip
if ((prm%totalNtwin > 0_pInt) .or. (prm%totalNtrans > 0_pInt)) then if ((prm%totalNtwin > 0_pInt) .or. (prm%totalNtrans > 0_pInt)) then
state(instance)%mfp_slip(s,of) = & state(instance)%mfp_slip(s,of) = &
param(instance)%GrainSize/(1.0_pReal+param(instance)%GrainSize*& prm%GrainSize/(1.0_pReal+prm%GrainSize*&
(state(instance)%invLambdaSlip(s,of) + & (state(instance)%invLambdaSlip(s,of) + state(instance)%invLambdaSlipTwin(s,of) + state(instance)%invLambdaSlipTrans(s,of)))
state(instance)%invLambdaSlipTwin(s,of) + &
state(instance)%invLambdaSlipTrans(s,of)))
else else
state(instance)%mfp_slip(s,of) = & state(instance)%mfp_slip(s,of) = &
param(instance)%GrainSize/& prm%GrainSize/&
(1.0_pReal+param(instance)%GrainSize*(state(instance)%invLambdaSlip(s,of))) !!!!!! correct? (1.0_pReal+prm%GrainSize*(state(instance)%invLambdaSlip(s,of))) !!!!!! correct?
endif endif
enddo enddo
!* mean free path between 2 obstacles seen by a growing twin !* mean free path between 2 obstacles seen by a growing twin
state(instance)%mfp_twin(:,of) = prm%Cmfptwin*prm%GrainSize/&
state(instance)%mfp_twin(:,of) = & (1.0_pReal+prm%GrainSize*state(instance)%invLambdaTwin(:,of))
param(instance)%Cmfptwin*param(instance)%GrainSize/&
(1.0_pReal+param(instance)%GrainSize*state(ph)%invLambdaTwin(:,of))
!* mean free path between 2 obstacles seen by a growing martensite !* mean free path between 2 obstacles seen by a growing martensite
state(instance)%mfp_trans(:,of) = prm%Cmfptrans*prm%GrainSize/&
state(instance)%mfp_trans(:,of) = & (1.0_pReal+prm%GrainSize*state(instance)%invLambdaTrans(:,of))
param(instance)%Cmfptrans*param(instance)%GrainSize/&
(1.0_pReal+param(instance)%GrainSize*state(instance)%invLambdaTrans(:,of))
!* threshold stress for dislocation motion !* threshold stress for dislocation motion
forall (s = 1_pInt:prm%totalNslip) & forall (s = 1_pInt:prm%totalNslip) &
state(instance)%threshold_stress_slip(s,of) = & state(instance)%threshold_stress_slip(s,of) = &
lattice_mu(ph)*prm%burgers_slip(s)*& lattice_mu(ph)*prm%burgers_slip(s)*&
sqrt(dot_product((state(instance)%rhoEdge(1_pInt:prm%totalNslip,of)+state(instance)%rhoEdgeDip(1_pInt:prm%totalNslip,of)),& sqrt(dot_product(state(instance)%rhoEdge(1_pInt:prm%totalNslip,of)+state(instance)%rhoEdgeDip(1_pInt:prm%totalNslip,of),&
prm%interaction_SlipSlip(s,1:prm%totalNslip))) prm%interaction_SlipSlip(s,1:prm%totalNslip)))
!* threshold stress for growing twin !* threshold stress for growing twin
state(instance)%threshold_stress_twin(:,of) = & state(instance)%threshold_stress_twin(:,of) = &
param(instance)%Cthresholdtwin* & prm%Cthresholdtwin* &
(sfe/(3.0_pReal**prm%burgers_twin & (sfe/(3.0_pReal*prm%burgers_twin) &
+ 3.0_pReal*prm%burgers_twin*lattice_mu(ph)/& + 3.0_pReal*prm%burgers_twin*lattice_mu(ph)/&
(param(instance)%L0_twin*prm%burgers_slip)) & (param(instance)%L0_twin*prm%burgers_slip))
)
!* threshold stress for growing martensite !* threshold stress for growing martensite
state(instance)%threshold_stress_trans(:,of) = & state(instance)%threshold_stress_trans(:,of) = &
param(instance)%Cthresholdtrans* & prm%Cthresholdtrans* &
(sfe/(3.0_pReal*prm%burgers_trans) & (sfe/(3.0_pReal*prm%burgers_trans) &
+ 3.0_pReal*prm%burgers_trans*lattice_mu(ph)/& + 3.0_pReal*prm%burgers_trans*lattice_mu(ph)/&
(param(instance)%L0_trans*prm%burgers_slip)& (prm%L0_trans*prm%burgers_slip)&
+ param(instance)%transStackHeight*param(instance)%deltaG/ & + prm%transStackHeight*prm%deltaG/ &
(3.0_pReal*prm%burgers_trans) & (3.0_pReal*prm%burgers_trans) &
) )
!* final twin volume after growth !* final twin volume after growth
state(instance)%twinVolume(:,of) = & state(instance)%twinVolume(:,of) = (PI/4.0_pReal)*prm%twinsize*state(instance)%mfp_twin(:,of)**2.0_pReal
(pi/4.0_pReal)*prm%twinsize*&
state(instance)%mfp_twin(:,of)**2.0_pReal
!* final martensite volume after growth !* final martensite volume after growth
state(instance)%martensiteVolume(:,of) = & state(instance)%martensiteVolume(:,of) = (PI/4.0_pReal)*prm%lamellarsizePerTransSystem*state(instance)%mfp_trans(:,of)**2.0_pReal
(pi/4.0_pReal)*prm%lamellarsizePerTransSystem*&
state(instance)%mfp_trans(:,of)**(2.0_pReal)
!* equilibrium separation of partial dislocations (twin) !* equilibrium separation of partial dislocations (twin)
x0 = lattice_mu(ph)*prm%burgers_twin**(2.0_pReal)/& x0 = lattice_mu(ph)*prm%burgers_twin**2.0_pReal/(sfe*8.0_pReal*PI)*(2.0_pReal+lattice_nu(ph))/(1.0_pReal-lattice_nu(ph))
(sfe*8.0_pReal*pi)*(2.0_pReal+lattice_nu(ph))/(1.0_pReal-lattice_nu(ph)) tau_r_twin(:,instance)= lattice_mu(ph)*prm%burgers_twin/(2.0_pReal*PI)*&
tau_r_twin(:,instance)= & (1/(x0+prm%xc_twin)+cos(pi/3.0_pReal)/x0)
lattice_mu(ph)*prm%burgers_twin/(2.0_pReal*pi)*&
(1/(x0+param(instance)%xc_twin)+cos(pi/3.0_pReal)/x0)
!* equilibrium separation of partial dislocations (trans)
x0 = lattice_mu(ph)*prm%burgers_trans**(2.0_pReal)/& !* equilibrium separation of partial dislocations (trans)
(sfe*8.0_pReal*pi)*(2.0_pReal+lattice_nu(ph))/(1.0_pReal-lattice_nu(ph)) x0 = lattice_mu(ph)*prm%burgers_trans**2.0_pReal/(sfe*8.0_pReal*PI)*(2.0_pReal+lattice_nu(ph))/(1.0_pReal-lattice_nu(ph))
tau_r_trans(:,instance)= & tau_r_trans(:,instance)= lattice_mu(ph)*prm%burgers_trans/(2.0_pReal*PI)*&
lattice_mu(ph)*prm%burgers_trans/(2.0_pReal*pi)*& (1/(x0+prm%xc_trans)+cos(pi/3.0_pReal)/x0)
(1/(x0+param(instance)%xc_trans)+cos(pi/3.0_pReal)/x0)
end subroutine plastic_dislotwin_microstructure end subroutine plastic_dislotwin_microstructure
@ -1296,7 +1279,7 @@ prm => param(instance)
if((abs(tau_slip(j))-state(instance)%threshold_stress_slip(j,of)) > tol_math_check) then if((abs(tau_slip(j))-state(instance)%threshold_stress_slip(j,of)) > tol_math_check) then
!* Stress ratios !* Stress ratios
stressRatio =((abs(tau_slip(j))- state(instance)%threshold_stress_slip(j,of))/& stressRatio =((abs(tau_slip(j))- state(instance)%threshold_stress_slip(j,of))/&
(param(instance)%SolidSolutionStrength+prm%tau_peierlsPerSlipFamily(f))) (prm%SolidSolutionStrength+prm%tau_peierlsPerSlipFamily(f)))
StressRatio_p = stressRatio** prm%p(f) StressRatio_p = stressRatio** prm%p(f)
StressRatio_pminus1 = stressRatio**(prm%p(f)-1.0_pReal) StressRatio_pminus1 = stressRatio**(prm%p(f)-1.0_pReal)
!* Boltzmann ratio !* Boltzmann ratio
@ -1315,7 +1298,7 @@ prm => param(instance)
dgdot_dtauslip(j) = & dgdot_dtauslip(j) = &
abs(gdot_slip(j))*BoltzmannRatio*prm%p(f)& abs(gdot_slip(j))*BoltzmannRatio*prm%p(f)&
*prm%q(f)/& *prm%q(f)/&
(param(instance)%SolidSolutionStrength+prm%tau_peierlsPerSlipFamily(f))*& (prm%SolidSolutionStrength+prm%tau_peierlsPerSlipFamily(f))*&
StressRatio_pminus1*(1-StressRatio_p)**(prm%q(f)-1.0_pReal) StressRatio_pminus1*(1-StressRatio_p)**(prm%q(f)-1.0_pReal)
endif endif
@ -1434,7 +1417,7 @@ prm => param(instance)
gdot_twin(j) = & gdot_twin(j) = &
(1.0_pReal-sumf-sumftr)*lattice_shearTwin(index_myFamily+i,ph)*& (1.0_pReal-sumf-sumftr)*lattice_shearTwin(index_myFamily+i,ph)*&
state(instance)%twinVolume(j,of)*Ndot0_twin*exp(-StressRatio_r) state(instance)%twinVolume(j,of)*Ndot0_twin*exp(-StressRatio_r)
dgdot_dtautwin(j) = ((gdot_twin(j)*prm%r(f))/tau_twin(j))*StressRatio_r dgdot_dtautwin(j) = ((gdot_twin(j)*prm%r(j))/tau_twin(j))*StressRatio_r
endif endif
!* Plastic velocity gradient for mechanical twinning !* Plastic velocity gradient for mechanical twinning
@ -1677,7 +1660,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
!* Stress ratios !* Stress ratios
if (tau_twin(j) > tol_math_check) then if (tau_twin(j) > tol_math_check) then
StressRatio_r = (state(instance)%threshold_stress_twin(j,of)/& StressRatio_r = (state(instance)%threshold_stress_twin(j,of)/&
tau_twin(j))**prm%r(f) tau_twin(j))**prm%r(j)
!* Shear rates and their derivatives due to twin !* Shear rates and their derivatives due to twin
select case(lattice_structure(ph)) select case(lattice_structure(ph))
case (LATTICE_fcc_ID) case (LATTICE_fcc_ID)
@ -1746,7 +1729,8 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
enddo enddo
enddo enddo
if (el==12) write(6,*) Tstar_v
if (el==12) write(6,*) plasticState(ph)%dotState(:,of),ph,of,el
end subroutine plastic_dislotwin_dotState end subroutine plastic_dislotwin_dotState
@ -1998,7 +1982,7 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
Ndot0_twin=prm%Ndot0_twin(j) Ndot0_twin=prm%Ndot0_twin(j)
end select end select
StressRatio_r = (state(instance)%threshold_stress_twin(j,of)/tau) & StressRatio_r = (state(instance)%threshold_stress_twin(j,of)/tau) &
**prm%r(f) **prm%r(j)
plastic_dislotwin_postResults(c+j) = & plastic_dislotwin_postResults(c+j) = &
(param(instance)%MaxTwinFraction-sumf)*lattice_shearTwin(index_myFamily+i,ph)*& (param(instance)%MaxTwinFraction-sumf)*lattice_shearTwin(index_myFamily+i,ph)*&
state(instance)%twinVolume(j,of)*Ndot0_twin*exp(-StressRatio_r) state(instance)%twinVolume(j,of)*Ndot0_twin*exp(-StressRatio_r)