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

Merge remote-tracking branch 'origin/20-NewStyleDislotwin' into development

This commit is contained in:
Martin Diehl 2019-01-29 16:50:32 +01:00
parent 9b77bdd122 4967ac0132
commit 8764c61589
1 changed files with 92 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

194
src/plastic_dislotwin.f90
View File

@ -46,7 +46,6 @@ module plastic_dislotwin
real(pReal) :: &
mu, &
nu, &
CAtomicVolume, & !< atomic volume in Bugers vector unit
D0, & !< prefactor for self-diffusion coefficient
Qsd, & !< activation energy for dislocation climb
GrainSize, & !<grain size
@ -86,7 +85,8 @@ module plastic_dislotwin
Ndot0_trans, & !< trans nucleation rate [1/m³s] for each trans system
twinsize, & !< twin thickness [m] for each twin system
CLambdaSlip, & !< Adj. parameter for distance between 2 forest dislocations for each slip system
lamellarsize, & !< martensite lamellar thickness [m] for each trans system and instance
atomicVolume, &
lamellarsize, & !< martensite lamellar thickness [m] for each trans system and instance
p, & !< p-exponent in glide velocity
q, & !< q-exponent in glide velocity
r, & !< r-exponent in twin nucleation rate
@ -150,8 +150,8 @@ module plastic_dislotwin
threshold_stress_trans, &
twinVolume, &
martensiteVolume, &
tau_r_twin, & !< stress to bring partial close together for each twin system and instance
tau_r_trans !< stress to bring partial close together for each trans system and instance
tau_r_twin, & !< stress to bring partials close together (twin)
tau_r_trans !< stress to bring partials close together (trans)
end type tDislotwinMicrostructure
!--------------------------------------------------------------------------------------------------
@ -269,6 +269,10 @@ subroutine plastic_dislotwin_init
dst => microstructure(phase_plasticityInstance(p)), &
config => config_phase(p))
prm%aTolRho = config%getFloat('atol_rho', defaultVal=0.0_pReal)
prm%aTolTwinFrac = config%getFloat('atol_twinfrac', defaultVal=0.0_pReal)
prm%aTolTransFrac = config%getFloat('atol_transfrac', defaultVal=0.0_pReal)
! This data is read in already in lattice
prm%mu = lattice_mu(p)
prm%nu = lattice_nu(p)
@ -307,6 +311,9 @@ subroutine plastic_dislotwin_init
defaultVal=[(0.0_pReal, i=1,size(prm%Nslip))]) ! Deprecated
prm%CEdgeDipMinDistance = config%getFloat('cedgedipmindistance')
prm%D0 = config%getFloat('d0')
prm%Qsd = config%getFloat('qsd')
prm%atomicVolume = config%getFloat('catomicvolume') * prm%burgers_slip**3.0_pReal
! expand: family => system
prm%rho0 = math_expand(prm%rho0, prm%Nslip)
@ -318,19 +325,22 @@ subroutine plastic_dislotwin_init
prm%p = math_expand(prm%p, prm%Nslip)
prm%q = math_expand(prm%q, prm%Nslip)
prm%B = math_expand(prm%B, prm%Nslip)
prm%tau_peierls = math_expand(prm%tau_peierls, prm%Nslip)
prm%tau_peierls = math_expand(prm%tau_peierls, prm%Nslip)
prm%atomicVolume = math_expand(prm%atomicVolume,prm%Nslip)
! sanity checks
if (any(prm%rho0 < 0.0_pReal)) extmsg = trim(extmsg)//'rho0 '
if (any(prm%rhoDip0 < 0.0_pReal)) extmsg = trim(extmsg)//'rhoDip0 '
if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//'v0 '
if (any(prm%burgers_slip <= 0.0_pReal)) extmsg = trim(extmsg)//'burgers_slip '
if (any(prm%Qedge <= 0.0_pReal)) extmsg = trim(extmsg)//'Qedge '
if (any(prm%CLambdaSlip <= 0.0_pReal)) extmsg = trim(extmsg)//'CLambdaSlip '
if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//'B '
if (any(prm%tau_peierls < 0.0_pReal)) extmsg = trim(extmsg)//'tau_peierls '
if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//'p '
if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//'q '
if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' D0'
if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' Qsd'
if (any(prm%rho0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho0'
if (any(prm%rhoDip0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoDip0'
if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0'
if (any(prm%burgers_slip <= 0.0_pReal)) extmsg = trim(extmsg)//' burgers_slip'
if (any(prm%Qedge <= 0.0_pReal)) extmsg = trim(extmsg)//' Qedge'
if (any(prm%CLambdaSlip <= 0.0_pReal)) extmsg = trim(extmsg)//' CLambdaSlip'
if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B'
if (any(prm%tau_peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls'
if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p'
if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q'
else slipActive
allocate(prm%burgers_slip(0))
@ -443,67 +453,48 @@ subroutine plastic_dislotwin_init
config%getFloats('interaction_sliptrans'), &
config%getString('lattice_structure'))
if (prm%fccTwinTransNucleation .and. prm%totalNtrans > 12_pInt) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if ntrans is [6,6]
endif
prm%aTolRho = config%getFloat('atol_rho', defaultVal=0.0_pReal)
prm%aTolTwinFrac = config%getFloat('atol_twinfrac', defaultVal=0.0_pReal)
prm%aTolTransFrac = config%getFloat('atol_transfrac', defaultVal=0.0_pReal)
prm%CAtomicVolume = config%getFloat('catomicvolume')
prm%GrainSize = config%getFloat('grainsize')
prm%D0 = config%getFloat('d0')
prm%Qsd = config%getFloat('qsd')
prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! Deprecated
if (config%keyExists('dipoleformationfactor')) call IO_error(1,ext_msg='use /nodipoleformation/')
prm%dipoleformation = .not. config%keyExists('/nodipoleformation/')
prm%sbVelocity = config%getFloat('shearbandvelocity',defaultVal=0.0_pReal)
endif
!--------------------------------------------------------------------------------------------------
! shearband related parameters
prm%sbVelocity = config%getFloat('shearbandvelocity',defaultVal=0.0_pReal)
if (prm%sbVelocity > 0.0_pReal) then
prm%sbResistance = config%getFloat('shearbandresistance')
prm%sbQedge = config%getFloat('qedgepersbsystem')
prm%pShearBand = config%getFloat('p_shearband')
prm%qShearBand = config%getFloat('q_shearband')
prm%sbQedge = config%getFloat('qedgepersbsystem')
prm%pShearBand = config%getFloat('p_shearband')
prm%qShearBand = config%getFloat('q_shearband')
! sanity checks
if (prm%sbResistance < 0.0_pReal) extmsg = trim(extmsg)//' shearbandresistance'
if (prm%sbQedge < 0.0_pReal) extmsg = trim(extmsg)//' qedgepersbsystem'
if (prm%pShearBand <= 0.0_pReal) extmsg = trim(extmsg)//' p_shearband'
if (prm%qShearBand <= 0.0_pReal) extmsg = trim(extmsg)//' q_shearband'
endif
prm%GrainSize = config%getFloat('grainsize')
prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! Deprecated
if (config%keyExists('dipoleformationfactor')) call IO_error(1,ext_msg='use /nodipoleformation/')
prm%dipoleformation = .not. config%keyExists('/nodipoleformation/')
!if (Ndot0PerTwinFamily(f,p) < 0.0_pReal) &
! call IO_error(211_pInt,el=p,ext_msg='ndot0_twin ('//PLASTICITY_DISLOTWIN_label//')')
if (prm%CAtomicVolume <= 0.0_pReal) &
if (any(prm%atomicVolume <= 0.0_pReal)) &
call IO_error(211_pInt,el=p,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOTWIN_label//')')
if (prm%D0 <= 0.0_pReal) &
call IO_error(211_pInt,el=p,ext_msg='D0 ('//PLASTICITY_DISLOTWIN_label//')')
if (prm%Qsd <= 0.0_pReal) &
call IO_error(211_pInt,el=p,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')')
if (prm%totalNtwin > 0_pInt) then
if (dEq0(prm%SFE_0K) .and. &
dEq0(prm%dSFE_dT) .and. &
lattice_structure(p) == LATTICE_fcc_ID) &
call IO_error(211_pInt,el=p,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
if (prm%aTolRho <= 0.0_pReal) &
call IO_error(211_pInt,el=p,ext_msg='aTolRho ('//PLASTICITY_DISLOTWIN_label//')')
if (prm%aTolTwinFrac <= 0.0_pReal) &
call IO_error(211_pInt,el=p,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOTWIN_label//')')
endif
if (prm%totalNtrans > 0_pInt) then
if (dEq0(prm%SFE_0K) .and. &
dEq0(prm%dSFE_dT) .and. &
lattice_structure(p) == LATTICE_fcc_ID) &
call IO_error(211_pInt,el=p,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
if (prm%aTolTransFrac <= 0.0_pReal) &
call IO_error(211_pInt,el=p,ext_msg='aTolTransFrac ('//PLASTICITY_DISLOTWIN_label//')')
endif
!if (prm%sbResistance < 0.0_pReal) &
! call IO_error(211_pInt,el=p,ext_msg='sbResistance ('//PLASTICITY_DISLOTWIN_label//')')
!if (prm%sbVelocity < 0.0_pReal) &
! call IO_error(211_pInt,el=p,ext_msg='sbVelocity ('//PLASTICITY_DISLOTWIN_label//')')
!if (prm%sbVelocity > 0.0_pReal .and. &
! prm%pShearBand <= 0.0_pReal) &
! call IO_error(211_pInt,el=p,ext_msg='pShearBand ('//PLASTICITY_DISLOTWIN_label//')')
if (prm%sbVelocity > 0.0_pReal .and. &
prm%qShearBand <= 0.0_pReal) &
call IO_error(211_pInt,el=p,ext_msg='qShearBand ('//PLASTICITY_DISLOTWIN_label//')')
outputs = config%getStrings('(output)', defaultVal=emptyStringArray)
allocate(prm%outputID(0))
@ -598,13 +589,13 @@ subroutine plastic_dislotwin_init
plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:)
startIndex = endIndex + 1_pInt
endIndex=endIndex+prm%totalNtwin
endIndex = endIndex + prm%totalNtwin
stt%twinFraction=>plasticState(p)%state(startIndex:endIndex,:)
dot%twinFraction=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolTwinFrac
startIndex = endIndex + 1_pInt
endIndex=endIndex+prm%totalNtrans
endIndex = endIndex + prm%totalNtrans
stt%strainTransFraction=>plasticState(p)%state(startIndex:endIndex,:)
dot%strainTransFraction=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolTransFrac
@ -618,13 +609,13 @@ subroutine plastic_dislotwin_init
allocate(dst%invLambdaTwin (prm%totalNtwin, NipcMyPhase),source=0.0_pReal)
allocate(dst%mfp_twin (prm%totalNtwin, NipcMyPhase),source=0.0_pReal)
allocate(dst%threshold_stress_twin (prm%totalNtwin, NipcMyPhase),source=0.0_pReal)
allocate(dst%tau_r_twin (prm%totalNtwin, NipcMyPhase),source=0.0_pReal) !* equilibrium separation of partial dislocations (twin)
allocate(dst%tau_r_twin (prm%totalNtwin, NipcMyPhase),source=0.0_pReal)
allocate(dst%twinVolume (prm%totalNtwin, NipcMyPhase),source=0.0_pReal)
allocate(dst%invLambdaTrans (prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
allocate(dst%mfp_trans (prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
allocate(dst%threshold_stress_trans(prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
allocate(dst%tau_r_trans (prm%totalNtrans,NipcMyPhase),source=0.0_pReal) !* equilibrium separation of partial dislocations (trans)
allocate(dst%tau_r_trans (prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
allocate(dst%martensiteVolume (prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
@ -814,8 +805,9 @@ subroutine plastic_dislotwin_dotState(Mp,Temperature,instance,of)
tol_math_check, &
dEq0
use math, only: &
math_clip, &
math_mul33xx33, &
pi
PI
use material, only: &
plasticState
@ -830,11 +822,13 @@ subroutine plastic_dislotwin_dotState(Mp,Temperature,instance,of)
integer(pInt) :: i
real(pReal) :: f_unrotated,&
EdgeDipMinDistance,AtomicVolume,VacancyDiffusion,&
VacancyDiffusion,&
EdgeDipDistance, ClimbVelocity,DotRhoEdgeDipClimb,DotRhoEdgeDipAnnihilation, &
DotRhoDipFormation,DotRhoMultiplication,DotRhoEdgeEdgeAnnihilation, &
DotRhoDipFormation,DotRhoEdgeEdgeAnnihilation, &
tau
real(pReal), dimension(plasticState(instance)%Nslip) :: &
EdgeDipMinDistance, &
DotRhoMultiplication, &
gdot_slip
real(pReal), dimension(plasticState(instance)%Ntwin) :: &
gdot_twin
@ -847,55 +841,51 @@ subroutine plastic_dislotwin_dotState(Mp,Temperature,instance,of)
f_unrotated = 1.0_pReal &
- sum(stt%twinFraction(1_pInt:prm%totalNtwin,of)) &
- sum(stt%strainTransFraction(1_pInt:prm%totalNtrans,of))
VacancyDiffusion = prm%D0*exp(-prm%Qsd/(kB*Temperature))
call kinetics_slip(Mp,temperature,instance,of,gdot_slip)
dot%accshear_slip(:,of) = abs(gdot_slip)
DotRhoMultiplication = abs(gdot_slip)/(prm%burgers_slip*dst%mfp_slip(:,of))
EdgeDipMinDistance = prm%CEdgeDipMinDistance*prm%burgers_slip
slipState: do i = 1_pInt, prm%totalNslip
tau = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i))
DotRhoMultiplication = abs(gdot_slip(i))/(prm%burgers_slip(i)*dst%mfp_slip(i,of))
EdgeDipMinDistance = prm%CEdgeDipMinDistance*prm%burgers_slip(i)
significantSlipStress2: if (dEq0(tau)) then
significantSlipStress: if (dEq0(tau)) then
DotRhoDipFormation = 0.0_pReal
else significantSlipStress2
EdgeDipDistance = (3.0_pReal*prm%mu*prm%burgers_slip(i))/(16.0_pReal*PI*abs(tau))
if (EdgeDipDistance>dst%mfp_slip(i,of)) EdgeDipDistance = dst%mfp_slip(i,of)
if (EdgeDipDistance<EdgeDipMinDistance) EdgeDipDistance = EdgeDipMinDistance
DotRhoEdgeDipClimb = 0.0_pReal
else significantSlipStress
EdgeDipDistance = 3.0_pReal*prm%mu*prm%burgers_slip(i)/(16.0_pReal*PI*abs(tau))
EdgeDipDistance = math_clip(EdgeDipDistance, right = dst%mfp_slip(i,of))
EdgeDipDistance = math_clip(EdgeDipDistance, left = EdgeDipMinDistance(i))
if (prm%dipoleFormation) then
DotRhoDipFormation = ((2.0_pReal*(EdgeDipDistance-EdgeDipMinDistance))/prm%burgers_slip(i)) &
DotRhoDipFormation = 2.0_pReal*(EdgeDipDistance-EdgeDipMinDistance(i))/prm%burgers_slip(i) &
* stt%rhoEdge(i,of)*abs(gdot_slip(i))
else
DotRhoDipFormation = 0.0_pReal
endif
endif significantSlipStress2
!* Spontaneous annihilation of 2 single edge dislocations
DotRhoEdgeEdgeAnnihilation = 2.0_pReal*EdgeDipMinDistance/prm%burgers_slip(i) &
* stt%rhoEdge(i,of)*abs(gdot_slip(i))
!* Spontaneous annihilation of a single edge dislocation with a dipole constituent
DotRhoEdgeDipAnnihilation = 2.0_pReal*EdgeDipMinDistance/prm%burgers_slip(i) &
* stt%rhoEdgeDip(i,of)*abs(gdot_slip(i))
!* Dislocation dipole climb
AtomicVolume = prm%CAtomicVolume*prm%burgers_slip(i)**(3.0_pReal) ! no need to calculate this over and over again
VacancyDiffusion = prm%D0*exp(-prm%Qsd/(kB*Temperature))
if (dEq0(tau)) then
DotRhoEdgeDipClimb = 0.0_pReal
else
if (dEq0(EdgeDipDistance-EdgeDipMinDistance)) then
if (dEq0(EdgeDipDistance-EdgeDipMinDistance(i))) then
DotRhoEdgeDipClimb = 0.0_pReal
else
ClimbVelocity = 3.0_pReal*prm%mu*VacancyDiffusion*AtomicVolume/ &
(2.0_pReal*pi*kB*Temperature*(EdgeDipDistance+EdgeDipMinDistance))
DotRhoEdgeDipClimb = 4.0_pReal*ClimbVelocity*stt%rhoEdgeDip(i,of)/ &
(EdgeDipDistance-EdgeDipMinDistance)
ClimbVelocity = 3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume(i) &
/ (2.0_pReal*PI*kB*Temperature*(EdgeDipDistance+EdgeDipMinDistance(i)))
DotRhoEdgeDipClimb = 4.0_pReal*ClimbVelocity*stt%rhoEdgeDip(i,of) &
/ (EdgeDipDistance-EdgeDipMinDistance(i))
endif
endif
dot%rhoEdge(i,of) = DotRhoMultiplication-DotRhoDipFormation-DotRhoEdgeEdgeAnnihilation
endif significantSlipStress
!* Spontaneous annihilation of 2 single edge dislocations
DotRhoEdgeEdgeAnnihilation = 2.0_pReal*EdgeDipMinDistance(i)/prm%burgers_slip(i) &
* stt%rhoEdge(i,of)*abs(gdot_slip(i))
!* Spontaneous annihilation of a single edge dislocation with a dipole constituent
DotRhoEdgeDipAnnihilation = 2.0_pReal*EdgeDipMinDistance(i)/prm%burgers_slip(i) &
* stt%rhoEdgeDip(i,of)*abs(gdot_slip(i))
dot%rhoEdge(i,of) = DotRhoMultiplication(i)-DotRhoDipFormation-DotRhoEdgeEdgeAnnihilation
dot%rhoEdgeDip(i,of) = DotRhoDipFormation-DotRhoEdgeDipAnnihilation-DotRhoEdgeDipClimb
dot%accshear_slip(i,of) = abs(gdot_slip(i))
enddo slipState
call kinetics_twin(Mp,temperature,gdot_slip,instance,of,gdot_twin)
@ -918,10 +908,10 @@ subroutine plastic_dislotwin_dependentState(temperature,instance,of)
implicit none
integer(pInt), intent(in) :: &
instance, & !< component-ID of integration point
instance, &
of
real(pReal), intent(in) :: &
temperature !< temperature at IP
temperature
integer(pInt) :: &
i
@ -940,7 +930,7 @@ subroutine plastic_dislotwin_dependentState(temperature,instance,of)
sumf_twin = sum(stt%twinFraction(1:prm%totalNtwin,of))
sumf_trans = sum(stt%strainTransFraction(1:prm%totalNtrans,of))
sfe = prm%SFE_0K + prm%dSFE_dT * Temperature
SFE = prm%SFE_0K + prm%dSFE_dT * Temperature
!* rescaled volume fraction for topology
fOverStacksize = stt%twinFraction(1_pInt:prm%totalNtwin,of)/prm%twinsize !ToDo: this is per system
@ -999,11 +989,11 @@ subroutine plastic_dislotwin_dependentState(temperature,instance,of)
!* threshold stress for growing twin/martensite
if(prm%totalNtwin == prm%totalNslip) &
dst%threshold_stress_twin(:,of) = prm%Cthresholdtwin* &
(sfe/(3.0_pReal*prm%burgers_twin)+ 3.0_pReal*prm%burgers_twin*prm%mu/ &
(SFE/(3.0_pReal*prm%burgers_twin)+ 3.0_pReal*prm%burgers_twin*prm%mu/ &
(prm%L0_twin*prm%burgers_slip)) ! slip burgers here correct?
if(prm%totalNtrans == prm%totalNslip) &
dst%threshold_stress_trans(:,of) = prm%Cthresholdtrans* &
(sfe/(3.0_pReal*prm%burgers_trans) + 3.0_pReal*prm%burgers_trans*prm%mu/&
(SFE/(3.0_pReal*prm%burgers_trans) + 3.0_pReal*prm%burgers_trans*prm%mu/&
(prm%L0_trans*prm%burgers_slip) + prm%transStackHeight*prm%deltaG/ (3.0_pReal*prm%burgers_trans) )