started to introduce non-schmid behavior and disabled shear banding in dislokmc

This commit is contained in:
Martin Diehl 2014-10-27 09:19:36 +00:00
parent c5176f3bf7
commit bec0af7b06
1 changed files with 142 additions and 175 deletions

View File

@ -59,8 +59,8 @@ module constitutive_dislokmc
constitutive_dislokmc_D0, & !< prefactor for self-diffusion coefficient
constitutive_dislokmc_Qsd, & !< activation energy for dislocation climb
constitutive_dislokmc_GrainSize, & !< grain size
constitutive_dislokmc_pShearBand, & !< p-exponent in shearband velocity
constitutive_dislokmc_qShearBand, & !< q-exponent in shearband velocity
!constitutive_dislokmc_pShearBand, & !< p-exponent in shearband velocity
!constitutive_dislokmc_qShearBand, & !< q-exponent in shearband velocity
constitutive_dislokmc_MaxTwinFraction, & !< maximum allowed total twin volume fraction
constitutive_dislokmc_CEdgeDipMinDistance, & !<
constitutive_dislokmc_Cmfptwin, & !<
@ -69,9 +69,9 @@ module constitutive_dislokmc
constitutive_dislokmc_L0, & !< Length of twin nuclei in Burgers vectors
constitutive_dislokmc_xc, & !< critical distance for formation of twin nucleus
constitutive_dislokmc_VcrossSlip, & !< cross slip volume
constitutive_dislokmc_sbResistance, & !< value for shearband resistance (might become an internal state variable at some point)
constitutive_dislokmc_sbVelocity, & !< value for shearband velocity_0
constitutive_dislokmc_sbQedge, & !< value for shearband systems Qedge
!constitutive_dislokmc_sbResistance, & !< value for shearband resistance (might become an internal state variable at some point)
!constitutive_dislokmc_sbVelocity, & !< value for shearband velocity_0
!constitutive_dislokmc_sbQedge, & !< value for shearband systems Qedge
constitutive_dislokmc_SFE_0K, & !< stacking fault energy at zero K
constitutive_dislokmc_dSFE_dT, & !< temperature dependance of stacking fault energy
constitutive_dislokmc_dipoleFormationFactor, & !< scaling factor for dipole formation: 0: off, 1: on. other values not useful
@ -110,15 +110,14 @@ module constitutive_dislokmc
constitutive_dislokmc_uPerSlipFamily, & !< u-exponent in glide velocity
constitutive_dislokmc_vPerSlipFamily, & !< v-exponent in glide velocity
constitutive_dislokmc_sPerSlipFamily, & !< self-hardening in glide velocity
constitutive_dislokmc_rPerTwinFamily !< r-exponent in twin nucleation rate
constitutive_dislokmc_rPerTwinFamily, & !< r-exponent in twin nucleation rate
constitutive_dislokmc_nonSchmidCoeff !< non-Schmid coefficients (bcc)
real(pReal), dimension(:,:,:), allocatable, private :: &
constitutive_dislokmc_interactionMatrix_SlipSlip, & !< interaction matrix of the different slip systems for each instance
constitutive_dislokmc_interactionMatrix_SlipTwin, & !< interaction matrix of slip systems with twin systems for each instance
constitutive_dislokmc_interactionMatrix_TwinSlip, & !< interaction matrix of twin systems with slip systems for each instance
constitutive_dislokmc_interactionMatrix_TwinTwin, & !< interaction matrix of the different twin systems for each instance
constitutive_dislokmc_forestProjectionEdge !< matrix of forest projections of edge dislocations for each instance
real(pReal), dimension(:,:,:,:,:), allocatable, private :: &
constitutive_dislokmc_sbSv
enum, bind(c)
enumerator :: undefined_ID, &
@ -136,11 +135,11 @@ module constitutive_dislokmc
accumulated_shear_twin_ID, &
mfp_twin_ID, &
resolved_stress_twin_ID, &
threshold_stress_twin_ID, &
resolved_stress_shearband_ID, &
shear_rate_shearband_ID, &
sb_eigenvalues_ID, &
sb_eigenvectors_ID
threshold_stress_twin_ID
!resolved_stress_shearband_ID, &
!shear_rate_shearband_ID, &
!sb_eigenvalues_ID, &
!sb_eigenvectors_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
constitutive_dislokmc_outputID !< ID of each post result output
@ -215,7 +214,7 @@ subroutine constitutive_dislokmc_init(fileUnit)
integer(pInt) :: maxNinstance,mySize=0_pInt,phase,maxTotalNslip,maxTotalNtwin,&
f,instance,j,k,l,m,n,o,p,q,r,s,ns,nt, &
Nchunks_SlipSlip, Nchunks_SlipTwin, Nchunks_TwinSlip, Nchunks_TwinTwin, &
Nchunks_SlipFamilies, Nchunks_TwinFamilies, &
Nchunks_SlipFamilies, Nchunks_TwinFamilies, Nchunks_nonSchmid, &
index_myFamily, index_otherFamily
integer(pInt) :: sizeState, sizeDotState
integer(pInt) :: NofMyPhase
@ -251,8 +250,8 @@ subroutine constitutive_dislokmc_init(fileUnit)
allocate(constitutive_dislokmc_D0(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_Qsd(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_GrainSize(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_pShearBand(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_qShearBand(maxNinstance), source=0.0_pReal)
!allocate(constitutive_dislokmc_pShearBand(maxNinstance), source=0.0_pReal)
!allocate(constitutive_dislokmc_qShearBand(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_MaxTwinFraction(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_CEdgeDipMinDistance(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_Cmfptwin(maxNinstance), source=0.0_pReal)
@ -263,9 +262,9 @@ subroutine constitutive_dislokmc_init(fileUnit)
allocate(constitutive_dislokmc_VcrossSlip(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_aTolRho(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_aTolTwinFrac(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_sbResistance(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_sbVelocity(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_sbQedge(maxNinstance), source=0.0_pReal)
!allocate(constitutive_dislokmc_sbResistance(maxNinstance), source=0.0_pReal)
!allocate(constitutive_dislokmc_sbVelocity(maxNinstance), source=0.0_pReal)
!allocate(constitutive_dislokmc_sbQedge(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_SFE_0K(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_dSFE_dT(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislokmc_dipoleFormationFactor(maxNinstance), source=1.0_pReal) !should be on by default
@ -301,8 +300,8 @@ subroutine constitutive_dislokmc_init(fileUnit)
source=0.0_pReal)
allocate(constitutive_dislokmc_interaction_TwinTwin(lattice_maxNinteraction,maxNinstance), &
source=0.0_pReal)
allocate(constitutive_dislokmc_sbSv(6,6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), &
source=0.0_pReal)
allocate(constitutive_dislokmc_nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstance), &
source=0.0_pReal)
rewind(fileUnit)
@ -327,6 +326,7 @@ subroutine constitutive_dislokmc_init(fileUnit)
Nchunks_SlipTwin = maxval(lattice_interactionSlipTwin(:,:,phase))
Nchunks_TwinSlip = maxval(lattice_interactionTwinSlip(:,:,phase))
Nchunks_TwinTwin = maxval(lattice_interactionTwinTwin(:,:,phase))
Nchunks_nonSchmid = lattice_NnonSchmid(phase)
if(allocated(tempPerSlip)) deallocate(tempPerSlip)
if(allocated(tempPerTwin)) deallocate(tempPerTwin)
allocate(tempPerSlip(Nchunks_SlipFamilies))
@ -411,26 +411,26 @@ subroutine constitutive_dislokmc_init(fileUnit)
constitutive_dislokmc_outputID(constitutive_dislokmc_Noutput(instance),instance) = threshold_stress_twin_ID
constitutive_dislokmc_output(constitutive_dislokmc_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,positions,2_pInt))
case ('resolved_stress_shearband')
constitutive_dislokmc_Noutput(instance) = constitutive_dislokmc_Noutput(instance) + 1_pInt
constitutive_dislokmc_outputID(constitutive_dislokmc_Noutput(instance),instance) = resolved_stress_shearband_ID
constitutive_dislokmc_output(constitutive_dislokmc_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,positions,2_pInt))
case ('shear_rate_shearband')
constitutive_dislokmc_Noutput(instance) = constitutive_dislokmc_Noutput(instance) + 1_pInt
constitutive_dislokmc_outputID(constitutive_dislokmc_Noutput(instance),instance) = shear_rate_shearband_ID
constitutive_dislokmc_output(constitutive_dislokmc_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,positions,2_pInt))
case ('sb_eigenvalues')
constitutive_dislokmc_Noutput(instance) = constitutive_dislokmc_Noutput(instance) + 1_pInt
constitutive_dislokmc_outputID(constitutive_dislokmc_Noutput(instance),instance) = sb_eigenvalues_ID
constitutive_dislokmc_output(constitutive_dislokmc_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,positions,2_pInt))
case ('sb_eigenvectors')
constitutive_dislokmc_Noutput(instance) = constitutive_dislokmc_Noutput(instance) + 1_pInt
constitutive_dislokmc_outputID(constitutive_dislokmc_Noutput(instance),instance) = sb_eigenvectors_ID
constitutive_dislokmc_output(constitutive_dislokmc_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,positions,2_pInt))
! case ('resolved_stress_shearband')
! constitutive_dislokmc_Noutput(instance) = constitutive_dislokmc_Noutput(instance) + 1_pInt
! constitutive_dislokmc_outputID(constitutive_dislokmc_Noutput(instance),instance) = resolved_stress_shearband_ID
! constitutive_dislokmc_output(constitutive_dislokmc_Noutput(instance),instance) = &
! IO_lc(IO_stringValue(line,positions,2_pInt))
! case ('shear_rate_shearband')
! constitutive_dislokmc_Noutput(instance) = constitutive_dislokmc_Noutput(instance) + 1_pInt
! constitutive_dislokmc_outputID(constitutive_dislokmc_Noutput(instance),instance) = shear_rate_shearband_ID
! constitutive_dislokmc_output(constitutive_dislokmc_Noutput(instance),instance) = &
! IO_lc(IO_stringValue(line,positions,2_pInt))
! case ('sb_eigenvalues')
! constitutive_dislokmc_Noutput(instance) = constitutive_dislokmc_Noutput(instance) + 1_pInt
! constitutive_dislokmc_outputID(constitutive_dislokmc_Noutput(instance),instance) = sb_eigenvalues_ID
! constitutive_dislokmc_output(constitutive_dislokmc_Noutput(instance),instance) = &
! IO_lc(IO_stringValue(line,positions,2_pInt))
! case ('sb_eigenvectors')
! constitutive_dislokmc_Noutput(instance) = constitutive_dislokmc_Noutput(instance) + 1_pInt
! constitutive_dislokmc_outputID(constitutive_dislokmc_Noutput(instance),instance) = sb_eigenvectors_ID
! constitutive_dislokmc_output(constitutive_dislokmc_Noutput(instance),instance) = &
! IO_lc(IO_stringValue(line,positions,2_pInt))
end select
!--------------------------------------------------------------------------------------------------
! parameters depending on number of slip system families
@ -529,16 +529,22 @@ subroutine constitutive_dislokmc_init(fileUnit)
do j = 1_pInt, Nchunks_TwinTwin
constitutive_dislokmc_interaction_TwinTwin(j,instance) = IO_floatValue(line,positions,1_pInt+j)
enddo
case ('nonschmid_coefficients')
if (positions(1) < 1_pInt + Nchunks_nonSchmid) &
call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOKMC_label//')')
do j = 1_pInt,Nchunks_nonSchmid
constitutive_dislokmc_nonSchmidCoeff(j,instance) = IO_floatValue(line,positions,1_pInt+j)
enddo
!--------------------------------------------------------------------------------------------------
! parameters independent of number of slip/twin systems
case ('grainsize')
constitutive_dislokmc_GrainSize(instance) = IO_floatValue(line,positions,2_pInt)
case ('maxtwinfraction')
constitutive_dislokmc_MaxTwinFraction(instance) = IO_floatValue(line,positions,2_pInt)
case ('p_shearband')
constitutive_dislokmc_pShearBand(instance) = IO_floatValue(line,positions,2_pInt)
case ('q_shearband')
constitutive_dislokmc_qShearBand(instance) = IO_floatValue(line,positions,2_pInt)
! case ('p_shearband')
! constitutive_dislokmc_pShearBand(instance) = IO_floatValue(line,positions,2_pInt)
! case ('q_shearband')
! constitutive_dislokmc_qShearBand(instance) = IO_floatValue(line,positions,2_pInt)
case ('d0')
constitutive_dislokmc_D0(instance) = IO_floatValue(line,positions,2_pInt)
case ('qsd')
@ -569,12 +575,12 @@ subroutine constitutive_dislokmc_init(fileUnit)
constitutive_dislokmc_dSFE_dT(instance) = IO_floatValue(line,positions,2_pInt)
case ('dipoleformationfactor')
constitutive_dislokmc_dipoleFormationFactor(instance) = IO_floatValue(line,positions,2_pInt)
case ('shearbandresistance')
constitutive_dislokmc_sbResistance(instance) = IO_floatValue(line,positions,2_pInt)
case ('shearbandvelocity')
constitutive_dislokmc_sbVelocity(instance) = IO_floatValue(line,positions,2_pInt)
case ('qedgepersbsystem')
constitutive_dislokmc_sbQedge(instance) = IO_floatValue(line,positions,2_pInt)
!case ('shearbandresistance')
! constitutive_dislokmc_sbResistance(instance) = IO_floatValue(line,positions,2_pInt)
!case ('shearbandvelocity')
! constitutive_dislokmc_sbVelocity(instance) = IO_floatValue(line,positions,2_pInt)
!case ('qedgepersbsystem')
! constitutive_dislokmc_sbQedge(instance) = IO_floatValue(line,positions,2_pInt)
end select
endif; endif
enddo parsingFile
@ -624,19 +630,19 @@ subroutine constitutive_dislokmc_init(fileUnit)
if (constitutive_dislokmc_aTolTwinFrac(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOKMC_label//')')
endif
if (constitutive_dislokmc_sbResistance(instance) < 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='sbResistance ('//PLASTICITY_DISLOKMC_label//')')
if (constitutive_dislokmc_sbVelocity(instance) < 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='sbVelocity ('//PLASTICITY_DISLOKMC_label//')')
if (constitutive_dislokmc_sbVelocity(instance) > 0.0_pReal .and. &
constitutive_dislokmc_pShearBand(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='pShearBand ('//PLASTICITY_DISLOKMC_label//')')
! if (constitutive_dislokmc_sbResistance(instance) < 0.0_pReal) &
! call IO_error(211_pInt,el=instance,ext_msg='sbResistance ('//PLASTICITY_DISLOKMC_label//')')
! if (constitutive_dislokmc_sbVelocity(instance) < 0.0_pReal) &
! call IO_error(211_pInt,el=instance,ext_msg='sbVelocity ('//PLASTICITY_DISLOKMC_label//')')
! if (constitutive_dislokmc_sbVelocity(instance) > 0.0_pReal .and. &
! constitutive_dislokmc_pShearBand(instance) <= 0.0_pReal) &
! call IO_error(211_pInt,el=instance,ext_msg='pShearBand ('//PLASTICITY_DISLOKMC_label//')')
if (constitutive_dislokmc_dipoleFormationFactor(instance) /= 0.0_pReal .and. &
constitutive_dislokmc_dipoleFormationFactor(instance) /= 1.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='dipoleFormationFactor ('//PLASTICITY_DISLOKMC_label//')')
if (constitutive_dislokmc_sbVelocity(instance) > 0.0_pReal .and. &
constitutive_dislokmc_qShearBand(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='qShearBand ('//PLASTICITY_DISLOKMC_label//')')
! if (constitutive_dislokmc_sbVelocity(instance) > 0.0_pReal .and. &
! constitutive_dislokmc_qShearBand(instance) <= 0.0_pReal) &
! call IO_error(211_pInt,el=instance,ext_msg='qShearBand ('//PLASTICITY_DISLOKMC_label//')')
!--------------------------------------------------------------------------------------------------
! Determine total number of active slip or twin systems
@ -709,14 +715,14 @@ subroutine constitutive_dislokmc_init(fileUnit)
threshold_stress_twin_ID &
)
mySize = nt
case(resolved_stress_shearband_ID, &
shear_rate_shearband_ID &
)
mySize = 6_pInt
case(sb_eigenvalues_ID)
mySize = 3_pInt
case(sb_eigenvectors_ID)
mySize = 9_pInt
! case(resolved_stress_shearband_ID, &
! shear_rate_shearband_ID &
! )
! mySize = 6_pInt
! case(sb_eigenvalues_ID)
! mySize = 3_pInt
! case(sb_eigenvectors_ID)
! mySize = 9_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
@ -1216,6 +1222,7 @@ subroutine constitutive_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperatu
lattice_maxNtwinFamily, &
lattice_NslipSystem, &
lattice_NtwinSystem, &
lattice_NnonSchmid, &
lattice_shearTwin, &
lattice_structure, &
lattice_fcc_twinNucleationSlipPair, &
@ -1231,34 +1238,15 @@ subroutine constitutive_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperatu
integer(pInt) :: instance,ph,of,ns,nt,f,i,j,k,l,m,n,index_myFamily,s1,s2
real(pReal) :: sumf,StressRatio_p,StressRatio_pminus1,StressRatio_r,BoltzmannRatio,DotGamma0,Ndot0, &
StressRatio_u,StressRatio_uminus1
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333
real(pReal), dimension(3,3,2) :: &
nonSchmid_tensor
real(pReal), dimension(3,3,3,3) :: &
dLp_dTstar3333
real(pReal), dimension(constitutive_dislokmc_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip,dgdot_dtauslip,tau_slip,vel_slip,dvel_slip
gdot_slip_pos,dgdot_dtauslip,tau_slip_pos,tau_slip_neg,vel_slip,dvel_slip
real(pReal), dimension(constitutive_dislokmc_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_twin,dgdot_dtautwin,tau_twin
real(pReal), dimension(6) :: gdot_sb,dgdot_dtausb,tau_sb
real(pReal), dimension(3,3) :: eigVectors, sb_Smatrix
real(pReal), dimension(3) :: eigValues, sb_s, sb_m
real(pReal), dimension(3,6), parameter :: &
sb_sComposition = &
reshape(real([&
1, 0, 1, &
1, 0,-1, &
1, 1, 0, &
1,-1, 0, &
0, 1, 1, &
0, 1,-1 &
],pReal),[ 3,6]), &
sb_mComposition = &
reshape(real([&
1, 0,-1, &
1, 0,+1, &
1,-1, 0, &
1, 1, 0, &
0, 1,-1, &
0, 1, 1 &
],pReal),[ 3,6])
logical error
!* Shortened notation
of = mappingConstitutive(1,ipc,ip,el)
@ -1275,31 +1263,45 @@ subroutine constitutive_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperatu
dLp_dTstar = 0.0_pReal
!* Dislocation glide part
gdot_slip = 0.0_pReal
gdot_slip_pos = 0.0_pReal
dgdot_dtauslip = 0.0_pReal
j = 0_pInt
slipFamiliesLoop: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family
slipSystemsLoop: do i = 1_pInt,constitutive_dislokmc_Nslip(f,instance)
j = j+1_pInt
j = j+1_pInt
!* Calculation of Lp
!* Resolved shear stress on slip system
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))
!--------------------------------------------------------------------------------------------------
! Calculation of Lp
tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph))
tau_slip_neg(j) = tau_slip_pos(j)
nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1)
if((abs(tau_slip(j))-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then
do k = 1,lattice_NnonSchmid(ph)
tau_slip_pos(j) = tau_slip_pos(j) + constitutive_dislokmc_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph))
tau_slip_neg(j) = tau_slip_neg(j) + constitutive_dislokmc_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + constitutive_dislokmc_nonSchmidCoeff(k,instance)*&
lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,ph)
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + constitutive_dislokmc_nonSchmidCoeff(k,instance)*&
lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph)
enddo
if((abs(tau_slip_pos(j))-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then
!* Stress ratios
StressRatio_p = ((abs(tau_slip(j))- plasticState(ph)%state(6*ns+4*nt+j, of))/&
StressRatio_p = ((abs(tau_slip_pos(j))- plasticState(ph)%state(6*ns+4*nt+j, of))/&
(constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
**constitutive_dislokmc_pPerSlipFamily(f,instance)
StressRatio_pminus1 = ((abs(tau_slip(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
StressRatio_pminus1 = ((abs(tau_slip_pos(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
(constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
**(constitutive_dislokmc_pPerSlipFamily(f,instance)-1.0_pReal)
StressRatio_u = ((abs(tau_slip(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
StressRatio_u = ((abs(tau_slip_pos(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
(constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
**constitutive_dislokmc_uPerSlipFamily(f,instance)
StressRatio_uminus1 = ((abs(tau_slip(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
StressRatio_uminus1 = ((abs(tau_slip_pos(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
(constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
**(constitutive_dislokmc_uPerSlipFamily(f,instance)-1.0_pReal)
@ -1315,9 +1317,9 @@ subroutine constitutive_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperatu
* (1.0_pReal-constitutive_dislokmc_sPerSlipFamily(f,instance) &
* exp(-BoltzmannRatio*(1-StressRatio_p) ** constitutive_dislokmc_qPerSlipFamily(f,instance)))
gdot_slip(j) = (1.0_pReal - sumf) * DotGamma0 &
gdot_slip_pos(j) = (1.0_pReal - sumf) * DotGamma0 &
* StressRatio_u * vel_slip(j) &
* sign(1.0_pReal,tau_slip(j))
* sign(1.0_pReal,tau_slip_pos(j))
!* Derivatives of shear rates
@ -1340,7 +1342,7 @@ subroutine constitutive_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperatu
endif
!* Plastic velocity gradient for dislocation glide
Lp = Lp + gdot_slip(j)*lattice_Sslip(:,:,1,index_myFamily+i,ph)
Lp = Lp + gdot_slip_pos(j)*lattice_Sslip(:,:,1,index_myFamily+i,ph)
!* Calculation of the tangent of Lp
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
@ -1352,6 +1354,7 @@ subroutine constitutive_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperatu
enddo slipFamiliesLoop
!* Shear banding (shearband) part
#ifdef Shearband
if(constitutive_dislokmc_sbVelocity(instance) /= 0.0_pReal .and. &
constitutive_dislokmc_sbResistance(instance) /= 0.0_pReal) then
gdot_sb = 0.0_pReal
@ -1405,6 +1408,7 @@ subroutine constitutive_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperatu
sb_Smatrix(m,n)
enddo
end if
#endif
!* Mechanical twinning part
gdot_twin = 0.0_pReal
@ -1429,8 +1433,8 @@ subroutine constitutive_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperatu
s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i)
s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i)
if (tau_twin(j) < constitutive_dislokmc_tau_r(j,instance)) then
Ndot0=(abs(gdot_slip(s1))*(plasticState(ph)%state(s2,of)+plasticState(ph)%state(ns+s2, of))+&
abs(gdot_slip(s2))*(plasticState(ph)%state(s1, of)+plasticState(ph)%state(ns+s1, of)))/&
Ndot0=(abs(gdot_slip_pos(s1))*(plasticState(ph)%state(s2,of)+plasticState(ph)%state(ns+s2, of))+&
abs(gdot_slip_pos(s2))*(plasticState(ph)%state(s1, of)+plasticState(ph)%state(ns+s1, of)))/&
(constitutive_dislokmc_L0(instance)*constitutive_dislokmc_burgersPerSlipSystem(j,instance))*&
(1.0_pReal-exp(-constitutive_dislokmc_VcrossSlip(instance)/(kB*Temperature)*&
(constitutive_dislokmc_tau_r(j,instance)-tau_twin(j))))
@ -1510,7 +1514,7 @@ subroutine constitutive_dislokmc_dotState(Tstar_v,Temperature,ipc,ip,el)
real(pReal) :: sumf,StressRatio_p,StressRatio_pminus1,BoltzmannRatio,DotGamma0,StressRatio_u,StressRatio_uminus1,&
EdgeDipMinDistance,AtomicVolume,VacancyDiffusion,StressRatio_r,Ndot0
real(pReal), dimension(constitutive_dislokmc_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip,tau_slip,DotRhoMultiplication,EdgeDipDistance,DotRhoEdgeEdgeAnnihilation,DotRhoEdgeDipAnnihilation,&
gdot_slip_pos,tau_slip_pos,DotRhoMultiplication,EdgeDipDistance,DotRhoEdgeEdgeAnnihilation,DotRhoEdgeDipAnnihilation,&
ClimbVelocity,DotRhoEdgeDipClimb,DotRhoDipFormation
real(pReal), dimension(constitutive_dislokmc_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
tau_twin
@ -1527,7 +1531,7 @@ subroutine constitutive_dislokmc_dotState(Tstar_v,Temperature,ipc,ip,el)
plasticState(ph)%dotState(:,of) = 0.0_pReal
!* Dislocation density evolution
gdot_slip = 0.0_pReal
gdot_slip_pos = 0.0_pReal
j = 0_pInt
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family
@ -1536,21 +1540,21 @@ subroutine constitutive_dislokmc_dotState(Tstar_v,Temperature,ipc,ip,el)
!* Resolved shear stress on slip system
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))
tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))
if((abs(tau_slip(j))-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then
if((abs(tau_slip_pos(j))-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then
!* Stress ratios
StressRatio_p = ((abs(tau_slip(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
StressRatio_p = ((abs(tau_slip_pos(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
(constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
**constitutive_dislokmc_pPerSlipFamily(f,instance)
StressRatio_pminus1 = ((abs(tau_slip(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
StressRatio_pminus1 = ((abs(tau_slip_pos(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
(constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
**(constitutive_dislokmc_pPerSlipFamily(f,instance)-1.0_pReal)
StressRatio_u = ((abs(tau_slip(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
StressRatio_u = ((abs(tau_slip_pos(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
(constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
**constitutive_dislokmc_uPerSlipFamily(f,instance)
StressRatio_uminus1 = ((abs(tau_slip(j))-plasticState(ph)%state(6*ns+4*nt+j,of))/&
StressRatio_uminus1 = ((abs(tau_slip_pos(j))-plasticState(ph)%state(6*ns+4*nt+j,of))/&
(constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
**(constitutive_dislokmc_uPerSlipFamily(f,instance)-1.0_pReal)
@ -1561,50 +1565,50 @@ subroutine constitutive_dislokmc_dotState(Tstar_v,Temperature,ipc,ip,el)
plasticState(ph)%state(j, of)*constitutive_dislokmc_burgersPerSlipSystem(j,instance)*&
constitutive_dislokmc_v0PerSlipSystem(j,instance)
!* Shear rates due to slip
gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)** &
constitutive_dislokmc_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau_slip(j)) &
gdot_slip_pos(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)** &
constitutive_dislokmc_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau_slip_pos(j)) &
* (1_pInt-constitutive_dislokmc_sPerSlipFamily(f,instance) &
* exp(-BoltzmannRatio*(1_pInt-StressRatio_p) ** constitutive_dislokmc_qPerSlipFamily(f,instance))) &
* StressRatio_u
endif
!* Multiplication
DotRhoMultiplication(j) = abs(gdot_slip(j))/&
DotRhoMultiplication(j) = abs(gdot_slip_pos(j))/&
(constitutive_dislokmc_burgersPerSlipSystem(j,instance)* &
plasticState(ph)%state(5*ns+3*nt+j, of))
!* Dipole formation
EdgeDipMinDistance = &
constitutive_dislokmc_CEdgeDipMinDistance(instance)*constitutive_dislokmc_burgersPerSlipSystem(j,instance)
if (tau_slip(j) == 0.0_pReal) then
if (tau_slip_pos(j) == 0.0_pReal) then
DotRhoDipFormation(j) = 0.0_pReal
else
EdgeDipDistance(j) = &
(3.0_pReal*lattice_mu(ph)*constitutive_dislokmc_burgersPerSlipSystem(j,instance))/&
(16.0_pReal*pi*abs(tau_slip(j)))
(16.0_pReal*pi*abs(tau_slip_pos(j)))
if (EdgeDipDistance(j)>plasticState(ph)%state(5*ns+3*nt+j, of)) EdgeDipDistance(j)=plasticState(ph)%state(5*ns+3*nt+j, of)
if (EdgeDipDistance(j)<EdgeDipMinDistance) EdgeDipDistance(j)=EdgeDipMinDistance
DotRhoDipFormation(j) = &
((2.0_pReal*EdgeDipDistance(j))/constitutive_dislokmc_burgersPerSlipSystem(j,instance))*&
plasticState(ph)%state(j, of)*abs(gdot_slip(j))*constitutive_dislokmc_dipoleFormationFactor(instance)
plasticState(ph)%state(j, of)*abs(gdot_slip_pos(j))*constitutive_dislokmc_dipoleFormationFactor(instance)
endif
!* Spontaneous annihilation of 2 single edge dislocations
DotRhoEdgeEdgeAnnihilation(j) = &
((2.0_pReal*EdgeDipMinDistance)/constitutive_dislokmc_burgersPerSlipSystem(j,instance))*&
plasticState(ph)%state(j, of)*abs(gdot_slip(j))
plasticState(ph)%state(j, of)*abs(gdot_slip_pos(j))
!* Spontaneous annihilation of a single edge dislocation with a dipole constituent
DotRhoEdgeDipAnnihilation(j) = &
((2.0_pReal*EdgeDipMinDistance)/constitutive_dislokmc_burgersPerSlipSystem(j,instance))*&
plasticState(ph)%state(ns+j, of)*abs(gdot_slip(j))
plasticState(ph)%state(ns+j, of)*abs(gdot_slip_pos(j))
!* Dislocation dipole climb
AtomicVolume = &
constitutive_dislokmc_CAtomicVolume(instance)*constitutive_dislokmc_burgersPerSlipSystem(j,instance)**(3.0_pReal)
VacancyDiffusion = &
constitutive_dislokmc_D0(instance)*exp(-constitutive_dislokmc_Qsd(instance)/(kB*Temperature))
if (tau_slip(j) == 0.0_pReal) then
if (tau_slip_pos(j) == 0.0_pReal) then
DotRhoEdgeDipClimb(j) = 0.0_pReal
else
ClimbVelocity(j) = &
@ -1623,7 +1627,7 @@ subroutine constitutive_dislokmc_dotState(Tstar_v,Temperature,ipc,ip,el)
DotRhoDipFormation(j)-DotRhoEdgeDipAnnihilation(j)-DotRhoEdgeDipClimb(j)
!* Dotstate for accumulated shear due to slip
plasticState(ph)%dotState(2_pInt*ns+j, of) = gdot_slip(j)
plasticState(ph)%dotState(2_pInt*ns+j, of) = gdot_slip_pos(j)
enddo
enddo
@ -1647,8 +1651,8 @@ subroutine constitutive_dislokmc_dotState(Tstar_v,Temperature,ipc,ip,el)
s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i)
s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i)
if (tau_twin(j) < constitutive_dislokmc_tau_r(j,instance)) then
Ndot0=(abs(gdot_slip(s1))*(plasticState(ph)%state(s2, of)+plasticState(ph)%state(ns+s2, of))+&
abs(gdot_slip(s2))*(plasticState(ph)%state(s1, of)+plasticState(ph)%state(ns+s1, of)))/&
Ndot0=(abs(gdot_slip_pos(s1))*(plasticState(ph)%state(s2, of)+plasticState(ph)%state(ns+s2, of))+&
abs(gdot_slip_pos(s2))*(plasticState(ph)%state(s1, of)+plasticState(ph)%state(ns+s1, of)))/&
(constitutive_dislokmc_L0(instance)*constitutive_dislokmc_burgersPerSlipSystem(j,instance))*&
(1.0_pReal-exp(-constitutive_dislokmc_VcrossSlip(instance)/(kB*Temperature)*&
(constitutive_dislokmc_tau_r(j,instance)-tau_twin(j))))
@ -1818,7 +1822,7 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
real(pReal) :: dvel_slip
real(pReal) :: StressRatio_u,StressRatio_uminus1
real(preal), dimension(constitutive_dislokmc_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip, vel_slip
gdot_slip_pos, vel_slip
real(pReal), dimension(3,3) :: eigVectors
real(pReal), dimension (3) :: eigValues
logical :: error
@ -1930,36 +1934,6 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
plasticState(ph)%state(5*ns+3*nt+j, of))
enddo; enddo
c = c + ns
case (resolved_stress_shearband_ID)
do j = 1_pInt,6_pInt ! loop over all shearband families
constitutive_dislokmc_postResults(c+j) = dot_product(Tstar_v, &
constitutive_dislokmc_sbSv(1:6,j,ipc,ip,el))
enddo
c = c + 6_pInt
case (shear_rate_shearband_ID)
do j = 1_pInt,6_pInt ! loop over all shearbands
!* Resolved shear stress on shearband system
tau = dot_product(Tstar_v,constitutive_dislokmc_sbSv(1:6,j,ipc,ip,el))
!* Stress ratios
if (abs(tau) < tol_math_check) then
StressRatio_p = 0.0_pReal
StressRatio_pminus1 = 0.0_pReal
else
StressRatio_p = (abs(tau)/constitutive_dislokmc_sbResistance(instance))**&
constitutive_dislokmc_pShearBand(instance)
StressRatio_pminus1 = (abs(tau)/constitutive_dislokmc_sbResistance(instance))**&
(constitutive_dislokmc_pShearBand(instance)-1.0_pReal)
endif
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislokmc_sbQedge(instance)/(kB*Temperature)
!* Initial shear rates
DotGamma0 = constitutive_dislokmc_sbVelocity(instance)
! Shear rate due to shear band
constitutive_dislokmc_postResults(c+j) = &
DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislokmc_qShearBand(instance))*&
sign(1.0_pReal,tau)
enddo
c = c + 6_pInt
case (twin_fraction_ID)
constitutive_dislokmc_postResults(c+1_pInt:c+nt) = plasticState(ph)%state((3_pInt*ns+1_pInt):(3_pInt*ns+nt), of)
c = c + nt
@ -2001,13 +1975,13 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
constitutive_dislokmc_v0PerSlipSystem(j,instance)
!* Shear rates due to slip
gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**&
gdot_slip_pos(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**&
constitutive_dislokmc_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau) &
* (1.0_pReal-constitutive_dislokmc_sPerSlipFamily(f,instance) &
* exp(-BoltzmannRatio*(1_pInt-StressRatio_p) ** constitutive_dislokmc_qPerSlipFamily(f,instance)) ) &
* StressRatio_u
else
gdot_slip(j) = 0.0_pReal
gdot_slip_pos(j) = 0.0_pReal
endif
enddo;enddo
@ -2030,8 +2004,8 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i)
s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i)
if (tau < constitutive_dislokmc_tau_r(j,instance)) then
Ndot0=(abs(gdot_slip(s1))*(plasticState(ph)%state(s2, of)+plasticState(ph)%state(ns+s2, of))+&
abs(gdot_slip(s2))*(plasticState(ph)%state(s1, of)+plasticState(ph)%state(ns+s1, of)))/&
Ndot0=(abs(gdot_slip_pos(s1))*(plasticState(ph)%state(s2, of)+plasticState(ph)%state(ns+s2, of))+&
abs(gdot_slip_pos(s2))*(plasticState(ph)%state(s1, of)+plasticState(ph)%state(ns+s1, of)))/&
(constitutive_dislokmc_L0(instance)*&
constitutive_dislokmc_burgersPerSlipSystem(j,instance))*&
(1.0_pReal-exp(-constitutive_dislokmc_VcrossSlip(instance)/(kB*Temperature)*&
@ -2111,7 +2085,7 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
* (1.0_pReal-constitutive_dislokmc_sPerSlipFamily(f,instance) &
* exp(-BoltzmannRatio*(1-StressRatio_p) ** constitutive_dislokmc_qPerSlipFamily(f,instance)))
gdot_slip(j) = (1.0_pReal - sumf) * DotGamma0 &
gdot_slip_pos(j) = (1.0_pReal - sumf) * DotGamma0 &
* StressRatio_u * vel_slip(j) &
* sign(1.0_pReal,tau)
!* Derivatives of shear rates
@ -2133,25 +2107,18 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
else
gdot_slip(j) = 0.0_pReal
gdot_slip_pos(j) = 0.0_pReal
dgdot_dtauslip = 0.0_pReal
endif
!* Stress exponent
if (gdot_slip(j)==0.0_pReal) then
if (gdot_slip_pos(j)==0.0_pReal) then
constitutive_dislokmc_postResults(c+j) = 0.0_pReal
else
constitutive_dislokmc_postResults(c+j) = (tau/gdot_slip(j))*dgdot_dtauslip
constitutive_dislokmc_postResults(c+j) = (tau/gdot_slip_pos(j))*dgdot_dtauslip
endif
enddo ; enddo
c = c + ns
case (sb_eigenvalues_ID)
forall (j = 1_pInt:3_pInt) &
constitutive_dislokmc_postResults(c+j) = eigValues(j)
c = c + 3_pInt
case (sb_eigenvectors_ID)
constitutive_dislokmc_postResults(c+1_pInt:c+9_pInt) = reshape(eigVectors,[9])
c = c + 9_pInt
end select
enddo
end function constitutive_dislokmc_postResults