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

fixed bug introduced with lattice_structure change nonlocal, but with DEBUG=ON OPTIMIZATION=OFF there is an FPE. Division by zero? Marked in the code 

forgot to commit dislotwin last time, now seems to work
This commit is contained in:
Martin Diehl 2014-03-11 23:55:40 +00:00
parent ef8fbf4dda
commit ff1b1c1a50
3 changed files with 312 additions and 182 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

422
code/constitutive_dislotwin.f90
View File

@ -77,10 +77,9 @@ module constitutive_dislotwin
constitutive_dislotwin_D0, & !< prefactor for self-diffusion coefficient constitutive_dislotwin_D0, & !< prefactor for self-diffusion coefficient
constitutive_dislotwin_Qsd, & !< activation energy for dislocation climb constitutive_dislotwin_Qsd, & !< activation energy for dislocation climb
constitutive_dislotwin_GrainSize, & !< grain size constitutive_dislotwin_GrainSize, & !< grain size
constitutive_dislotwin_p, & !< p-exponent in glide velocity constitutive_dislotwin_pShearBand, & !< p-exponent in shearband velocity
constitutive_dislotwin_q, & !< q-exponent in glide velocity constitutive_dislotwin_qShearBand, & !< q-exponent in shearband velocity
constitutive_dislotwin_MaxTwinFraction, & !< maximum allowed total twin volume fraction constitutive_dislotwin_MaxTwinFraction, & !< maximum allowed total twin volume fraction
constitutive_dislotwin_r, & !< r-exponent in twin nucleation rate
constitutive_dislotwin_CEdgeDipMinDistance, & !< constitutive_dislotwin_CEdgeDipMinDistance, & !<
constitutive_dislotwin_Cmfptwin, & !< constitutive_dislotwin_Cmfptwin, & !<
constitutive_dislotwin_Cthresholdtwin, & !< constitutive_dislotwin_Cthresholdtwin, & !<
@ -97,11 +96,9 @@ module constitutive_dislotwin
constitutive_dislotwin_aTolTwinFrac !< absolute tolerance for integration of twin volume fraction constitutive_dislotwin_aTolTwinFrac !< absolute tolerance for integration of twin volume fraction
real(pReal), dimension(:,:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:,:), allocatable, private :: &
constitutive_dislotwin_Ctwin_66 !< twin elasticity matrix in Mandel notation for each instance constitutive_dislotwin_Ctwin66 !< twin elasticity matrix in Mandel notation for each instance
real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: &
constitutive_dislotwin_Ctwin_3333 !< twin elasticity matrix for each instance constitutive_dislotwin_Ctwin3333 !< twin elasticity matrix for each instance
real(pReal), dimension(:,:), allocatable, private :: & real(pReal), dimension(:,:), allocatable, private :: &
constitutive_dislotwin_rhoEdge0, & !< initial edge dislocation density per slip system for each family and instance constitutive_dislotwin_rhoEdge0, & !< initial edge dislocation density per slip system for each family and instance
constitutive_dislotwin_rhoEdgeDip0, & !< initial edge dipole density per slip system for each family and instance constitutive_dislotwin_rhoEdgeDip0, & !< initial edge dipole density per slip system for each family and instance
@ -113,6 +110,7 @@ module constitutive_dislotwin
constitutive_dislotwin_QedgePerSlipSystem, & !< activation energy for glide [J] for each slip system and instance constitutive_dislotwin_QedgePerSlipSystem, & !< activation energy for glide [J] for each slip system and instance
constitutive_dislotwin_v0PerSlipFamily, & !< dislocation velocity prefactor [m/s] for each family and instance constitutive_dislotwin_v0PerSlipFamily, & !< dislocation velocity prefactor [m/s] for each family and instance
constitutive_dislotwin_v0PerSlipSystem, & !< dislocation velocity prefactor [m/s] for each slip system and instance constitutive_dislotwin_v0PerSlipSystem, & !< dislocation velocity prefactor [m/s] for each slip system and instance
constitutive_dislotwin_tau_peierlsPerSlipFamily, & !< Peierls stress [Pa] for each family and instance
constitutive_dislotwin_Ndot0PerTwinFamily, & !< twin nucleation rate [1/m³s] for each twin family and instance constitutive_dislotwin_Ndot0PerTwinFamily, & !< twin nucleation rate [1/m³s] for each twin family and instance
constitutive_dislotwin_Ndot0PerTwinSystem, & !< twin nucleation rate [1/m³s] for each twin system and instance constitutive_dislotwin_Ndot0PerTwinSystem, & !< twin nucleation rate [1/m³s] for each twin system and instance
constitutive_dislotwin_tau_r, & !< stress to bring partial close together for each twin system and instance constitutive_dislotwin_tau_r, & !< stress to bring partial close together for each twin system and instance
@ -123,7 +121,10 @@ module constitutive_dislotwin
constitutive_dislotwin_interaction_SlipSlip, & !< coefficients for slip-slip interaction for each interaction type and instance constitutive_dislotwin_interaction_SlipSlip, & !< coefficients for slip-slip interaction for each interaction type and instance
constitutive_dislotwin_interaction_SlipTwin, & !< coefficients for slip-twin interaction for each interaction type and instance constitutive_dislotwin_interaction_SlipTwin, & !< coefficients for slip-twin interaction for each interaction type and instance
constitutive_dislotwin_interaction_TwinSlip, & !< coefficients for twin-slip interaction for each interaction type and instance constitutive_dislotwin_interaction_TwinSlip, & !< coefficients for twin-slip interaction for each interaction type and instance
constitutive_dislotwin_interaction_TwinTwin !< coefficients for twin-twin interaction for each interaction type and instance constitutive_dislotwin_interaction_TwinTwin, & !< coefficients for twin-twin interaction for each interaction type and instance
constitutive_dislotwin_pPerSlipFamily, & !< p-exponent in glide velocity
constitutive_dislotwin_qPerSlipFamily, & !< q-exponent in glide velocity
constitutive_dislotwin_rPerTwinFamily !< r-exponent in twin nucleation rate
real(pReal), dimension(:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:), allocatable, private :: &
constitutive_dislotwin_interactionMatrix_SlipSlip, & !< interaction matrix of the different slip systems for each instance constitutive_dislotwin_interactionMatrix_SlipSlip, & !< interaction matrix of the different slip systems for each instance
constitutive_dislotwin_interactionMatrix_SlipTwin, & !< interaction matrix of slip systems with twin systems for each instance constitutive_dislotwin_interactionMatrix_SlipTwin, & !< interaction matrix of slip systems with twin systems for each instance
@ -225,6 +226,7 @@ subroutine constitutive_dislotwin_init(fileUnit)
character(len=65536) :: & character(len=65536) :: &
tag = '', & tag = '', &
line = '' line = ''
real(pReal), dimension(:), allocatable :: tempPerSlip, tempPerTwin
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>' write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
@ -253,10 +255,9 @@ subroutine constitutive_dislotwin_init(fileUnit)
allocate(constitutive_dislotwin_D0(maxNinstance), source=0.0_pReal) allocate(constitutive_dislotwin_D0(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_Qsd(maxNinstance), source=0.0_pReal) allocate(constitutive_dislotwin_Qsd(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_GrainSize(maxNinstance), source=0.0_pReal) allocate(constitutive_dislotwin_GrainSize(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_p(maxNinstance), source=0.0_pReal) allocate(constitutive_dislotwin_pShearBand(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_q(maxNinstance), source=0.0_pReal) allocate(constitutive_dislotwin_qShearBand(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_MaxTwinFraction(maxNinstance), source=0.0_pReal) allocate(constitutive_dislotwin_MaxTwinFraction(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_r(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_CEdgeDipMinDistance(maxNinstance), source=0.0_pReal) allocate(constitutive_dislotwin_CEdgeDipMinDistance(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_Cmfptwin(maxNinstance), source=0.0_pReal) allocate(constitutive_dislotwin_Cmfptwin(maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_Cthresholdtwin(maxNinstance), source=0.0_pReal) allocate(constitutive_dislotwin_Cthresholdtwin(maxNinstance), source=0.0_pReal)
@ -281,12 +282,17 @@ subroutine constitutive_dislotwin_init(fileUnit)
source=0.0_pReal) source=0.0_pReal)
allocate(constitutive_dislotwin_v0PerSlipFamily(lattice_maxNslipFamily,maxNinstance), & allocate(constitutive_dislotwin_v0PerSlipFamily(lattice_maxNslipFamily,maxNinstance), &
source=0.0_pReal) source=0.0_pReal)
allocate(constitutive_dislotwin_tau_peierlsPerSlipFamily(lattice_maxNslipFamily,maxNinstance), &
source=0.0_pReal)
allocate(constitutive_dislotwin_pPerSlipFamily(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal)
allocate(constitutive_dislotwin_qPerSlipFamily(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal)
allocate(constitutive_dislotwin_Ndot0PerTwinFamily(lattice_maxNtwinFamily,maxNinstance), & allocate(constitutive_dislotwin_Ndot0PerTwinFamily(lattice_maxNtwinFamily,maxNinstance), &
source=0.0_pReal) source=0.0_pReal)
allocate(constitutive_dislotwin_twinsizePerTwinFamily(lattice_maxNtwinFamily,maxNinstance), & allocate(constitutive_dislotwin_twinsizePerTwinFamily(lattice_maxNtwinFamily,maxNinstance), &
source=0.0_pReal) source=0.0_pReal)
allocate(constitutive_dislotwin_CLambdaSlipPerSlipFamily(lattice_maxNslipFamily,maxNinstance), & allocate(constitutive_dislotwin_CLambdaSlipPerSlipFamily(lattice_maxNslipFamily,maxNinstance), &
source=0.0_pReal) source=0.0_pReal)
allocate(constitutive_dislotwin_rPerTwinFamily(lattice_maxNtwinFamily,maxNinstance),source=0.0_pReal)
allocate(constitutive_dislotwin_interaction_SlipSlip(lattice_maxNinteraction,maxNinstance), & allocate(constitutive_dislotwin_interaction_SlipSlip(lattice_maxNinteraction,maxNinstance), &
source=0.0_pReal) source=0.0_pReal)
allocate(constitutive_dislotwin_interaction_SlipTwin(lattice_maxNinteraction,maxNinstance), & allocate(constitutive_dislotwin_interaction_SlipTwin(lattice_maxNinteraction,maxNinstance), &
@ -321,6 +327,10 @@ subroutine constitutive_dislotwin_init(fileUnit)
Nchunks_SlipTwin = maxval(lattice_interactionSlipTwin(:,:,phase)) Nchunks_SlipTwin = maxval(lattice_interactionSlipTwin(:,:,phase))
Nchunks_TwinSlip = maxval(lattice_interactionTwinSlip(:,:,phase)) Nchunks_TwinSlip = maxval(lattice_interactionTwinSlip(:,:,phase))
Nchunks_TwinTwin = maxval(lattice_interactionTwinTwin(:,:,phase)) Nchunks_TwinTwin = maxval(lattice_interactionTwinTwin(:,:,phase))
if(allocated(tempPerSlip)) deallocate(tempPerSlip)
if(allocated(tempPerTwin)) deallocate(tempPerTwin)
allocate(tempPerSlip(Nchunks_SlipFamilies))
allocate(tempPerTwin(Nchunks_TwinFamilies))
endif endif
cycle ! skip to next line cycle ! skip to next line
endif endif
@ -377,66 +387,78 @@ subroutine constitutive_dislotwin_init(fileUnit)
case default case default
call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//' ('//PLASTICITY_DISLOTWIN_label//')') call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//' ('//PLASTICITY_DISLOTWIN_label//')')
end select end select
!--------------------------------------------------------------------------------------------------
! parameters depending on slip number of slip system families
case ('nslip') case ('nslip')
if (positions(1) < 1_pInt + Nchunks_SlipFamilies) & if (positions(1) < Nchunks_SlipFamilies + 1_pInt) &
call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')')
if (positions(1) > Nchunks_SlipFamilies + 1_pInt) &
call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')')
Nchunks_SlipFamilies = positions(1) - 1_pInt Nchunks_SlipFamilies = positions(1) - 1_pInt
do j = 1_pInt, Nchunks_SlipFamilies do j = 1_pInt, Nchunks_SlipFamilies
constitutive_dislotwin_Nslip(j,instance) = IO_intValue(line,positions,1_pInt+j) constitutive_dislotwin_Nslip(j,instance) = IO_intValue(line,positions,1_pInt+j)
enddo enddo
case ('rhoedge0','rhoedgedip0','slipburgers','qedge','v0','clambdaslip','tau_peierls','p_slip','q_slip')
do j = 1_pInt, Nchunks_SlipFamilies
tempPerSlip(j) = IO_floatValue(line,positions,1_pInt+j)
enddo
select case(tag)
case ('rhoedge0')
constitutive_dislotwin_rhoEdge0(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('rhoedgedip0')
constitutive_dislotwin_rhoEdgeDip0(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('slipburgers')
constitutive_dislotwin_burgersPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('qedge')
constitutive_dislotwin_QedgePerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('v0')
constitutive_dislotwin_v0PerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('clambdaslip')
constitutive_dislotwin_CLambdaSlipPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('tau_peierls')
if (lattice_structure(phase) /= LATTICE_bcc_ID) &
call IO_warning(42_pInt,ext_msg=trim(tag)//' for non-bcc ('//PLASTICITY_DISLOTWIN_label//')')
constitutive_dislotwin_tau_peierlsPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('p_slip')
constitutive_dislotwin_pPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('q_slip')
constitutive_dislotwin_qPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies)
end select
!--------------------------------------------------------------------------------------------------
! parameters depending on slip number of twin families
case ('ntwin') case ('ntwin')
if (positions(1) < 1_pInt + Nchunks_TwinFamilies) & if (positions(1) < Nchunks_TwinFamilies + 1_pInt) &
call IO_warning(51_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')') call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')')
if (positions(1) > Nchunks_TwinFamilies + 1_pInt) &
call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOTWIN_label//')')
Nchunks_TwinFamilies = positions(1) - 1_pInt Nchunks_TwinFamilies = positions(1) - 1_pInt
do j = 1_pInt, Nchunks_TwinFamilies do j = 1_pInt, Nchunks_TwinFamilies
constitutive_dislotwin_Ntwin(j,instance) = IO_intValue(line,positions,1_pInt+j) constitutive_dislotwin_Ntwin(j,instance) = IO_intValue(line,positions,1_pInt+j)
enddo enddo
case ('rhoedge0') case ('ndot0','twinsize','twinburgers','r_twin')
do j = 1_pInt, Nchunks_SlipFamilies
constitutive_dislotwin_rhoEdge0(j,instance) = IO_floatValue(line,positions,1_pInt+j)
enddo
case ('rhoedgedip0')
do j = 1_pInt, Nchunks_SlipFamilies
constitutive_dislotwin_rhoEdgeDip0(j,instance) = IO_floatValue(line,positions,1_pInt+j)
enddo
case ('slipburgers')
do j = 1_pInt, Nchunks_SlipFamilies
constitutive_dislotwin_burgersPerSlipFamily(j,instance) = IO_floatValue(line,positions,1_pInt+j)
enddo
case ('twinburgers')
do j = 1_pInt, Nchunks_TwinFamilies
constitutive_dislotwin_burgersPerTwinFamily(j,instance) = IO_floatValue(line,positions,1_pInt+j)
enddo
case ('qedge')
do j = 1_pInt, Nchunks_SlipFamilies
constitutive_dislotwin_QedgePerSlipFamily(j,instance) = IO_floatValue(line,positions,1_pInt+j)
enddo
case ('v0')
do j = 1_pInt, Nchunks_SlipFamilies
constitutive_dislotwin_v0PerSlipFamily(j,instance) = IO_floatValue(line,positions,1_pInt+j)
enddo
case ('ndot0')
do j = 1_pInt, Nchunks_TwinFamilies do j = 1_pInt, Nchunks_TwinFamilies
constitutive_dislotwin_Ndot0PerTwinFamily(j,instance) = IO_floatValue(line,positions,1_pInt+j) tempPerTwin(j) = IO_floatValue(line,positions,1_pInt+j)
enddo
case ('twinsize')
do j = 1_pInt, Nchunks_TwinFamilies
constitutive_dislotwin_twinsizePerTwinFamily(j,instance) = IO_floatValue(line,positions,1_pInt+j)
enddo
case ('clambdaslip')
do j = 1_pInt, Nchunks_SlipFamilies
constitutive_dislotwin_CLambdaSlipPerSlipFamily(j,instance) = IO_floatValue(line,positions,1_pInt+j)
enddo enddo
select case(tag)
case ('ndot0')
if (lattice_structure(phase) == LATTICE_fcc_ID) &
call IO_warning(42_pInt,ext_msg=trim(tag)//' for fcc ('//PLASTICITY_DISLOTWIN_label//')')
constitutive_dislotwin_Ndot0PerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies)
case ('twinsize')
constitutive_dislotwin_twinsizePerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies)
case ('twinburgers')
constitutive_dislotwin_burgersPerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies)
case ('r_twin')
constitutive_dislotwin_rPerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies)
end select
case ('grainsize') case ('grainsize')
constitutive_dislotwin_GrainSize(instance) = IO_floatValue(line,positions,2_pInt) constitutive_dislotwin_GrainSize(instance) = IO_floatValue(line,positions,2_pInt)
case ('maxtwinfraction') case ('maxtwinfraction')
constitutive_dislotwin_MaxTwinFraction(instance) = IO_floatValue(line,positions,2_pInt) constitutive_dislotwin_MaxTwinFraction(instance) = IO_floatValue(line,positions,2_pInt)
case ('pexponent') case ('p_shearband')
constitutive_dislotwin_p(instance) = IO_floatValue(line,positions,2_pInt) constitutive_dislotwin_pShearBand(instance) = IO_floatValue(line,positions,2_pInt)
case ('qexponent') case ('q_shearband')
constitutive_dislotwin_q(instance) = IO_floatValue(line,positions,2_pInt) constitutive_dislotwin_qShearBand(instance) = IO_floatValue(line,positions,2_pInt)
case ('rexponent')
constitutive_dislotwin_r(instance) = IO_floatValue(line,positions,2_pInt)
case ('d0') case ('d0')
constitutive_dislotwin_D0(instance) = IO_floatValue(line,positions,2_pInt) constitutive_dislotwin_D0(instance) = IO_floatValue(line,positions,2_pInt)
case ('qsd') case ('qsd')
@ -518,6 +540,8 @@ subroutine constitutive_dislotwin_init(fileUnit)
call IO_error(211_pInt,el=instance,ext_msg='slipBurgers ('//PLASTICITY_DISLOTWIN_label//')') call IO_error(211_pInt,el=instance,ext_msg='slipBurgers ('//PLASTICITY_DISLOTWIN_label//')')
if (constitutive_dislotwin_v0PerSlipFamily(f,instance) <= 0.0_pReal) & if (constitutive_dislotwin_v0PerSlipFamily(f,instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='v0 ('//PLASTICITY_DISLOTWIN_label//')') call IO_error(211_pInt,el=instance,ext_msg='v0 ('//PLASTICITY_DISLOTWIN_label//')')
if (constitutive_dislotwin_tau_peierlsPerSlipFamily(f,instance) < 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='tau_peierls ('//PLASTICITY_DISLOTWIN_label//')')
endif endif
enddo enddo
do f = 1_pInt,lattice_maxNtwinFamily do f = 1_pInt,lattice_maxNtwinFamily
@ -534,16 +558,27 @@ subroutine constitutive_dislotwin_init(fileUnit)
call IO_error(211_pInt,el=instance,ext_msg='D0 ('//PLASTICITY_DISLOTWIN_label//')') call IO_error(211_pInt,el=instance,ext_msg='D0 ('//PLASTICITY_DISLOTWIN_label//')')
if (constitutive_dislotwin_Qsd(instance) <= 0.0_pReal) & if (constitutive_dislotwin_Qsd(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')') call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')')
if (constitutive_dislotwin_SFE_0K(instance) == 0.0_pReal .and. constitutive_dislotwin_dSFE_dT(instance) == 0.0_pReal) & if (sum(constitutive_dislotwin_Ntwin(:,instance)) > 0_pInt) then
call IO_error(211_pInt,el=instance,ext_msg='SFE ('//PLASTICITY_DISLOTWIN_label//')') if (constitutive_dislotwin_SFE_0K(instance) == 0.0_pReal .and. &
if (constitutive_dislotwin_aTolRho(instance) <= 0.0_pReal) & constitutive_dislotwin_dSFE_dT(instance) == 0.0_pReal .and. &
call IO_error(211_pInt,el=instance,ext_msg='aTolRho ('//PLASTICITY_DISLOTWIN_label//')') lattice_structure(phase) == LATTICE_fcc_ID) &
if (constitutive_dislotwin_aTolTwinFrac(instance) <= 0.0_pReal) & call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
call IO_error(211_pInt,el=instance,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOTWIN_label//')') if (constitutive_dislotwin_aTolRho(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='aTolRho ('//PLASTICITY_DISLOTWIN_label//')')
if (constitutive_dislotwin_aTolTwinFrac(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOTWIN_label//')')
endif
if (constitutive_dislotwin_sbResistance(instance) < 0.0_pReal) & if (constitutive_dislotwin_sbResistance(instance) < 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='sbResistance ('//PLASTICITY_DISLOTWIN_label//')') call IO_error(211_pInt,el=instance,ext_msg='sbResistance ('//PLASTICITY_DISLOTWIN_label//')')
if (constitutive_dislotwin_sbVelocity(instance) < 0.0_pReal) & if (constitutive_dislotwin_sbVelocity(instance) < 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='sbVelocity ('//PLASTICITY_DISLOTWIN_label//')') call IO_error(211_pInt,el=instance,ext_msg='sbVelocity ('//PLASTICITY_DISLOTWIN_label//')')
if (constitutive_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
constitutive_dislotwin_pShearBand(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='pShearBand ('//PLASTICITY_DISLOTWIN_label//')')
if (constitutive_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
constitutive_dislotwin_qShearBand(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='qShearBand ('//PLASTICITY_DISLOTWIN_label//')')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Determine total number of active slip or twin systems ! Determine total number of active slip or twin systems
constitutive_dislotwin_Nslip(:,instance) = min(lattice_NslipSystem(:,phase),constitutive_dislotwin_Nslip(:,instance)) constitutive_dislotwin_Nslip(:,instance) = min(lattice_NslipSystem(:,phase),constitutive_dislotwin_Nslip(:,instance))
@ -567,16 +602,23 @@ subroutine constitutive_dislotwin_init(fileUnit)
allocate(constitutive_dislotwin_twinsizePerTwinSystem(maxTotalNtwin, maxNinstance), source=0.0_pReal) allocate(constitutive_dislotwin_twinsizePerTwinSystem(maxTotalNtwin, maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_CLambdaSlipPerSlipSystem(maxTotalNslip, maxNinstance),source=0.0_pReal) allocate(constitutive_dislotwin_CLambdaSlipPerSlipSystem(maxTotalNslip, maxNinstance),source=0.0_pReal)
allocate(constitutive_dislotwin_interactionMatrix_SlipSlip(maxTotalNslip,maxTotalNslip,maxNinstance), & allocate(constitutive_dislotwin_interactionMatrix_SlipSlip(maxval(constitutive_dislotwin_totalNslip),& ! slip resistance from slip activity
source=0.0_pReal) maxval(constitutive_dislotwin_totalNslip),&
allocate(constitutive_dislotwin_interactionMatrix_SlipTwin(maxTotalNslip,maxTotalNtwin,maxNinstance), & maxNinstance), source=0.0_pReal)
source=0.0_pReal) allocate(constitutive_dislotwin_interactionMatrix_SlipTwin(maxval(constitutive_dislotwin_totalNslip),& ! slip resistance from twin activity
allocate(constitutive_dislotwin_interactionMatrix_TwinSlip(maxTotalNtwin,maxTotalNslip,maxNinstance), & maxval(constitutive_dislotwin_totalNtwin),&
source=0.0_pReal) maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_interactionMatrix_TwinTwin(maxTotalNtwin,maxTotalNtwin,maxNinstance), & allocate(constitutive_dislotwin_interactionMatrix_TwinSlip(maxval(constitutive_dislotwin_totalNtwin),& ! twin resistance from slip activity
source=0.0_pReal) maxval(constitutive_dislotwin_totalNslip),&
maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_interactionMatrix_TwinTwin(maxval(constitutive_dislotwin_totalNtwin),& ! twin resistance from twin activity
maxval(constitutive_dislotwin_totalNtwin),&
maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_forestProjectionEdge(maxTotalNslip,maxTotalNslip,maxNinstance), & allocate(constitutive_dislotwin_forestProjectionEdge(maxTotalNslip,maxTotalNslip,maxNinstance), &
source=0.0_pReal) source=0.0_pReal)
allocate(constitutive_dislotwin_Ctwin66(6,6,maxTotalNtwin,maxNinstance), source=0.0_pReal)
allocate(constitutive_dislotwin_Ctwin3333(3,3,3,3,maxTotalNtwin,maxNinstance), source=0.0_pReal)
initializeInstances: do phase = 1_pInt, size(phase_plasticity) initializeInstances: do phase = 1_pInt, size(phase_plasticity)
if (phase_plasticity(phase) == PLASTICITY_dislotwin_ID) then if (phase_plasticity(phase) == PLASTICITY_dislotwin_ID) then
instance = phase_plasticityInstance(phase) instance = phase_plasticityInstance(phase)
@ -705,19 +747,19 @@ subroutine constitutive_dislotwin_init(fileUnit)
!* Rotate twin elasticity matrices !* Rotate twin elasticity matrices
index_otherFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! index in full lattice twin list index_otherFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! index in full lattice twin list
do l = 1_pInt,3_pInt ; do m = 1_pInt,3_pInt ; do n = 1_pInt,3_pInt ; do o = 1_pInt,3_pInt do l = 1_pInt,3_pInt; do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt; do o = 1_pInt,3_pInt
do p = 1_pInt,3_pInt ; do q = 1_pInt,3_pInt ; do r = 1_pInt,3_pInt ; do s = 1_pInt,3_pInt do p = 1_pInt,3_pInt; do q = 1_pInt,3_pInt; do r = 1_pInt,3_pInt; do s = 1_pInt,3_pInt
constitutive_dislotwin_Ctwin_3333(l,m,n,o,index_myFamily+j,instance) = & constitutive_dislotwin_Ctwin3333(l,m,n,o,index_myFamily+j,instance) = &
constitutive_dislotwin_Ctwin_3333(l,m,n,o,index_myFamily+j,instance) + & constitutive_dislotwin_Ctwin3333(l,m,n,o,index_myFamily+j,instance) + &
lattice_C3333(p,q,r,s,instance) * & lattice_C3333(p,q,r,s,instance) * &
lattice_Qtwin(l,p,index_otherFamily+j,phase) * & lattice_Qtwin(l,p,index_otherFamily+j,phase) * &
lattice_Qtwin(m,q,index_otherFamily+j,phase) * & lattice_Qtwin(m,q,index_otherFamily+j,phase) * &
lattice_Qtwin(n,r,index_otherFamily+j,phase) * & lattice_Qtwin(n,r,index_otherFamily+j,phase) * &
lattice_Qtwin(o,s,index_otherFamily+j,phase) lattice_Qtwin(o,s,index_otherFamily+j,phase)
enddo ; enddo ; enddo ; enddo enddo; enddo; enddo; enddo
enddo ; enddo ; enddo ; enddo enddo; enddo; enddo; enddo
constitutive_dislotwin_Ctwin_66(1:6,1:6,index_myFamily+j,instance) = & constitutive_dislotwin_Ctwin66(1:6,1:6,index_myFamily+j,instance) = &
math_Mandel3333to66(constitutive_dislotwin_Ctwin_3333(1:3,1:3,1:3,1:3,index_myFamily+j,instance)) math_Mandel3333to66(constitutive_dislotwin_Ctwin3333(1:3,1:3,1:3,1:3,index_myFamily+j,instance))
!* Interaction matrices !* Interaction matrices
do o = 1_pInt,lattice_maxNslipFamily do o = 1_pInt,lattice_maxNslipFamily
@ -757,7 +799,9 @@ function constitutive_dislotwin_stateInit(instance,phase)
pi pi
use lattice, only: & use lattice, only: &
lattice_maxNslipFamily, & lattice_maxNslipFamily, &
lattice_mu lattice_structure, &
lattice_mu, &
lattice_bcc_ID
implicit none implicit none
integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity
@ -790,7 +834,7 @@ function constitutive_dislotwin_stateInit(instance,phase)
index_myFamily+constitutive_dislotwin_Nslip(f,instance)) = & index_myFamily+constitutive_dislotwin_Nslip(f,instance)) = &
constitutive_dislotwin_rhoEdgeDip0(f,instance) constitutive_dislotwin_rhoEdgeDip0(f,instance)
enddo enddo
constitutive_dislotwin_stateInit(1_pInt:ns) = rhoEdge0 constitutive_dislotwin_stateInit(1_pInt:ns) = rhoEdge0
constitutive_dislotwin_stateInit(ns+1_pInt:2_pInt*ns) = rhoEdgeDip0 constitutive_dislotwin_stateInit(ns+1_pInt:2_pInt*ns) = rhoEdgeDip0
@ -806,11 +850,23 @@ function constitutive_dislotwin_stateInit(instance,phase)
constitutive_dislotwin_GrainSize(instance)/(1.0_pReal+invLambdaSlip0(i)*constitutive_dislotwin_GrainSize(instance)) constitutive_dislotwin_GrainSize(instance)/(1.0_pReal+invLambdaSlip0(i)*constitutive_dislotwin_GrainSize(instance))
constitutive_dislotwin_stateInit(5_pInt*ns+3_pInt*nt+1:6_pInt*ns+3_pInt*nt) = MeanFreePathSlip0 constitutive_dislotwin_stateInit(5_pInt*ns+3_pInt*nt+1:6_pInt*ns+3_pInt*nt) = MeanFreePathSlip0
forall (i = 1_pInt:ns) & if (lattice_structure(phase) == lattice_bcc_ID) then
tauSlipThreshold0(i) = constitutive_dislotwin_SolidSolutionStrength(instance) + & j = 0_pInt
lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(i,instance) * & slipFamiliesLoop: do f = 1_pInt,lattice_maxNslipFamily
sqrt(dot_product((rhoEdge0+rhoEdgeDip0),constitutive_dislotwin_interactionMatrix_SlipSlip(i,1:ns,instance))) slipSystemsLoop: do i = 1_pInt,constitutive_dislotwin_Nslip(f,instance)
constitutive_dislotwin_stateInit(6_pInt*ns+4_pInt*nt+1:7_pInt*ns+4_pInt*nt) = tauSlipThreshold0 j = j+1_pInt
tauSlipThreshold0(i) = constitutive_dislotwin_tau_peierlsPerSlipFamily(f,instance)
enddo slipSystemsLoop
enddo slipFamiliesLoop
else
forall (i = 1_pInt:ns) &
tauSlipThreshold0(i) = constitutive_dislotwin_SolidSolutionStrength(instance) + &
lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(i,instance) * &
sqrt(dot_product((rhoEdge0+rhoEdgeDip0),constitutive_dislotwin_interactionMatrix_SlipSlip(i,1:ns,instance)))
endif
constitutive_dislotwin_stateInit(6_pInt*ns+4_pInt*nt+1:7_pInt*ns+4_pInt*nt) = tauSlipThreshold0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize dependent twin microstructural variables ! initialize dependent twin microstructural variables
@ -822,7 +878,7 @@ function constitutive_dislotwin_stateInit(instance,phase)
TwinVolume0(j) = & TwinVolume0(j) = &
(pi/4.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(j,instance)*MeanFreePathTwin0(j)**(2.0_pReal) (pi/4.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(j,instance)*MeanFreePathTwin0(j)**(2.0_pReal)
constitutive_dislotwin_stateInit(7_pInt*ns+5_pInt*nt+1_pInt:7_pInt*ns+6_pInt*nt) = TwinVolume0 constitutive_dislotwin_stateInit(7_pInt*ns+5_pInt*nt+1_pInt:7_pInt*ns+6_pInt*nt) = TwinVolume0
end function constitutive_dislotwin_stateInit end function constitutive_dislotwin_stateInit
@ -923,8 +979,11 @@ subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el)
material_phase, & material_phase, &
phase_plasticityInstance phase_plasticityInstance
use lattice, only: & use lattice, only: &
lattice_structure, &
lattice_mu, & lattice_mu, &
lattice_nu lattice_nu, &
lattice_bcc_ID
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -1008,7 +1067,7 @@ subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el)
(1.0_pReal+constitutive_dislotwin_GrainSize(instance)*(state(ipc,ip,el)%p(3_pInt*ns+s))) (1.0_pReal+constitutive_dislotwin_GrainSize(instance)*(state(ipc,ip,el)%p(3_pInt*ns+s)))
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
forall (t = 1_pInt:nt) & forall (t = 1_pInt:nt) &
state(ipc,ip,el)%p(6_pInt*ns+3_pInt*nt+t) = & state(ipc,ip,el)%p(6_pInt*ns+3_pInt*nt+t) = &
@ -1016,11 +1075,13 @@ subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el)
(1.0_pReal+constitutive_dislotwin_GrainSize(instance)*state(ipc,ip,el)%p(5_pInt*ns+2_pInt*nt+t)) (1.0_pReal+constitutive_dislotwin_GrainSize(instance)*state(ipc,ip,el)%p(5_pInt*ns+2_pInt*nt+t))
!* threshold stress for dislocation motion !* threshold stress for dislocation motion
forall (s = 1_pInt:ns) & if(lattice_structure(phase) /= LATTICE_BCC_ID) then ! bcc value remains constant
state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+s) = constitutive_dislotwin_SolidSolutionStrength(instance)+ & forall (s = 1_pInt:ns) &
lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(s,instance)*& state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+s) = constitutive_dislotwin_SolidSolutionStrength(instance)+ &
sqrt(dot_product((state(ipc,ip,el)%p(1:ns)+state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)),& lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(s,instance)*&
constitutive_dislotwin_interactionMatrix_SlipSlip(s,1:ns,instance))) sqrt(dot_product((state(ipc,ip,el)%p(1:ns)+state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)),&
constitutive_dislotwin_interactionMatrix_SlipSlip(s,1:ns,instance)))
endif
!* threshold stress for growing twin !* threshold stress for growing twin
forall (t = 1_pInt:nt) & forall (t = 1_pInt:nt) &
@ -1052,7 +1113,8 @@ end subroutine constitutive_dislotwin_microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el) subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec, &
tol_math_check
use math, only: & use math, only: &
math_Plain3333to99, & math_Plain3333to99, &
math_Mandel6to33, & math_Mandel6to33, &
@ -1146,10 +1208,17 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
!* Calculation of Lp !* Calculation of Lp
!* Resolved shear stress on slip system !* Resolved shear stress on slip system
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase)) tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase))
!* Stress ratios !* Stress ratios
StressRatio_p = (abs(tau_slip(j))/state(ipc,ip,el)%p(6*ns+4*nt+j))**constitutive_dislotwin_p(instance) if (abs(tau_slip(j)) < tol_math_check) then
StressRatio_pminus1 = (abs(tau_slip(j))/state(ipc,ip,el)%p(6*ns+4*nt+j))**(constitutive_dislotwin_p(instance)-1.0_pReal) StressRatio_p = 0.0_pReal
StressRatio_pminus1 = 0.0_pReal
else
StressRatio_p = (abs(tau_slip(j))/state(ipc,ip,el)%p(6*ns+4*nt+j))&
**constitutive_dislotwin_pPerSlipFamily(f,instance)
StressRatio_pminus1 = (abs(tau_slip(j))/state(ipc,ip,el)%p(6*ns+4*nt+j))&
**(constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal)
endif
!* Boltzmann ratio !* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature)
!* Initial shear rates !* Initial shear rates
@ -1159,14 +1228,14 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
!* Shear rates due to slip !* Shear rates due to slip
gdot_slip(j) = (1.0_pReal - sumf) * DotGamma0 & gdot_slip(j) = (1.0_pReal - sumf) * DotGamma0 &
* exp(-BoltzmannRatio*(1-StressRatio_p) ** constitutive_dislotwin_q(instance)) & * exp(-BoltzmannRatio*(1-StressRatio_p) ** constitutive_dislotwin_qPerSlipFamily(f,instance)) &
* sign(1.0_pReal,tau_slip(j)) * sign(1.0_pReal,tau_slip(j))
!* Derivatives of shear rates !* Derivatives of shear rates
dgdot_dtauslip(j) = & dgdot_dtauslip(j) = &
((abs(gdot_slip(j))*BoltzmannRatio*& ((abs(gdot_slip(j))*BoltzmannRatio*constitutive_dislotwin_pPerSlipFamily(f,instance)&
constitutive_dislotwin_p(instance)*constitutive_dislotwin_q(instance))/state(ipc,ip,el)%p(6*ns+4*nt+j))*& *constitutive_dislotwin_qPerSlipFamily(f,instance))/state(ipc,ip,el)%p(6*ns+4*nt+j))*&
StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_q(instance)-1.0_pReal) StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_qPerSlipFamily(f,instance)-1.0_pReal)
!* Plastic velocity gradient for dislocation glide !* Plastic velocity gradient for dislocation glide
Lp = Lp + gdot_slip(j)*lattice_Sslip(:,:,1,index_myFamily+i,phase) Lp = Lp + gdot_slip(j)*lattice_Sslip(:,:,1,index_myFamily+i,phase)
@ -1197,23 +1266,31 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
tau_sb(j) = dot_product(Tstar_v,constitutive_dislotwin_sbSv(1:6,j,ipc,ip,el)) tau_sb(j) = dot_product(Tstar_v,constitutive_dislotwin_sbSv(1:6,j,ipc,ip,el))
!* Stress ratios !* Stress ratios
StressRatio_p = (abs(tau_sb(j))/constitutive_dislotwin_sbResistance(instance))**constitutive_dislotwin_p(instance) if (abs(tau_sb(j)) < tol_math_check) then
StressRatio_pminus1 = (abs(tau_sb(j))/constitutive_dislotwin_sbResistance(instance))& StressRatio_p = 0.0_pReal
**(constitutive_dislotwin_p(instance)-1.0_pReal) StressRatio_pminus1 = 0.0_pReal
else
StressRatio_p = (abs(tau_sb(j))/constitutive_dislotwin_sbResistance(instance))&
**constitutive_dislotwin_pShearBand(instance)
StressRatio_pminus1 = (abs(tau_sb(j))/constitutive_dislotwin_sbResistance(instance))&
**(constitutive_dislotwin_pShearBand(instance)-1.0_pReal)
endif
!* Boltzmann ratio !* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_sbQedge(instance)/(kB*Temperature) BoltzmannRatio = constitutive_dislotwin_sbQedge(instance)/(kB*Temperature)
!* Initial shear rates !* Initial shear rates
DotGamma0 = constitutive_dislotwin_sbVelocity(instance) DotGamma0 = constitutive_dislotwin_sbVelocity(instance)
!* Shear rates due to shearband !* Shear rates due to shearband
gdot_sb(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislotwin_q(instance))*& gdot_sb(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**&
sign(1.0_pReal,tau_sb(j)) constitutive_dislotwin_qShearBand(instance))*sign(1.0_pReal,tau_sb(j))
!* Derivatives of shear rates !* Derivatives of shear rates
dgdot_dtausb(j) = & dgdot_dtausb(j) = &
((abs(gdot_sb(j))*BoltzmannRatio*& ((abs(gdot_sb(j))*BoltzmannRatio*&
constitutive_dislotwin_p(instance)*constitutive_dislotwin_q(instance))/constitutive_dislotwin_sbResistance(instance))*& constitutive_dislotwin_pShearBand(instance)*constitutive_dislotwin_qShearBand(instance))/&
StressRatio_pminus1*(1_pInt-StressRatio_p)**(constitutive_dislotwin_q(instance)-1.0_pReal) constitutive_dislotwin_sbResistance(instance))*&
StressRatio_pminus1*(1_pInt-StressRatio_p)**(constitutive_dislotwin_qShearBand(instance)-1.0_pReal)
!* Plastic velocity gradient for shear banding !* Plastic velocity gradient for shear banding
Lp = Lp + gdot_sb(j)*sb_Smatrix Lp = Lp + gdot_sb(j)*sb_Smatrix
@ -1238,13 +1315,13 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
!* Calculation of Lp !* Calculation of Lp
!* Resolved shear stress on twin system !* Resolved shear stress on twin system
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,phase)) tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,phase))
!* Stress ratios !* Stress ratios
StressRatio_r = (state(ipc,ip,el)%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_r(instance) if (tau_twin(j) > tol_math_check) then
StressRatio_r = (state(ipc,ip,el)%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_rPerTwinFamily(f,instance)
!* Shear rates and their derivatives due to twin !* Shear rates and their derivatives due to twin
if ( tau_twin(j) > 0.0_pReal ) then
select case(lattice_structure(phase)) select case(lattice_structure(phase))
case (LATTICE_fcc_ID) case (LATTICE_fcc_ID)
s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i)
@ -1264,7 +1341,7 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
gdot_twin(j) = & gdot_twin(j) = &
(constitutive_dislotwin_MaxTwinFraction(instance)-sumf)*lattice_shearTwin(index_myFamily+i,phase)*& (constitutive_dislotwin_MaxTwinFraction(instance)-sumf)*lattice_shearTwin(index_myFamily+i,phase)*&
state(ipc,ip,el)%p(7*ns+5*nt+j)*Ndot0*exp(-StressRatio_r) state(ipc,ip,el)%p(7*ns+5*nt+j)*Ndot0*exp(-StressRatio_r)
dgdot_dtautwin(j) = ((gdot_twin(j)*constitutive_dislotwin_r(instance))/tau_twin(j))*StressRatio_r dgdot_dtautwin(j) = ((gdot_twin(j)*constitutive_dislotwin_rPerTwinFamily(f,instance))/tau_twin(j))*StressRatio_r
endif endif
!* Plastic velocity gradient for mechanical twinning !* Plastic velocity gradient for mechanical twinning
@ -1289,7 +1366,8 @@ end subroutine constitutive_dislotwin_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,el) pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec, &
tol_math_check
use math, only: & use math, only: &
pi pi
use mesh, only: & use mesh, only: &
@ -1310,7 +1388,8 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
lattice_mu, & lattice_mu, &
lattice_structure, & lattice_structure, &
lattice_fcc_twinNucleationSlipPair, & lattice_fcc_twinNucleationSlipPair, &
LATTICE_fcc_ID LATTICE_fcc_ID, &
LATTICE_bcc_ID
implicit none implicit none
real(pReal), dimension(6), intent(in):: & real(pReal), dimension(6), intent(in):: &
@ -1357,11 +1436,16 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
!* Resolved shear stress on slip system !* Resolved shear stress on slip system
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase)) tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase))
!* Stress ratios
StressRatio_p = (abs(tau_slip(j))/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& if (abs(tau_slip(j)) < tol_math_check) then
constitutive_dislotwin_p(instance) StressRatio_p = 0.0_pReal
StressRatio_pminus1 = (abs(tau_slip(j))/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& StressRatio_pminus1 = 0.0_pReal
(constitutive_dislotwin_p(instance)-1.0_pReal) else
StressRatio_p = (abs(tau_slip(j))/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
constitutive_dislotwin_pPerSlipFamily(f,instance)
StressRatio_pminus1 = (abs(tau_slip(j))/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
(constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal)
endif
!* Boltzmann ratio !* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature)
!* Initial shear rates !* Initial shear rates
@ -1370,8 +1454,8 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
constitutive_dislotwin_v0PerSlipSystem(j,instance) constitutive_dislotwin_v0PerSlipSystem(j,instance)
!* Shear rates due to slip !* Shear rates due to slip
gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislotwin_q(instance))*& gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)** &
sign(1.0_pReal,tau_slip(j)) constitutive_dislotwin_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau_slip(j))
!* Multiplication !* Multiplication
DotRhoMultiplication(j) = abs(gdot_slip(j))/& DotRhoMultiplication(j) = abs(gdot_slip(j))/&
@ -1442,10 +1526,10 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
!* Resolved shear stress on twin system !* Resolved shear stress on twin system
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,phase)) tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,phase))
!* Stress ratios !* Stress ratios
StressRatio_r = (state(ipc,ip,el)%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_r(instance) if (tau_twin(j) > tol_math_check) then
StressRatio_r = (state(ipc,ip,el)%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_rPerTwinFamily(f,instance)
!* Shear rates and their derivatives due to twin !* Shear rates and their derivatives due to twin
if ( tau_twin(j) > 0.0_pReal ) then
select case(lattice_structure(phase)) select case(lattice_structure(phase))
case (LATTICE_fcc_ID) case (LATTICE_fcc_ID)
s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i)
@ -1466,7 +1550,7 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e
(constitutive_dislotwin_MaxTwinFraction(instance)-sumf)*& (constitutive_dislotwin_MaxTwinFraction(instance)-sumf)*&
state(ipc,ip,el)%p(7_pInt*ns+5_pInt*nt+j)*Ndot0*exp(-StressRatio_r) state(ipc,ip,el)%p(7_pInt*ns+5_pInt*nt+j)*Ndot0*exp(-StressRatio_r)
!* Dotstate for accumulated shear due to twin !* Dotstate for accumulated shear due to twin
constitutive_dislotwin_dotstate(3_pInt*ns+nt+j) = constitutive_dislotwin_dotState(3_pInt*ns+j) * & constitutive_dislotwin_dotState(3_pInt*ns+nt+j) = constitutive_dislotwin_dotState(3_pInt*ns+j) * &
lattice_sheartwin(index_myfamily+i,phase) lattice_sheartwin(index_myfamily+i,phase)
endif endif
enddo enddo
@ -1480,7 +1564,8 @@ end function constitutive_dislotwin_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec, &
tol_math_check
use math, only: & use math, only: &
pi, & pi, &
math_Mandel6to33, & math_Mandel6to33, &
@ -1568,10 +1653,15 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
!* Resolved shear stress on slip system !* Resolved shear stress on slip system
tau = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase)) tau = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase))
!* Stress ratios !* Stress ratios
StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& if (abs(tau) < tol_math_check) then
constitutive_dislotwin_p(instance) StressRatio_p = 0.0_pReal
StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& StressRatio_pminus1 = 0.0_pReal
(constitutive_dislotwin_p(instance)-1.0_pReal) else
StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
constitutive_dislotwin_pPerSlipFamily(f,instance)
StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
(constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal)
endif
!* Boltzmann ratio !* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature)
!* Initial shear rates !* Initial shear rates
@ -1582,7 +1672,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
!* Shear rates due to slip !* Shear rates due to slip
constitutive_dislotwin_postResults(c+j) = & constitutive_dislotwin_postResults(c+j) = &
DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**& DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**&
constitutive_dislotwin_q(instance))*sign(1.0_pReal,tau) constitutive_dislotwin_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau)
enddo ; enddo enddo ; enddo
c = c + ns c = c + ns
case (accumulated_shear_slip_ID) case (accumulated_shear_slip_ID)
@ -1616,8 +1706,8 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
constitutive_dislotwin_postResults(c+j) = & constitutive_dislotwin_postResults(c+j) = &
(3.0_pReal*lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(j,instance))/& (3.0_pReal*lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(j,instance))/&
(16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase)))) (16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase))))
constitutive_dislotwin_postResults(c+j) = min(constitutive_dislotwin_postResults(c+j),state(ipc,ip,el)%p(5*ns+3*nt+j)) constitutive_dislotwin_postResults(c+j)=min(constitutive_dislotwin_postResults(c+j),state(ipc,ip,el)%p(5*ns+3*nt+j))
! constitutive_dislotwin_postResults(c+j) = max(constitutive_dislotwin_postResults(c+j),state(ipc,ip,el)%p(4*ns+2*nt+j)) ! constitutive_dislotwin_postResults(c+j)=max(constitutive_dislotwin_postResults(c+j),state(ipc,ip,el)%p(4*ns+2*nt+j))
enddo; enddo enddo; enddo
c = c + ns c = c + ns
case (resolved_stress_shearband_ID) case (resolved_stress_shearband_ID)
@ -1630,17 +1720,23 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
!* Resolved shear stress on shearband system !* Resolved shear stress on shearband system
tau = dot_product(Tstar_v,constitutive_dislotwin_sbSv(1:6,j,ipc,ip,el)) tau = dot_product(Tstar_v,constitutive_dislotwin_sbSv(1:6,j,ipc,ip,el))
!* Stress ratios !* Stress ratios
StressRatio_p = (abs(tau)/constitutive_dislotwin_sbResistance(instance))**constitutive_dislotwin_p(instance) if (abs(tau) < tol_math_check) then
StressRatio_pminus1 = (abs(tau)/constitutive_dislotwin_sbResistance(instance))& StressRatio_p = 0.0_pReal
**(constitutive_dislotwin_p(instance)-1.0_pReal) StressRatio_pminus1 = 0.0_pReal
else
StressRatio_p = (abs(tau)/constitutive_dislotwin_sbResistance(instance))**&
constitutive_dislotwin_pShearBand(instance)
StressRatio_pminus1 = (abs(tau)/constitutive_dislotwin_sbResistance(instance))**&
(constitutive_dislotwin_pShearBand(instance)-1.0_pReal)
endif
!* Boltzmann ratio !* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_sbQedge(instance)/(kB*Temperature) BoltzmannRatio = constitutive_dislotwin_sbQedge(instance)/(kB*Temperature)
!* Initial shear rates !* Initial shear rates
DotGamma0 = constitutive_dislotwin_sbVelocity(instance) DotGamma0 = constitutive_dislotwin_sbVelocity(instance)
! Shear rate due to shear band
!* Shear rates due to slip
constitutive_dislotwin_postResults(c+j) = & constitutive_dislotwin_postResults(c+j) = &
DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislotwin_q(instance))*sign(1.0_pReal,tau) DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislotwin_qShearBand(instance))*&
sign(1.0_pReal,tau)
enddo enddo
c = c + 6_pInt c = c + 6_pInt
case (twin_fraction_ID) case (twin_fraction_ID)
@ -1658,10 +1754,15 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
!* Resolved shear stress on slip system !* Resolved shear stress on slip system
tau = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase)) tau = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase))
!* Stress ratios !* Stress ratios
StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**& if (abs(tau) < tol_math_check) then
constitutive_dislotwin_p(instance) StressRatio_p = 0.0_pReal
StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**& StressRatio_pminus1 = 0.0_pReal
(constitutive_dislotwin_p(instance)-1.0_pReal) else
StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
constitutive_dislotwin_pPerSlipFamily(f,instance)
StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
(constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal)
endif
!* Boltzmann ratio !* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature)
!* Initial shear rates !* Initial shear rates
@ -1671,7 +1772,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
!* Shear rates due to slip !* Shear rates due to slip
gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**& gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**&
constitutive_dislotwin_q(instance))*sign(1.0_pReal,tau) constitutive_dislotwin_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau)
enddo;enddo enddo;enddo
j = 0_pInt j = 0_pInt
@ -1683,7 +1784,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
!* Resolved shear stress on twin system !* Resolved shear stress on twin system
tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,phase)) tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,phase))
!* Stress ratios !* Stress ratios
StressRatio_r = (state(ipc,ip,el)%p(7_pInt*ns+4_pInt*nt+j)/tau)**constitutive_dislotwin_r(instance) StressRatio_r = (state(ipc,ip,el)%p(7_pInt*ns+4_pInt*nt+j)/tau)**constitutive_dislotwin_rPerTwinFamily(f,instance)
!* Shear rates due to twin !* Shear rates due to twin
if ( tau > 0.0_pReal ) then if ( tau > 0.0_pReal ) then
@ -1738,31 +1839,36 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,constitutive_dislotwin_Nslip(f,instance) ! process each (active) slip system in family do i = 1_pInt,constitutive_dislotwin_Nslip(f,instance) ! process each (active) slip system in family
j = j + 1_pInt j = j + 1_pInt
!* Resolved shear stress on slip system !* Resolved shear stress on slip system
tau = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase)) tau = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase))
!* Stress ratios if (abs(tau) < tol_math_check) then
StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& StressRatio_p = 0.0_pReal
constitutive_dislotwin_p(instance) StressRatio_pminus1 = 0.0_pReal
StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& else
(constitutive_dislotwin_p(instance)-1.0_pReal) StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
constitutive_dislotwin_pPerSlipFamily(f,instance)
StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**&
(constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal)
endif
!* Boltzmann ratio !* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature)
!* Initial shear rates !* Initial shear rates
DotGamma0 = & DotGamma0 = &
state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* & state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* &
constitutive_dislotwin_v0PerSlipSystem(j,instance) constitutive_dislotwin_v0PerSlipSystem(j,instance)
!* Shear rates due to slip !* Shear rates due to slip
gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**& gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**&
constitutive_dislotwin_q(instance))*sign(1.0_pReal,tau) constitutive_dislotwin_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau)
!* Derivatives of shear rates !* Derivatives of shear rates
dgdot_dtauslip = & dgdot_dtauslip = &
((abs(gdot_slip(j))*BoltzmannRatio*& ((abs(gdot_slip(j))*BoltzmannRatio*&
constitutive_dislotwin_p(instance)*constitutive_dislotwin_q(instance))/state(ipc,ip,el)%p(6*ns+4*nt+j))*& constitutive_dislotwin_pPerSlipFamily(f,instance)*constitutive_dislotwin_qPerSlipFamily(f,instance))/&
StressRatio_pminus1*(1_pInt-StressRatio_p)**(constitutive_dislotwin_q(instance)-1.0_pReal) state(ipc,ip,el)%p(6*ns+4*nt+j))*StressRatio_pminus1*(1_pInt-StressRatio_p)**&
(constitutive_dislotwin_qPerSlipFamily(f,instance)-1.0_pReal)
!* Stress exponent !* Stress exponent
if (gdot_slip(j)==0.0_pReal) then if (gdot_slip(j)==0.0_pReal) then
constitutive_dislotwin_postResults(c+j) = 0.0_pReal constitutive_dislotwin_postResults(c+j) = 0.0_pReal

47
code/constitutive_nonlocal.f90
View File

@ -91,9 +91,6 @@ iRhoD, & !< state in
iV, & !< state indices for dislcation velocities iV, & !< state indices for dislcation velocities
iD !< state indices for stable dipole height iD !< state indices for stable dipole height
integer(pInt), dimension(:), allocatable, public :: &
constitutive_nonlocal_structure !< number representing the kind of lattice structure
integer(pInt), dimension(:), allocatable, private :: & integer(pInt), dimension(:), allocatable, private :: &
totalNslip !< total number of active slip systems for each instance totalNslip !< total number of active slip systems for each instance
@ -752,6 +749,7 @@ allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstances), s
sanityChecks: do phase = 1_pInt, size(phase_plasticity) sanityChecks: do phase = 1_pInt, size(phase_plasticity)
myPhase: if (phase_plasticity(phase) == PLASTICITY_NONLOCAL_ID) then myPhase: if (phase_plasticity(phase) == PLASTICITY_NONLOCAL_ID) then
instance = phase_plasticityInstance(phase)
if (sum(Nslip(:,instance)) <= 0_pInt) & if (sum(Nslip(:,instance)) <= 0_pInt) &
call IO_error(211_pInt,ext_msg='Nslip ('//PLASTICITY_NONLOCAL_label//')') call IO_error(211_pInt,ext_msg='Nslip ('//PLASTICITY_NONLOCAL_label//')')
do o = 1_pInt,maxval(phase_Noutput) do o = 1_pInt,maxval(phase_Noutput)
@ -842,7 +840,7 @@ allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstances), s
!*** determine total number of active slip systems !*** determine total number of active slip systems
Nslip(1:lattice_maxNslipFamily,instance) = min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase), & Nslip(1:lattice_maxNslipFamily,instance) = min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase), &
Nslip(1:lattice_maxNslipFamily,instance) ) ! we can't use more slip systems per family than specified in lattice Nslip(1:lattice_maxNslipFamily,instance) ) ! we can't use more slip systems per family than specified in lattice
totalNslip(instance) = sum(Nslip(1:lattice_maxNslipFamily,instance)) totalNslip(instance) = sum(Nslip(1:lattice_maxNslipFamily,instance))
endif myPhase endif myPhase
enddo sanityChecks enddo sanityChecks
@ -1194,8 +1192,6 @@ maxNinstances = int(count(phase_plasticity == PLASTICITY_NONLOCAL_ID),pInt)
do e = 1_pInt,mesh_NcpElems do e = 1_pInt,mesh_NcpElems
do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e)))
if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e))) & if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e))) &
write(6,*) shape(state(1,i,e)%p)
flush(6)
state(1,i,e)%p = 0.0_pReal state(1,i,e)%p = 0.0_pReal
enddo enddo
enddo enddo
@ -2111,6 +2107,45 @@ dLower = minDipoleHeight(1:ns,1:2,instance)
dUpper(1:ns,1) = lattice_mu(phase) * burgers(1:ns,instance) & dUpper(1:ns,1) = lattice_mu(phase) * burgers(1:ns,instance) &
/ (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau)) / (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau))
dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) / (4.0_pReal * pi * abs(tau)) dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) / (4.0_pReal * pi * abs(tau))
!in the test there is an FPE exception. Divistion by zero?
forall (c = 1_pInt:2_pInt) & forall (c = 1_pInt:2_pInt) &
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &

25
code/crystallite.f90
View File

@ -3293,8 +3293,7 @@ subroutine crystallite_orientations
use material, only: & use material, only: &
material_phase, & material_phase, &
homogenization_Ngrains, & homogenization_Ngrains, &
phase_localPlasticity, & phase_localPlasticity
phase_plasticityInstance
use mesh, only: & use mesh, only: &
mesh_element, & mesh_element, &
mesh_ipNeighborhood, & mesh_ipNeighborhood, &
@ -3305,7 +3304,6 @@ subroutine crystallite_orientations
lattice_qDisorientation, & lattice_qDisorientation, &
lattice_structure lattice_structure
use constitutive_nonlocal, only: & use constitutive_nonlocal, only: &
constitutive_nonlocal_structure, &
constitutive_nonlocal_updateCompatibility constitutive_nonlocal_updateCompatibility
@ -3318,11 +3316,7 @@ subroutine crystallite_orientations
neighboring_e, & ! element index of my neighbor neighboring_e, & ! element index of my neighbor
neighboring_i, & ! integration point index of my neighbor neighboring_i, & ! integration point index of my neighbor
myPhase, & ! phase myPhase, & ! phase
neighboringPhase, & neighboringPhase
myInstance, & ! instance of plasticity
neighboringInstance, &
myStructure, & ! lattice structure
neighboringStructure
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
U, & U, &
R R
@ -3357,16 +3351,13 @@ subroutine crystallite_orientations
! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL --- ! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL ---
! --- we use crystallite_orientation from above, so need a seperate loop ! --- we use crystallite_orientation from above, so need a separate loop
!$OMP PARALLEL DO PRIVATE(myPhase,myInstance,myStructure,neighboring_e,neighboring_i,neighboringPhase,& !$OMP PARALLEL DO PRIVATE(myPhase,lattice_structure,neighboring_e,neighboring_i,neighboringPhase)
!$OMP neighboringInstance,neighboringStructure)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
myPhase = material_phase(1,i,e) ! get my phase myPhase = material_phase(1,i,e) ! get my phase
if (.not. phase_localPlasticity(myPhase)) then ! if nonlocal model if (.not. phase_localPlasticity(myPhase)) then ! if nonlocal model
myInstance = phase_plasticityInstance(myPhase)
! --- calculate disorientation between me and my neighbor --- ! --- calculate disorientation between me and my neighbor ---
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) ! loop through my neighbors do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) ! loop through my neighbors
@ -3375,13 +3366,11 @@ subroutine crystallite_orientations
if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists
neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase
if (.not. phase_localPlasticity(neighboringPhase)) then ! neighbor got also nonlocal plasticity if (.not. phase_localPlasticity(neighboringPhase)) then ! neighbor got also nonlocal plasticity
neighboringInstance = phase_plasticityInstance(neighboringPhase) if (lattice_structure(myPhase) == lattice_structure(neighboringPhase)) then ! if my neighbor has same crystal structure like me
neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure
if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me
crystallite_disorientation(:,n,1,i,e) = & crystallite_disorientation(:,n,1,i,e) = &
lattice_qDisorientation( crystallite_orientation(1:4,1,i,e), & lattice_qDisorientation( crystallite_orientation(1:4,1,i,e), &
crystallite_orientation(1:4,1,neighboring_i,neighboring_e), & crystallite_orientation(1:4,1,neighboring_i,neighboring_e), &
lattice_structure(myPhase)) ! calculate disorientation for given symmetry lattice_structure(myPhase)) ! calculate disorientation for given symmetry
else ! for neighbor with different phase else ! for neighbor with different phase
crystallite_disorientation(:,n,1,i,e) = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal] ! 180 degree rotation about 100 axis crystallite_disorientation(:,n,1,i,e) = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal] ! 180 degree rotation about 100 axis
endif endif