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

Merge remote-tracking branch 'origin/development' into YAML-improvements

This commit is contained in:
Sharan Roongta 2021-07-23 17:55:03 +02:00
parent ca8775e0d4 6c53615d11
commit c347688410
24 changed files with 735 additions and 813 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 8fec909d1931b092b223b0560dd30c3339c6e5a7 Subproject commit 174ecac2d3ab7596bdb60184d6bb9e1a52cb7378

1
examples/config/numerics.yaml
View File

@ -79,6 +79,5 @@ commercialFEM:
unitlength: 1 # physical length of one computational length unit unitlength: 1 # physical length of one computational length unit
generic: generic:
charLength: 1.0 # characteristic length scale for gradient problems.
random_seed: 0 # fixed seeding for pseudo-random number generator, Default 0: use random seed. random_seed: 0 # fixed seeding for pseudo-random number generator, Default 0: use random seed.
residualStiffness: 1.0e-6 # non-zero residual damage. residualStiffness: 1.0e-6 # non-zero residual damage.

2
examples/config/phase/damage/anisobrittle_cubic.yaml
View File

@ -7,5 +7,5 @@ q: 20
output: [f_phi] output: [f_phi]
K_11: 1.0 D_11: 1.0
mu: 0.001 mu: 0.001

2
examples/config/phase/damage/isobrittle_generic.yaml
View File

@ -3,5 +3,5 @@ W_crit: 1400000.0
output: [f_phi] output: [f_phi]
K_11: 1.0 D_11: 1.0
mu: 0.001 mu: 0.001

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha4-112-gb36da7378 v3.0.0-alpha4-137-gb69b85754

9
src/IO.f90
View File

@ -81,12 +81,7 @@ function IO_readlines(fileName) result(fileContent)
rawData = IO_read(fileName) rawData = IO_read(fileName)
!-------------------------------------------------------------------------------------------------- N_lines = count([(rawData(l:l) == IO_EOL,l=1,len(rawData))])
! count lines to allocate string array
N_lines = 0
do l=1, len(rawData)
if (rawData(l:l) == IO_EOL) N_lines = N_lines+1
enddo
allocate(fileContent(N_lines)) allocate(fileContent(N_lines))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -261,7 +256,7 @@ pure function IO_lc(string)
do i=1,len(string) do i=1,len(string)
n = index(UPPER,string(i:i)) n = index(UPPER,string(i:i))
if(n/=0) then if (n/=0) then
IO_lc(i:i) = LOWER(n:n) IO_lc(i:i) = LOWER(n:n)
else else
IO_lc(i:i) = string(i:i) IO_lc(i:i) = string(i:i)

9
src/grid/DAMASK_grid.f90
View File

@ -482,20 +482,13 @@ subroutine getMaskedTensor(values,mask,tensor)
values = 0.0 values = 0.0
if (tensor%length == 9) then ! temporary support for deprecated 1D tensor
do i = 1,9
mask((i-1)/3+1,mod(i-1,3)+1) = tensor%get_asString(i) /= 'x'
if (mask((i-1)/3+1,mod(i-1,3)+1)) values((i-1)/3+1,mod(i-1,3)+1) = tensor%get_asFloat(i)
enddo
else
do i = 1,3 do i = 1,3
row => tensor%get(i) row => tensor%get(i)
do j = 1,3 do j = 1,3
mask(i,j) = row%get_asString(j) /= 'x' ! ToDo change to np.masked behavior mask(i,j) = row%get_asString(j) /= 'x'
if (mask(i,j)) values(i,j) = row%get_asFloat(j) if (mask(i,j)) values(i,j) = row%get_asFloat(j)
enddo enddo
enddo enddo
endif
end subroutine end subroutine

9
src/homogenization.f90
View File

@ -41,15 +41,6 @@ module homogenization
integer(kind(DAMAGE_none_ID)), dimension(:), allocatable :: & integer(kind(DAMAGE_none_ID)), dimension(:), allocatable :: &
damage_type !< nonlocal damage model damage_type !< nonlocal damage model
type, private :: tNumerics_damage
real(pReal) :: &
charLength !< characteristic length scale for gradient problems
end type tNumerics_damage
type(tNumerics_damage), private :: &
num_damage
logical, public :: & logical, public :: &
terminallyIll = .false. !< at least one material point is terminally ill terminallyIll = .false. !< at least one material point is terminally ill

11
src/homogenization_damage.f90
View File

@ -37,8 +37,7 @@ module subroutine damage_init()
class(tNode), pointer :: & class(tNode), pointer :: &
configHomogenizations, & configHomogenizations, &
configHomogenization, & configHomogenization, &
configHomogenizationDamage, & configHomogenizationDamage
num_generic
integer :: ho,Nmembers integer :: ho,Nmembers
@ -70,11 +69,6 @@ module subroutine damage_init()
end associate end associate
enddo enddo
!------------------------------------------------------------------------------------
! read numerics parameter
num_generic => config_numerics%get('generic',defaultVal= emptyDict)
num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
call pass_init() call pass_init()
end subroutine damage_init end subroutine damage_init
@ -119,8 +113,7 @@ module function homogenization_K_phi(ce) result(K)
real(pReal), dimension(3,3) :: K real(pReal), dimension(3,3) :: K
K = phase_K_phi(1,ce) & K = phase_K_phi(1,ce)
* num_damage%charLength**2
end function homogenization_K_phi end function homogenization_K_phi

340
src/lattice.f90
View File

@ -3,7 +3,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief contains lattice structure definitions including Schmid matrices for slip, twin, trans, !> @brief contains lattice definitions including Schmid matrices for slip, twin, trans,
! and cleavage as well as interaction among the various systems ! and cleavage as well as interaction among the various systems
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module lattice module lattice
@ -365,12 +365,6 @@ module lattice
1, 1, 1, 1,-2, 1 & 1, 1, 1, 1,-2, 1 &
],pReal),shape(BCT_SYSTEMSLIP)) !< bct slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x) ],pReal),shape(BCT_SYSTEMSLIP)) !< bct slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x)
! SHOULD NOT BE PART OF LATTICE BEGIN
real(pReal), dimension(:), allocatable, public, protected :: &
lattice_mu, lattice_nu
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_C66
! SHOULD NOT BE PART OF LATTICE END
interface lattice_forestProjection_edge interface lattice_forestProjection_edge
module procedure slipProjection_transverse module procedure slipProjection_transverse
@ -384,7 +378,8 @@ module lattice
lattice_init, & lattice_init, &
lattice_equivalent_nu, & lattice_equivalent_nu, &
lattice_equivalent_mu, & lattice_equivalent_mu, &
lattice_applyLatticeSymmetry33, & lattice_symmetrize_33, &
lattice_symmetrize_C66, &
lattice_SchmidMatrix_slip, & lattice_SchmidMatrix_slip, &
lattice_SchmidMatrix_twin, & lattice_SchmidMatrix_twin, &
lattice_SchmidMatrix_trans, & lattice_SchmidMatrix_trans, &
@ -414,64 +409,20 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine lattice_init subroutine lattice_init
integer :: Nphases, ph,i
class(tNode), pointer :: &
phases, &
phase, &
mech, &
elasticity
print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT)
! SHOULD NOT BE PART OF LATTICE BEGIN
phases => config_material%get('phase')
Nphases = phases%length
allocate(lattice_C66(6,6,Nphases), source=0.0_pReal)
allocate(lattice_mu, lattice_nu,&
source=[(0.0_pReal,i=1,Nphases)])
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanical')
elasticity => mech%get('elastic')
lattice_C66(1,1,ph) = elasticity%get_asFloat('C_11')
lattice_C66(1,2,ph) = elasticity%get_asFloat('C_12')
lattice_C66(4,4,ph) = elasticity%get_asFloat('C_44')
lattice_C66(1,3,ph) = elasticity%get_asFloat('C_13',defaultVal=0.0_pReal)
lattice_C66(2,3,ph) = elasticity%get_asFloat('C_23',defaultVal=0.0_pReal)
lattice_C66(3,3,ph) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal)
lattice_C66(6,6,ph) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal)
lattice_C66(1:6,1:6,ph) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,ph),phase%get_asString('lattice'))
lattice_nu(ph) = lattice_equivalent_nu(lattice_C66(1:6,1:6,ph),'voigt')
lattice_mu(ph) = lattice_equivalent_mu(lattice_C66(1:6,1:6,ph),'voigt')
lattice_C66(1:6,1:6,ph) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,ph))) ! Literature data is in Voigt notation
do i = 1, 6
if (abs(lattice_C66(i,i,ph))<tol_math_check) &
call IO_error(135,el=i,ip=ph,ext_msg='matrix diagonal "el"ement of phase "ip"')
enddo
! SHOULD NOT BE PART OF LATTICE END
call selfTest call selfTest
enddo
end subroutine lattice_init end subroutine lattice_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Characteristic shear for twinning !> @brief Characteristic shear for twinning
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(characteristicShear) function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(sum(Ntwin)) :: characteristicShear real(pReal), dimension(sum(Ntwin)) :: characteristicShear
@ -513,7 +464,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
myFamilies: do f = 1,size(Ntwin,1) myFamilies: do f = 1,size(Ntwin,1)
mySystems: do s = 1,Ntwin(f) mySystems: do s = 1,Ntwin(f)
a = a + 1 a = a + 1
select case(structure) select case(lattice)
case('cF','cI') case('cF','cI')
characteristicShear(a) = 0.5_pReal*sqrt(2.0_pReal) characteristicShear(a) = 0.5_pReal*sqrt(2.0_pReal)
case('hP') case('hP')
@ -531,7 +482,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
characteristicShear(a) = 2.0_pReal*(cOverA**2.0_pReal-2.0_pReal)/3.0_pReal/cOverA characteristicShear(a) = 2.0_pReal*(cOverA**2.0_pReal-2.0_pReal)/3.0_pReal/cOverA
end select end select
case default case default
call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(structure)) call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(lattice))
end select end select
enddo mySystems enddo mySystems
enddo myFamilies enddo myFamilies
@ -542,10 +493,10 @@ end function lattice_characteristicShear_Twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Rotated elasticity matrices for twinning in 66-vector notation !> @brief Rotated elasticity matrices for twinning in 66-vector notation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_C66_twin(Ntwin,C66,structure,CoverA) function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), dimension(6,6), intent(in) :: C66 !< unrotated parent stiffness matrix real(pReal), dimension(6,6), intent(in) :: C66 !< unrotated parent stiffness matrix
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin
@ -554,18 +505,18 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA)
type(rotation) :: R type(rotation) :: R
integer :: i integer :: i
select case(structure) select case(lattice)
case('cF') case('cF')
coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,& coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,&
structure,0.0_pReal) lattice,0.0_pReal)
case('cI') case('cI')
coordinateSystem = buildCoordinateSystem(Ntwin,BCC_NSLIPSYSTEM,BCC_SYSTEMTWIN,& coordinateSystem = buildCoordinateSystem(Ntwin,BCC_NSLIPSYSTEM,BCC_SYSTEMTWIN,&
structure,0.0_pReal) lattice,0.0_pReal)
case('hP') case('hP')
coordinateSystem = buildCoordinateSystem(Ntwin,HEX_NSLIPSYSTEM,HEX_SYSTEMTWIN,& coordinateSystem = buildCoordinateSystem(Ntwin,HEX_NSLIPSYSTEM,HEX_SYSTEMTWIN,&
structure,cOverA) lattice,cOverA)
case default case default
call IO_error(137,ext_msg='lattice_C66_twin: '//trim(structure)) call IO_error(137,ext_msg='lattice_C66_twin: '//trim(lattice))
end select end select
do i = 1, sum(Ntwin) do i = 1, sum(Ntwin)
@ -579,11 +530,11 @@ end function lattice_C66_twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Rotated elasticity matrices for transformation in 66-vector notation !> @brief Rotated elasticity matrices for transformation in 66-vector notation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_C66_trans(Ntrans,C_parent66,structure_target, & function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
cOverA_trans,a_bcc,a_fcc) cOverA_trans,a_bcc,a_fcc)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=*), intent(in) :: structure_target !< lattice structure character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol)
real(pReal), dimension(6,6), intent(in) :: C_parent66 real(pReal), dimension(6,6), intent(in) :: C_parent66
real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans
@ -595,9 +546,9 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! elasticity matrix of the target phase in cube orientation ! elasticity matrix of the target phase in cube orientation
if (structure_target == 'hP') then if (lattice_target == 'hP') then
if (cOverA_trans < 1.0_pReal .or. cOverA_trans > 2.0_pReal) & if (cOverA_trans < 1.0_pReal .or. cOverA_trans > 2.0_pReal) &
call IO_error(131,ext_msg='lattice_C66_trans: '//trim(structure_target)) call IO_error(131,ext_msg='lattice_C66_trans: '//trim(lattice_target))
C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pReal*C_parent66(4,4))/2.0_pReal C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pReal*C_parent66(4,4))/2.0_pReal
C_bar66(1,2) = (C_parent66(1,1) + 5.0_pReal*C_parent66(1,2) - 2.0_pReal*C_parent66(4,4))/6.0_pReal C_bar66(1,2) = (C_parent66(1,1) + 5.0_pReal*C_parent66(1,2) - 2.0_pReal*C_parent66(4,4))/6.0_pReal
C_bar66(3,3) = (C_parent66(1,1) + 2.0_pReal*C_parent66(1,2) + 4.0_pReal*C_parent66(4,4))/3.0_pReal C_bar66(3,3) = (C_parent66(1,1) + 2.0_pReal*C_parent66(1,2) + 4.0_pReal*C_parent66(4,4))/3.0_pReal
@ -611,13 +562,13 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
C_target_unrotated66(1,3) = C_bar66(1,3) C_target_unrotated66(1,3) = C_bar66(1,3)
C_target_unrotated66(3,3) = C_bar66(3,3) C_target_unrotated66(3,3) = C_bar66(3,3)
C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2.0_pReal/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2))) C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2.0_pReal/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2)))
C_target_unrotated66 = applyLatticeSymmetryC66(C_target_unrotated66,'hP') C_target_unrotated66 = lattice_symmetrize_C66(C_target_unrotated66,'hP')
elseif (structure_target == 'cI') then elseif (lattice_target == 'cI') then
if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) & if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) &
call IO_error(134,ext_msg='lattice_C66_trans: '//trim(structure_target)) call IO_error(134,ext_msg='lattice_C66_trans: '//trim(lattice_target))
C_target_unrotated66 = C_parent66 C_target_unrotated66 = C_parent66
else else
call IO_error(137,ext_msg='lattice_C66_trans : '//trim(structure_target)) call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target))
endif endif
do i = 1, 6 do i = 1, 6
@ -686,11 +637,11 @@ end function lattice_nonSchmidMatrix
!> @brief Slip-slip interaction matrix !> @brief Slip-slip interaction matrix
!> details only active slip systems are considered !> details only active slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) result(interactionMatrix) function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), dimension(sum(Nslip),sum(Nslip)) :: interactionMatrix real(pReal), dimension(sum(Nslip),sum(Nslip)) :: interactionMatrix
integer, dimension(:), allocatable :: NslipMax integer, dimension(:), allocatable :: NslipMax
@ -908,7 +859,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
],shape(BCT_INTERACTIONSLIPSLIP)) ],shape(BCT_INTERACTIONSLIPSLIP))
select case(structure) select case(lattice)
case('cF') case('cF')
interactionTypes = FCC_INTERACTIONSLIPSLIP interactionTypes = FCC_INTERACTIONSLIPSLIP
NslipMax = FCC_NSLIPSYSTEM NslipMax = FCC_NSLIPSYSTEM
@ -922,7 +873,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
interactionTypes = BCT_INTERACTIONSLIPSLIP interactionTypes = BCT_INTERACTIONSLIPSLIP
NslipMax = BCT_NSLIPSYSTEM NslipMax = BCT_NSLIPSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(structure)) call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes)
@ -934,11 +885,11 @@ end function lattice_interaction_SlipBySlip
!> @brief Twin-twin interaction matrix !> @brief Twin-twin interaction matrix
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) result(interactionMatrix) function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: interactionMatrix real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: interactionMatrix
integer, dimension(:), allocatable :: NtwinMax integer, dimension(:), allocatable :: NtwinMax
@ -1009,7 +960,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16 & 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16 &
],shape(HEX_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for hex ],shape(HEX_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for hex
select case(structure) select case(lattice)
case('cF') case('cF')
interactionTypes = FCC_INTERACTIONTWINTWIN interactionTypes = FCC_INTERACTIONTWINTWIN
NtwinMax = FCC_NTWINSYSTEM NtwinMax = FCC_NTWINSYSTEM
@ -1020,7 +971,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul
interactionTypes = HEX_INTERACTIONTWINTWIN interactionTypes = HEX_INTERACTIONTWINTWIN
NtwinMax = HEX_NTWINSYSTEM NtwinMax = HEX_NTWINSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(structure)) call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes)
@ -1032,11 +983,11 @@ end function lattice_interaction_TwinByTwin
!> @brief Trans-trans interaction matrix !> @brief Trans-trans interaction matrix
!> details only active trans systems are considered !> details only active trans systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TransByTrans(Ntrans,interactionValues,structure) result(interactionMatrix) function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction
character(len=*), intent(in) :: structure !< lattice structure (parent crystal) character(len=2), intent(in) :: lattice !<Bravais lattice (Pearson symbol) (parent crystal)
real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: interactionMatrix real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: interactionMatrix
integer, dimension(:), allocatable :: NtransMax integer, dimension(:), allocatable :: NtransMax
@ -1058,11 +1009,11 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,structure) re
2,2,2,2,2,2,2,2,2,1,1,1 & 2,2,2,2,2,2,2,2,2,1,1,1 &
],shape(FCC_INTERACTIONTRANSTRANS)) !< Trans-trans interaction types for fcc ],shape(FCC_INTERACTIONTRANSTRANS)) !< Trans-trans interaction types for fcc
if(structure == 'cF') then if(lattice == 'cF') then
interactionTypes = FCC_INTERACTIONTRANSTRANS interactionTypes = FCC_INTERACTIONTRANSTRANS
NtransMax = FCC_NTRANSSYSTEM NtransMax = FCC_NTRANSSYSTEM
else else
call IO_error(137,ext_msg='lattice_interaction_TransByTrans: '//trim(structure)) call IO_error(137,ext_msg='lattice_interaction_TransByTrans: '//trim(lattice))
end if end if
interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes)
@ -1074,12 +1025,12 @@ end function lattice_interaction_TransByTrans
!> @brief Slip-twin interaction matrix !> @brief Slip-twin interaction matrix
!> details only active slip and twin systems are considered !> details only active slip and twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) result(interactionMatrix) function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntwin !< number of active twin systems per family Ntwin !< number of active twin systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-twin interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-twin interaction
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: interactionMatrix real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: interactionMatrix
integer, dimension(:), allocatable :: NslipMax, & integer, dimension(:), allocatable :: NslipMax, &
@ -1211,7 +1162,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure)
21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24 & 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24 &
],shape(HEX_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for hex ],shape(HEX_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for hex
select case(structure) select case(lattice)
case('cF') case('cF')
interactionTypes = FCC_INTERACTIONSLIPTWIN interactionTypes = FCC_INTERACTIONSLIPTWIN
NslipMax = FCC_NSLIPSYSTEM NslipMax = FCC_NSLIPSYSTEM
@ -1225,7 +1176,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure)
NslipMax = HEX_NSLIPSYSTEM NslipMax = HEX_NSLIPSYSTEM
NtwinMax = HEX_NTWINSYSTEM NtwinMax = HEX_NTWINSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(structure)) call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes)
@ -1237,12 +1188,12 @@ end function lattice_interaction_SlipByTwin
!> @brief Slip-trans interaction matrix !> @brief Slip-trans interaction matrix
!> details only active slip and trans systems are considered !> details only active slip and trans systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structure) result(interactionMatrix) function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntrans !< number of active trans systems per family Ntrans !< number of active trans systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-trans interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-trans interaction
character(len=*), intent(in) :: structure !< lattice structure (parent crystal) character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) (parent crystal)
real(pReal), dimension(sum(Nslip),sum(Ntrans)) :: interactionMatrix real(pReal), dimension(sum(Nslip),sum(Ntrans)) :: interactionMatrix
integer, dimension(:), allocatable :: NslipMax, & integer, dimension(:), allocatable :: NslipMax, &
@ -1272,13 +1223,13 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur
4,4,4,4,4,4,4,4,4,4,4,4 & 4,4,4,4,4,4,4,4,4,4,4,4 &
],shape(FCC_INTERACTIONSLIPTRANS)) !< Slip-trans interaction types for fcc ],shape(FCC_INTERACTIONSLIPTRANS)) !< Slip-trans interaction types for fcc
select case(structure) select case(lattice)
case('cF') case('cF')
interactionTypes = FCC_INTERACTIONSLIPTRANS interactionTypes = FCC_INTERACTIONSLIPTRANS
NslipMax = FCC_NSLIPSYSTEM NslipMax = FCC_NSLIPSYSTEM
NtransMax = FCC_NTRANSSYSTEM NtransMax = FCC_NTRANSSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(structure)) call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Nslip,Ntrans,NslipMax,NtransMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Ntrans,NslipMax,NtransMax,interactionValues,interactionTypes)
@ -1290,12 +1241,12 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur
!> @brief Twin-slip interaction matrix !> @brief Twin-slip interaction matrix
!> details only active twin and slip systems are considered !> details only active twin and slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) result(interactionMatrix) function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family
Nslip !< number of active slip systems per family Nslip !< number of active slip systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: interactionMatrix real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: interactionMatrix
integer, dimension(:), allocatable :: NtwinMax, & integer, dimension(:), allocatable :: NtwinMax, &
@ -1339,7 +1290,7 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure)
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24 & 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24 &
],shape(HEX_INTERACTIONTWINSLIP)) !< Twin-slip interaction types for hex ],shape(HEX_INTERACTIONTWINSLIP)) !< Twin-slip interaction types for hex
select case(structure) select case(lattice)
case('cF') case('cF')
interactionTypes = FCC_INTERACTIONTWINSLIP interactionTypes = FCC_INTERACTIONTWINSLIP
NtwinMax = FCC_NTWINSYSTEM NtwinMax = FCC_NTWINSYSTEM
@ -1353,7 +1304,7 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure)
NtwinMax = HEX_NTWINSYSTEM NtwinMax = HEX_NTWINSYSTEM
NslipMax = HEX_NSLIPSYSTEM NslipMax = HEX_NSLIPSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(structure)) call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes)
@ -1365,10 +1316,10 @@ end function lattice_interaction_TwinBySlip
!> @brief Schmid matrix for slip !> @brief Schmid matrix for slip
!> details only active slip systems are considered !> details only active slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA real(pReal), intent(in) :: cOverA
real(pReal), dimension(3,3,sum(Nslip)) :: SchmidMatrix real(pReal), dimension(3,3,sum(Nslip)) :: SchmidMatrix
@ -1377,7 +1328,7 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix)
integer, dimension(:), allocatable :: NslipMax integer, dimension(:), allocatable :: NslipMax
integer :: i integer :: i
select case(structure) select case(lattice)
case('cF') case('cF')
NslipMax = FCC_NSLIPSYSTEM NslipMax = FCC_NSLIPSYSTEM
slipSystems = FCC_SYSTEMSLIP slipSystems = FCC_SYSTEMSLIP
@ -1392,15 +1343,15 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix)
slipSystems = BCT_SYSTEMSLIP slipSystems = BCT_SYSTEMSLIP
case default case default
allocate(NslipMax(0)) allocate(NslipMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(structure)) call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(lattice))
end select end select
if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) &
call IO_error(145,ext_msg='Nslip '//trim(structure)) call IO_error(145,ext_msg='Nslip '//trim(lattice))
if (any(Nslip < 0)) & if (any(Nslip < 0)) &
call IO_error(144,ext_msg='Nslip '//trim(structure)) call IO_error(144,ext_msg='Nslip '//trim(lattice))
coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA) coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,lattice,cOverA)
do i = 1, sum(Nslip) do i = 1, sum(Nslip)
SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
@ -1415,10 +1366,10 @@ end function lattice_SchmidMatrix_slip
!> @brief Schmid matrix for twinning !> @brief Schmid matrix for twinning
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,3,sum(Ntwin)) :: SchmidMatrix real(pReal), dimension(3,3,sum(Ntwin)) :: SchmidMatrix
@ -1427,7 +1378,7 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
integer, dimension(:), allocatable :: NtwinMax integer, dimension(:), allocatable :: NtwinMax
integer :: i integer :: i
select case(structure) select case(lattice)
case('cF') case('cF')
NtwinMax = FCC_NTWINSYSTEM NtwinMax = FCC_NTWINSYSTEM
twinSystems = FCC_SYSTEMTWIN twinSystems = FCC_SYSTEMTWIN
@ -1439,15 +1390,15 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
twinSystems = HEX_SYSTEMTWIN twinSystems = HEX_SYSTEMTWIN
case default case default
allocate(NtwinMax(0)) allocate(NtwinMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(structure)) call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(lattice))
end select end select
if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) &
call IO_error(145,ext_msg='Ntwin '//trim(structure)) call IO_error(145,ext_msg='Ntwin '//trim(lattice))
if (any(Ntwin < 0)) & if (any(Ntwin < 0)) &
call IO_error(144,ext_msg='Ntwin '//trim(structure)) call IO_error(144,ext_msg='Ntwin '//trim(lattice))
coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure,cOverA) coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,lattice,cOverA)
do i = 1, sum(Ntwin) do i = 1, sum(Ntwin)
SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
@ -1462,24 +1413,24 @@ end function lattice_SchmidMatrix_twin
!> @brief Schmid matrix for twinning !> @brief Schmid matrix for twinning
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_trans(Ntrans,structure_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix) function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=*), intent(in) :: structure_target !< lattice structure character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,3,sum(Ntrans)) :: SchmidMatrix real(pReal), dimension(3,3,sum(Ntrans)) :: SchmidMatrix
real(pReal), dimension(3,3,sum(Ntrans)) :: devNull real(pReal), dimension(3,3,sum(Ntrans)) :: devNull
real(pReal) :: a_bcc, a_fcc real(pReal) :: a_bcc, a_fcc
if (structure_target /= 'cI' .and. structure_target /= 'hP') & if (lattice_target /= 'cI' .and. lattice_target /= 'hP') &
call IO_error(137,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target)) call IO_error(137,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
if (structure_target == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & if (lattice_target == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) &
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target)) call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
if (structure_target == 'cI' .and. (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal)) & if (lattice_target == 'cI' .and. (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal)) &
call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target)) call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA,a_fcc,a_bcc) call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA,a_fcc,a_bcc)
@ -1490,10 +1441,10 @@ end function lattice_SchmidMatrix_trans
!> @brief Schmid matrix for cleavage !> @brief Schmid matrix for cleavage
!> details only active cleavage systems are considered !> details only active cleavage systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(SchmidMatrix) function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,3,3,sum(Ncleavage)) :: SchmidMatrix real(pReal), dimension(3,3,3,sum(Ncleavage)) :: SchmidMatrix
@ -1502,7 +1453,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid
integer, dimension(:), allocatable :: NcleavageMax integer, dimension(:), allocatable :: NcleavageMax
integer :: i integer :: i
select case(structure) select case(lattice)
case('cF') case('cF')
NcleavageMax = FCC_NCLEAVAGESYSTEM NcleavageMax = FCC_NCLEAVAGESYSTEM
cleavageSystems = FCC_SYSTEMCLEAVAGE cleavageSystems = FCC_SYSTEMCLEAVAGE
@ -1511,15 +1462,15 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid
cleavageSystems = BCC_SYSTEMCLEAVAGE cleavageSystems = BCC_SYSTEMCLEAVAGE
case default case default
allocate(NcleavageMax(0)) allocate(NcleavageMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(structure)) call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(lattice))
end select end select
if (any(NcleavageMax(1:size(Ncleavage)) - Ncleavage < 0)) & if (any(NcleavageMax(1:size(Ncleavage)) - Ncleavage < 0)) &
call IO_error(145,ext_msg='Ncleavage '//trim(structure)) call IO_error(145,ext_msg='Ncleavage '//trim(lattice))
if (any(Ncleavage < 0)) & if (any(Ncleavage < 0)) &
call IO_error(144,ext_msg='Ncleavage '//trim(structure)) call IO_error(144,ext_msg='Ncleavage '//trim(lattice))
coordinateSystem = buildCoordinateSystem(Ncleavage,NcleavageMax,cleavageSystems,structure,cOverA) coordinateSystem = buildCoordinateSystem(Ncleavage,NcleavageMax,cleavageSystems,lattice,cOverA)
do i = 1, sum(Ncleavage) do i = 1, sum(Ncleavage)
SchmidMatrix(1:3,1:3,1,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,1,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
@ -1533,16 +1484,16 @@ end function lattice_SchmidMatrix_cleavage
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Slip direction of slip systems (|| b) !> @brief Slip direction of slip systems (|| b)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_slip_direction(Nslip,structure,cOverA) result(d) function lattice_slip_direction(Nslip,lattice,cOverA) result(d)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,sum(Nslip)) :: d real(pReal), dimension(3,sum(Nslip)) :: d
real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem
coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA) coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
d = coordinateSystem(1:3,1,1:sum(Nslip)) d = coordinateSystem(1:3,1,1:sum(Nslip))
end function lattice_slip_direction end function lattice_slip_direction
@ -1551,16 +1502,16 @@ end function lattice_slip_direction
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Normal direction of slip systems (|| n) !> @brief Normal direction of slip systems (|| n)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_slip_normal(Nslip,structure,cOverA) result(n) function lattice_slip_normal(Nslip,lattice,cOverA) result(n)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,sum(Nslip)) :: n real(pReal), dimension(3,sum(Nslip)) :: n
real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem
coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA) coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
n = coordinateSystem(1:3,2,1:sum(Nslip)) n = coordinateSystem(1:3,2,1:sum(Nslip))
end function lattice_slip_normal end function lattice_slip_normal
@ -1569,16 +1520,16 @@ end function lattice_slip_normal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Transverse direction of slip systems (|| t = b x n) !> @brief Transverse direction of slip systems (|| t = b x n)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_slip_transverse(Nslip,structure,cOverA) result(t) function lattice_slip_transverse(Nslip,lattice,cOverA) result(t)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,sum(Nslip)) :: t real(pReal), dimension(3,sum(Nslip)) :: t
real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem
coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA) coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
t = coordinateSystem(1:3,3,1:sum(Nslip)) t = coordinateSystem(1:3,3,1:sum(Nslip))
end function lattice_slip_transverse end function lattice_slip_transverse
@ -1588,17 +1539,17 @@ end function lattice_slip_transverse
!> @brief Labels for slip systems !> @brief Labels for slip systems
!> details only active slip systems are considered !> details only active slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_labels_slip(Nslip,structure) result(labels) function lattice_labels_slip(Nslip,lattice) result(labels)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
character(len=:), dimension(:), allocatable :: labels character(len=:), dimension(:), allocatable :: labels
real(pReal), dimension(:,:), allocatable :: slipSystems real(pReal), dimension(:,:), allocatable :: slipSystems
integer, dimension(:), allocatable :: NslipMax integer, dimension(:), allocatable :: NslipMax
select case(structure) select case(lattice)
case('cF') case('cF')
NslipMax = FCC_NSLIPSYSTEM NslipMax = FCC_NSLIPSYSTEM
slipSystems = FCC_SYSTEMSLIP slipSystems = FCC_SYSTEMSLIP
@ -1612,13 +1563,13 @@ function lattice_labels_slip(Nslip,structure) result(labels)
NslipMax = BCT_NSLIPSYSTEM NslipMax = BCT_NSLIPSYSTEM
slipSystems = BCT_SYSTEMSLIP slipSystems = BCT_SYSTEMSLIP
case default case default
call IO_error(137,ext_msg='lattice_labels_slip: '//trim(structure)) call IO_error(137,ext_msg='lattice_labels_slip: '//trim(lattice))
end select end select
if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) &
call IO_error(145,ext_msg='Nslip '//trim(structure)) call IO_error(145,ext_msg='Nslip '//trim(lattice))
if (any(Nslip < 0)) & if (any(Nslip < 0)) &
call IO_error(144,ext_msg='Nslip '//trim(structure)) call IO_error(144,ext_msg='Nslip '//trim(lattice))
labels = getLabels(Nslip,NslipMax,slipSystems) labels = getLabels(Nslip,NslipMax,slipSystems)
@ -1626,19 +1577,19 @@ end function lattice_labels_slip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Return 3x3 tensor with symmetry according to given crystal structure !> @brief Return 3x3 tensor with symmetry according to given Bravais lattice
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_applyLatticeSymmetry33(T,structure) result(T_sym) pure function lattice_symmetrize_33(T,lattice) result(T_sym)
real(pReal), dimension(3,3) :: T_sym real(pReal), dimension(3,3) :: T_sym
real(pReal), dimension(3,3), intent(in) :: T real(pReal), dimension(3,3), intent(in) :: T
character(len=*), intent(in) :: structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
T_sym = 0.0_pReal T_sym = 0.0_pReal
select case(structure) select case(lattice)
case('cF','cI') case('cF','cI')
T_sym(1,1) = T(1,1) T_sym(1,1) = T(1,1)
T_sym(2,2) = T(1,1) T_sym(2,2) = T(1,1)
@ -1649,26 +1600,26 @@ pure function lattice_applyLatticeSymmetry33(T,structure) result(T_sym)
T_sym(3,3) = T(3,3) T_sym(3,3) = T(3,3)
end select end select
end function lattice_applyLatticeSymmetry33 end function lattice_symmetrize_33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Return stiffness matrix in 6x6 notation with symmetry according to given crystal structure !> @brief Return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice
!> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 !> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function applyLatticeSymmetryC66(C66,structure) result(C66_sym) pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym)
real(pReal), dimension(6,6) :: C66_sym real(pReal), dimension(6,6) :: C66_sym
real(pReal), dimension(6,6), intent(in) :: C66 real(pReal), dimension(6,6), intent(in) :: C66
character(len=*), intent(in) :: structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
integer :: i,j integer :: i,j
C66_sym = 0.0_pReal C66_sym = 0.0_pReal
select case(structure) select case(lattice)
case ('cF','cI') case ('cF','cI')
C66_sym(1,1) = C66(1,1); C66_sym(2,2) = C66(1,1); C66_sym(3,3) = C66(1,1) C66_sym(1,1) = C66(1,1); C66_sym(2,2) = C66(1,1); C66_sym(3,3) = C66(1,1)
C66_sym(1,2) = C66(1,2); C66_sym(1,3) = C66(1,2); C66_sym(2,3) = C66(1,2) C66_sym(1,2) = C66(1,2); C66_sym(1,3) = C66(1,2); C66_sym(2,3) = C66(1,2)
@ -1695,24 +1646,24 @@ pure function applyLatticeSymmetryC66(C66,structure) result(C66_sym)
enddo enddo
enddo enddo
end function applyLatticeSymmetryC66 end function lattice_symmetrize_C66
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Labels for twin systems !> @brief Labels for twin systems
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_labels_twin(Ntwin,structure) result(labels) function lattice_labels_twin(Ntwin,lattice) result(labels)
integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
character(len=:), dimension(:), allocatable :: labels character(len=:), dimension(:), allocatable :: labels
real(pReal), dimension(:,:), allocatable :: twinSystems real(pReal), dimension(:,:), allocatable :: twinSystems
integer, dimension(:), allocatable :: NtwinMax integer, dimension(:), allocatable :: NtwinMax
select case(structure) select case(lattice)
case('cF') case('cF')
NtwinMax = FCC_NTWINSYSTEM NtwinMax = FCC_NTWINSYSTEM
twinSystems = FCC_SYSTEMTWIN twinSystems = FCC_SYSTEMTWIN
@ -1723,13 +1674,13 @@ function lattice_labels_twin(Ntwin,structure) result(labels)
NtwinMax = HEX_NTWINSYSTEM NtwinMax = HEX_NTWINSYSTEM
twinSystems = HEX_SYSTEMTWIN twinSystems = HEX_SYSTEMTWIN
case default case default
call IO_error(137,ext_msg='lattice_labels_twin: '//trim(structure)) call IO_error(137,ext_msg='lattice_labels_twin: '//trim(lattice))
end select end select
if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) &
call IO_error(145,ext_msg='Ntwin '//trim(structure)) call IO_error(145,ext_msg='Ntwin '//trim(lattice))
if (any(Ntwin < 0)) & if (any(Ntwin < 0)) &
call IO_error(144,ext_msg='Ntwin '//trim(structure)) call IO_error(144,ext_msg='Ntwin '//trim(lattice))
labels = getLabels(Ntwin,NtwinMax,twinSystems) labels = getLabels(Ntwin,NtwinMax,twinSystems)
@ -1740,18 +1691,18 @@ end function lattice_labels_twin
!> @brief Projection of the transverse direction onto the slip plane !> @brief Projection of the transverse direction onto the slip plane
!> @details: This projection is used to calculate forest hardening for edge dislocations !> @details: This projection is used to calculate forest hardening for edge dislocations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function slipProjection_transverse(Nslip,structure,cOverA) result(projection) function slipProjection_transverse(Nslip,lattice,cOverA) result(projection)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection
real(pReal), dimension(3,sum(Nslip)) :: n, t real(pReal), dimension(3,sum(Nslip)) :: n, t
integer :: i, j integer :: i, j
n = lattice_slip_normal (Nslip,structure,cOverA) n = lattice_slip_normal (Nslip,lattice,cOverA)
t = lattice_slip_transverse(Nslip,structure,cOverA) t = lattice_slip_transverse(Nslip,lattice,cOverA)
do i=1, sum(Nslip); do j=1, sum(Nslip) do i=1, sum(Nslip); do j=1, sum(Nslip)
projection(i,j) = abs(math_inner(n(:,i),t(:,j))) projection(i,j) = abs(math_inner(n(:,i),t(:,j)))
@ -1764,18 +1715,18 @@ end function slipProjection_transverse
!> @brief Projection of the slip direction onto the slip plane !> @brief Projection of the slip direction onto the slip plane
!> @details: This projection is used to calculate forest hardening for screw dislocations !> @details: This projection is used to calculate forest hardening for screw dislocations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function slipProjection_direction(Nslip,structure,cOverA) result(projection) function slipProjection_direction(Nslip,lattice,cOverA) result(projection)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection
real(pReal), dimension(3,sum(Nslip)) :: n, d real(pReal), dimension(3,sum(Nslip)) :: n, d
integer :: i, j integer :: i, j
n = lattice_slip_normal (Nslip,structure,cOverA) n = lattice_slip_normal (Nslip,lattice,cOverA)
d = lattice_slip_direction(Nslip,structure,cOverA) d = lattice_slip_direction(Nslip,lattice,cOverA)
do i=1, sum(Nslip); do j=1, sum(Nslip) do i=1, sum(Nslip); do j=1, sum(Nslip)
projection(i,j) = abs(math_inner(n(:,i),d(:,j))) projection(i,j) = abs(math_inner(n(:,i),d(:,j)))
@ -1788,17 +1739,17 @@ end function slipProjection_direction
!> @brief build a local coordinate system on slip systems !> @brief build a local coordinate system on slip systems
!> @details Order: Direction, plane (normal), and common perpendicular !> @details Order: Direction, plane (normal), and common perpendicular
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem) function coordinateSystem_slip(Nslip,lattice,cOverA) result(coordinateSystem)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem
real(pReal), dimension(:,:), allocatable :: slipSystems real(pReal), dimension(:,:), allocatable :: slipSystems
integer, dimension(:), allocatable :: NslipMax integer, dimension(:), allocatable :: NslipMax
select case(structure) select case(lattice)
case('cF') case('cF')
NslipMax = FCC_NSLIPSYSTEM NslipMax = FCC_NSLIPSYSTEM
slipSystems = FCC_SYSTEMSLIP slipSystems = FCC_SYSTEMSLIP
@ -1813,15 +1764,15 @@ function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem)
slipSystems = BCT_SYSTEMSLIP slipSystems = BCT_SYSTEMSLIP
case default case default
allocate(NslipMax(0)) allocate(NslipMax(0))
call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(structure)) call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(lattice))
end select end select
if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) &
call IO_error(145,ext_msg='Nslip '//trim(structure)) call IO_error(145,ext_msg='Nslip '//trim(lattice))
if (any(Nslip < 0)) & if (any(Nslip < 0)) &
call IO_error(144,ext_msg='Nslip '//trim(structure)) call IO_error(144,ext_msg='Nslip '//trim(lattice))
coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA) coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,lattice,cOverA)
end function coordinateSystem_slip end function coordinateSystem_slip
@ -1873,15 +1824,15 @@ end function buildInteraction
!> @brief Build a local coordinate system on slip, twin, trans, cleavage systems !> @brief Build a local coordinate system on slip, twin, trans, cleavage systems
!> @details Order: Direction, plane (normal), and common perpendicular !> @details Order: Direction, plane (normal), and common perpendicular
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function buildCoordinateSystem(active,potential,system,structure,cOverA) function buildCoordinateSystem(active,potential,system,lattice,cOverA)
integer, dimension(:), intent(in) :: & integer, dimension(:), intent(in) :: &
active, & !< # of active systems per family active, & !< # of active systems per family
potential !< # of potential systems per family potential !< # of potential systems per family
real(pReal), dimension(:,:), intent(in) :: & real(pReal), dimension(:,:), intent(in) :: &
system system
character(len=*), intent(in) :: & character(len=2), intent(in) :: &
structure !< lattice structure lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
cOverA cOverA
real(pReal), dimension(3,3,sum(active)) :: & real(pReal), dimension(3,3,sum(active)) :: &
@ -1895,10 +1846,10 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA)
f, & !< index of my family f, & !< index of my family
s !< index of my system in current family s !< index of my system in current family
if (structure == 'tI' .and. cOverA > 2.0_pReal) & if (lattice == 'tI' .and. cOverA > 2.0_pReal) &
call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure)) call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(lattice))
if (structure == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & if (lattice == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) &
call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure)) call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(lattice))
a = 0 a = 0
activeFamilies: do f = 1,size(active,1) activeFamilies: do f = 1,size(active,1)
@ -1906,7 +1857,7 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA)
a = a + 1 a = a + 1
p = sum(potential(1:f-1))+s p = sum(potential(1:f-1))+s
select case(structure) select case(lattice)
case ('cF','cI','tI') case ('cF','cI','tI')
direction = system(1:3,p) direction = system(1:3,p)
@ -1921,7 +1872,7 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA)
system(8,p)/cOverA ] ! plane (hkil)->(h (h+2k)/sqrt(3) l/(p/a)) system(8,p)/cOverA ] ! plane (hkil)->(h (h+2k)/sqrt(3) l/(p/a))
case default case default
call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(structure)) call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(lattice))
end select end select
@ -1950,9 +1901,9 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
Q, & !< Total rotation: Q = R*B Q, & !< Total rotation: Q = R*B
S !< Eigendeformation tensor for phase transformation S !< Eigendeformation tensor for phase transformation
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
cOverA, & !< c/a for target hex structure cOverA, & !< c/a for target hex lattice
a_bcc, & !< lattice parameter a for target bcc structure a_bcc, & !< lattice parameter a for bcc target lattice
a_fcc !< lattice parameter a for parent fcc structure a_fcc !< lattice parameter a for fcc parent lattice
type(rotation) :: & type(rotation) :: &
R, & !< Pitsch rotation R, & !< Pitsch rotation
@ -2126,9 +2077,10 @@ end function getlabels
function lattice_equivalent_nu(C,assumption) result(nu) function lattice_equivalent_nu(C,assumption) result(nu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
real(pReal) :: K, mu, nu real(pReal) :: nu
real(pReal) :: K, mu
logical :: error logical :: error
real(pReal), dimension(6,6) :: S real(pReal), dimension(6,6) :: S
@ -2158,7 +2110,7 @@ end function lattice_equivalent_nu
function lattice_equivalent_mu(C,assumption) result(mu) function lattice_equivalent_mu(C,assumption) result(mu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
real(pReal) :: mu real(pReal) :: mu
logical :: error logical :: error
@ -2203,10 +2155,10 @@ subroutine selfTest
do i = 1, 10 do i = 1, 10
call random_number(C) call random_number(C)
C_cF = applyLatticeSymmetryC66(C,'cI') C_cF = lattice_symmetrize_C66(C,'cI')
C_cI = applyLatticeSymmetryC66(C,'cF') C_cI = lattice_symmetrize_C66(C,'cF')
C_hP = applyLatticeSymmetryC66(C,'hP') C_hP = lattice_symmetrize_C66(C,'hP')
C_tI = applyLatticeSymmetryC66(C,'tI') C_tI = lattice_symmetrize_C66(C,'tI')
if (any(dNeq(C_cI,transpose(C_cF)))) error stop 'SymmetryC66/cI-cF' if (any(dNeq(C_cI,transpose(C_cF)))) error stop 'SymmetryC66/cI-cF'
if (any(dNeq(C_cF,transpose(C_cI)))) error stop 'SymmetryC66/cF-cI' if (any(dNeq(C_cF,transpose(C_cI)))) error stop 'SymmetryC66/cF-cI'
@ -2226,10 +2178,10 @@ subroutine selfTest
if (any(dNeq(C(4,4),[C_tI(4,4),C_tI(5,5)]))) error stop 'SymmetryC_44-55/tI' if (any(dNeq(C(4,4),[C_tI(4,4),C_tI(5,5)]))) error stop 'SymmetryC_44-55/tI'
call random_number(T) call random_number(T)
T_cF = lattice_applyLatticeSymmetry33(T,'cI') T_cF = lattice_symmetrize_33(T,'cI')
T_cI = lattice_applyLatticeSymmetry33(T,'cF') T_cI = lattice_symmetrize_33(T,'cF')
T_hP = lattice_applyLatticeSymmetry33(T,'hP') T_hP = lattice_symmetrize_33(T,'hP')
T_tI = lattice_applyLatticeSymmetry33(T,'tI') T_tI = lattice_symmetrize_33(T,'tI')
if (any(dNeq0(T_cF) .and. math_I3<1.0_pReal)) error stop 'Symmetry33/c' if (any(dNeq0(T_cF) .and. math_I3<1.0_pReal)) error stop 'Symmetry33/c'
if (any(dNeq0(T_hP) .and. math_I3<1.0_pReal)) error stop 'Symmetry33/hP' if (any(dNeq0(T_hP) .and. math_I3<1.0_pReal)) error stop 'Symmetry33/hP'
@ -2244,7 +2196,7 @@ subroutine selfTest
call random_number(C) call random_number(C)
C(1,1) = C(1,1) + C(1,2) + 0.1_pReal C(1,1) = C(1,1) + C(1,2) + 0.1_pReal
C(4,4) = 0.5_pReal * (C(1,1) - C(1,2)) C(4,4) = 0.5_pReal * (C(1,1) - C(1,2))
C = applyLatticeSymmetryC66(C,'cI') C = lattice_symmetrize_C66(C,'cI')
if(dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt' if(dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
if(dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss' if(dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'

16
src/math.f90
View File

@ -1068,6 +1068,7 @@ integer pure function math_factorial(n)
integer, intent(in) :: n integer, intent(in) :: n
math_factorial = product(math_range(n)) math_factorial = product(math_range(n))
end function math_factorial end function math_factorial
@ -1081,6 +1082,7 @@ integer pure function math_binomial(n,k)
integer, intent(in) :: n, k integer, intent(in) :: n, k
integer :: i, k_, n_ integer :: i, k_, n_
k_ = min(k,n-k) k_ = min(k,n-k)
n_ = n n_ = n
math_binomial = merge(1,0,k_>-1) ! handling special cases k < 0 or k > n math_binomial = merge(1,0,k_>-1) ! handling special cases k < 0 or k > n
@ -1095,15 +1097,13 @@ end function math_binomial
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief multinomial coefficient !> @brief multinomial coefficient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
integer pure function math_multinomial(alpha) integer pure function math_multinomial(k)
integer, intent(in), dimension(:) :: alpha integer, intent(in), dimension(:) :: k
integer :: i integer :: i
math_multinomial = 1
do i = 1, size(alpha) math_multinomial = product([(math_binomial(sum(k(1:i)),k(i)),i=1,size(k))])
math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i))
enddo
end function math_multinomial end function math_multinomial
@ -1116,6 +1116,7 @@ real(pReal) pure function math_volTetrahedron(v1,v2,v3,v4)
real(pReal), dimension (3), intent(in) :: v1,v2,v3,v4 real(pReal), dimension (3), intent(in) :: v1,v2,v3,v4
real(pReal), dimension (3,3) :: m real(pReal), dimension (3,3) :: m
m(1:3,1) = v1-v2 m(1:3,1) = v1-v2
m(1:3,2) = v1-v3 m(1:3,2) = v1-v3
m(1:3,3) = v1-v4 m(1:3,3) = v1-v4
@ -1289,6 +1290,9 @@ subroutine selfTest
if(math_binomial(49,6) /= 13983816) & if(math_binomial(49,6) /= 13983816) &
error stop 'math_binomial' error stop 'math_binomial'
if(math_multinomial([1,2,3,4]) /= 12600) &
error stop 'math_multinomial'
ijk = cshift([1,2,3],int(r*1.0e2_pReal)) ijk = cshift([1,2,3],int(r*1.0e2_pReal))
if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) & if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
error stop 'math_LeviCivita(even)' error stop 'math_LeviCivita(even)'

18
src/phase_damage.f90
View File

@ -5,7 +5,7 @@ submodule(phase) damage
type :: tDamageParameters type :: tDamageParameters
real(pReal) :: mu = 0.0_pReal !< viscosity real(pReal) :: mu = 0.0_pReal !< viscosity
real(pReal), dimension(3,3) :: K = 0.0_pReal !< conductivity/diffusivity real(pReal), dimension(3,3) :: D = 0.0_pReal !< conductivity/diffusivity
end type tDamageParameters end type tDamageParameters
enum, bind(c); enumerator :: & enum, bind(c); enumerator :: &
@ -18,7 +18,7 @@ submodule(phase) damage
type :: tDataContainer type :: tDataContainer
real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi real(pReal), dimension(:), allocatable :: phi
end type tDataContainer end type tDataContainer
integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: &
@ -97,7 +97,6 @@ module subroutine damage_init
Nmembers = count(material_phaseID == ph) Nmembers = count(material_phaseID == ph)
allocate(current(ph)%phi(Nmembers),source=1.0_pReal) allocate(current(ph)%phi(Nmembers),source=1.0_pReal)
allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal)
phase => phases%get(ph) phase => phases%get(ph)
sources => phase%get('damage',defaultVal=emptyList) sources => phase%get('damage',defaultVal=emptyList)
@ -105,10 +104,10 @@ module subroutine damage_init
if (sources%length == 1) then if (sources%length == 1) then
damage_active = .true. damage_active = .true.
source => sources%get(1) source => sources%get(1)
param(ph)%mu = source%get_asFloat('mu',defaultVal=0.0_pReal) ! ToDo: make mandatory? param(ph)%mu = source%get_asFloat('mu')
param(ph)%K(1,1) = source%get_asFloat('K_11',defaultVal=0.0_pReal) ! ToDo: make mandatory? param(ph)%D(1,1) = source%get_asFloat('D_11')
param(ph)%K(3,3) = source%get_asFloat('K_33',defaultVal=0.0_pReal) ! ToDo: depends on symmetry if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%D(3,3) = source%get_asFloat('D_33')
param(ph)%K = lattice_applyLatticeSymmetry33(param(ph)%K,phase_lattice(ph)) param(ph)%D = lattice_symmetrize_33(param(ph)%D,phase_lattice(ph))
endif endif
enddo enddo
@ -385,9 +384,10 @@ module function phase_K_phi(co,ce) result(K)
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: K real(pReal), dimension(3,3) :: K
real(pReal), parameter :: l = 1.0_pReal
K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%D) \
K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%K) * l**2.0_pReal
end function phase_K_phi end function phase_K_phi

3
src/phase_damage_anisobrittle.f90
View File

@ -73,8 +73,7 @@ module function anisobrittle_init() result(mySources)
prm%s_crit = src%get_as1dFloat('s_crit', requiredSize=size(N_cl)) prm%s_crit = src%get_as1dFloat('s_crit', requiredSize=size(N_cl))
prm%g_crit = src%get_as1dFloat('g_crit', requiredSize=size(N_cl)) prm%g_crit = src%get_as1dFloat('g_crit', requiredSize=size(N_cl))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),& prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
! expand: family => system ! expand: family => system
prm%s_crit = math_expand(prm%s_crit,N_cl) prm%s_crit = math_expand(prm%s_crit,N_cl)

14
src/phase_mechanical.f90
View File

@ -168,6 +168,20 @@ submodule(phase) mechanical
integer, intent(in) :: ph,en integer, intent(in) :: ph,en
end function plastic_dislotwin_homogenizedC end function plastic_dislotwin_homogenizedC
module function elastic_C66(ph) result(C66)
real(pReal), dimension(6,6) :: C66
integer, intent(in) :: ph
end function elastic_C66
module function elastic_mu(ph) result(mu)
real(pReal) :: mu
integer, intent(in) :: ph
end function elastic_mu
module function elastic_nu(ph) result(nu)
real(pReal) :: nu
integer, intent(in) :: ph
end function elastic_nu
end interface end interface
type :: tOutput !< new requested output (per phase) type :: tOutput !< new requested output (per phase)

18
src/phase_mechanical_eigen.f90
View File

@ -86,17 +86,17 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki
kinematics, & kinematics, &
kinematics_type, & kinematics_type, &
mechanics mechanics
integer :: p,k integer :: ph,k
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(active_kinematics(kinematics_length,phases%length), source = .false. ) allocate(active_kinematics(kinematics_length,phases%length), source = .false. )
do p = 1, phases%length do ph = 1, phases%length
phase => phases%get(p) phase => phases%get(ph)
mechanics => phase%get('mechanical') mechanics => phase%get('mechanical')
kinematics => mechanics%get('eigen',defaultVal=emptyList) kinematics => mechanics%get('eigen',defaultVal=emptyList)
do k = 1, kinematics%length do k = 1, kinematics%length
kinematics_type => kinematics%get(k) kinematics_type => kinematics%get(k)
active_kinematics(k,p) = kinematics_type%get_asString('type') == kinematics_label active_kinematics(k,ph) = kinematics_type%get_asString('type') == kinematics_label
enddo enddo
enddo enddo
@ -118,17 +118,17 @@ function kinematics_active2(kinematics_label) result(active_kinematics)
phase, & phase, &
kinematics, & kinematics, &
kinematics_type kinematics_type
integer :: p integer :: ph
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(active_kinematics(phases%length), source = .false. ) allocate(active_kinematics(phases%length), source = .false.)
do p = 1, phases%length do ph = 1, phases%length
phase => phases%get(p) phase => phases%get(ph)
kinematics => phase%get('damage',defaultVal=emptyList) kinematics => phase%get('damage',defaultVal=emptyList)
if(kinematics%length < 1) return if(kinematics%length < 1) return
kinematics_type => kinematics%get(1) kinematics_type => kinematics%get(1)
if (.not. kinematics_type%contains('type')) continue if (.not. kinematics_type%contains('type')) continue
active_kinematics(p) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label active_kinematics(ph) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label
enddo enddo

2
src/phase_mechanical_eigen_thermalexpansion.f90
View File

@ -69,7 +69,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal) prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal)
endif endif
do i=1, size(prm%A,3) do i=1, size(prm%A,3)
prm%A(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%A(1:3,1:3,i),phase_lattice(p)) prm%A(1:3,1:3,i) = lattice_symmetrize_33(prm%A(1:3,1:3,i),phase_lattice(p))
enddo enddo
end associate end associate
endif endif

44
src/phase_mechanical_elastic.f90
View File

@ -2,7 +2,10 @@ submodule(phase:mechanical) elastic
type :: tParameters type :: tParameters
real(pReal), dimension(6,6) :: & real(pReal), dimension(6,6) :: &
C66 !< Elastic constants in Voig notation C66 = 0.0_pReal !< Elastic constants in Voigt notation
real(pReal) :: &
mu, &
nu
end type tParameters end type tParameters
type(tParameters), allocatable, dimension(:) :: param type(tParameters), allocatable, dimension(:) :: param
@ -37,6 +40,7 @@ module subroutine elastic_init(phases)
if (elastic%get_asString('type') /= 'Hooke') call IO_error(200,ext_msg=elastic%get_asString('type')) if (elastic%get_asString('type') /= 'Hooke') call IO_error(200,ext_msg=elastic%get_asString('type'))
associate(prm => param(ph)) associate(prm => param(ph))
prm%C66(1,1) = elastic%get_asFloat('C_11') prm%C66(1,1) = elastic%get_asFloat('C_11')
prm%C66(1,2) = elastic%get_asFloat('C_12') prm%C66(1,2) = elastic%get_asFloat('C_12')
prm%C66(4,4) = elastic%get_asFloat('C_44') prm%C66(4,4) = elastic%get_asFloat('C_44')
@ -46,6 +50,14 @@ module subroutine elastic_init(phases)
prm%C66(3,3) = elastic%get_asFloat('C_33') prm%C66(3,3) = elastic%get_asFloat('C_33')
endif endif
if (phase_lattice(ph) == 'tI') prm%C66(6,6) = elastic%get_asFloat('C_66') if (phase_lattice(ph) == 'tI') prm%C66(6,6) = elastic%get_asFloat('C_66')
prm%C66 = lattice_symmetrize_C66(prm%C66,phase_lattice(ph))
prm%nu = lattice_equivalent_nu(prm%C66,'voigt')
prm%mu = lattice_equivalent_mu(prm%C66,'voigt')
prm%C66 = math_sym3333to66(math_Voigt66to3333(prm%C66)) ! Literature data is in Voigt notation
end associate end associate
enddo enddo
@ -104,10 +116,38 @@ module function phase_homogenizedC(ph,en) result(C)
case (PLASTICITY_DISLOTWIN_ID) plasticType case (PLASTICITY_DISLOTWIN_ID) plasticType
C = plastic_dislotwin_homogenizedC(ph,en) C = plastic_dislotwin_homogenizedC(ph,en)
case default plasticType case default plasticType
C = lattice_C66(1:6,1:6,ph) C = param(ph)%C66
end select plasticType end select plasticType
end function phase_homogenizedC end function phase_homogenizedC
module function elastic_C66(ph) result(C66)
real(pReal), dimension(6,6) :: C66
integer, intent(in) :: ph
C66 = param(ph)%C66
end function elastic_C66
module function elastic_mu(ph) result(mu)
real(pReal) :: mu
integer, intent(in) :: ph
mu = param(ph)%mu
end function elastic_mu
module function elastic_nu(ph) result(nu)
real(pReal) :: nu
integer, intent(in) :: ph
nu = param(ph)%mu
end function elastic_nu
end submodule elastic end submodule elastic

47
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -36,8 +36,8 @@ submodule(phase:plastic) dislotungsten
forestProjection forestProjection
real(pReal), allocatable, dimension(:,:,:) :: & real(pReal), allocatable, dimension(:,:,:) :: &
P_sl, & P_sl, &
nonSchmid_pos, & P_nS_pos, &
nonSchmid_neg P_nS_neg
integer :: & integer :: &
sum_N_sl !< total number of active slip system sum_N_sl !< total number of active slip system
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
@ -129,30 +129,28 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray) prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
#endif #endif
! This data is read in already in lattice prm%mu = elastic_mu(ph)
prm%mu = lattice_mu(ph)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),& prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(trim(phase%get_asString('lattice')) == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray) a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
else else
prm%nonSchmid_pos = prm%P_sl prm%P_nS_pos = prm%P_sl
prm%nonSchmid_neg = prm%P_sl prm%P_nS_neg = prm%P_sl
endif endif
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), & prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), &
phase%get_asString('lattice')) phase_lattice(ph))
prm%forestProjection = lattice_forestProjection_edge(N_sl,phase%get_asString('lattice'),& prm%forestProjection = lattice_forestProjection_edge(N_sl,phase_lattice(ph),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase_cOverA(ph))
prm%forestProjection = transpose(prm%forestProjection) prm%forestProjection = transpose(prm%forestProjection)
rho_mob_0 = pl%get_as1dFloat('rho_mob_0', requiredSize=size(N_sl)) rho_mob_0 = pl%get_as1dFloat('rho_mob_0', requiredSize=size(N_sl))
@ -163,10 +161,8 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(N_sl)) prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(N_sl))
prm%tau_Peierls = pl%get_as1dFloat('tau_Peierls', requiredSize=size(N_sl)) prm%tau_Peierls = pl%get_as1dFloat('tau_Peierls', requiredSize=size(N_sl))
prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl), & prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl))
defaultVal=[(1.0_pReal,i=1,size(N_sl))]) prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl))
prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl), &
defaultVal=[(1.0_pReal,i=1,size(N_sl))])
prm%h = pl%get_as1dFloat('h', requiredSize=size(N_sl)) prm%h = pl%get_as1dFloat('h', requiredSize=size(N_sl))
prm%w = pl%get_as1dFloat('w', requiredSize=size(N_sl)) prm%w = pl%get_as1dFloat('w', requiredSize=size(N_sl))
prm%omega = pl%get_as1dFloat('omega', requiredSize=size(N_sl)) prm%omega = pl%get_as1dFloat('omega', requiredSize=size(N_sl))
@ -298,8 +294,8 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i) Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_pos(i) * prm%P_sl(k,l,i) * prm%nonSchmid_pos(m,n,i) & + ddot_gamma_dtau_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) &
+ ddot_gamma_dtau_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i) + ddot_gamma_dtau_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i)
enddo enddo
end associate end associate
@ -323,7 +319,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
real(pReal) :: & real(pReal) :: &
VacancyDiffusion VacancyDiffusion
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos, gdot_neg,& dot_gamma_pos, dot_gamma_neg,&
tau_pos,& tau_pos,&
tau_neg, & tau_neg, &
v_cl, & v_cl, &
@ -335,10 +331,10 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
dot => dotState(ph), dst => dependentState(ph)) dot => dotState(ph), dst => dependentState(ph))
call kinetics(Mp,T,ph,en,& call kinetics(Mp,T,ph,en,&
gdot_pos,gdot_neg, & dot_gamma_pos,dot_gamma_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg) tau_pos_out = tau_pos,tau_neg_out = tau_neg)
dot%gamma_sl(:,en) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs dot%gamma_sl(:,en) = (dot_gamma_pos+dot_gamma_neg) ! ToDo: needs to be abs
VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T))
where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg
@ -380,6 +376,7 @@ module subroutine dislotungsten_dependentState(ph,en)
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
dislocationSpacing dislocationSpacing
associate(prm => param(ph), stt => state(ph),dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph),dst => dependentState(ph))
dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
@ -466,8 +463,8 @@ pure subroutine kinetics(Mp,T,ph,en, &
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
do j = 1, prm%sum_N_sl do j = 1, prm%sum_N_sl
tau_pos(j) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,j)) tau_pos(j) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,j))
tau_neg(j) = math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,j)) tau_neg(j) = math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,j))
enddo enddo

101
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -184,27 +184,23 @@ module function plastic_dislotwin_init() result(myPlasticity)
#endif #endif
! This data is read in already in lattice ! This data is read in already in lattice
prm%mu = lattice_mu(ph) prm%mu = elastic_mu(ph)
prm%nu = lattice_nu(ph) prm%nu = elastic_nu(ph)
prm%C66 = lattice_C66(1:6,1:6,ph) prm%C66 = elastic_C66(ph)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),& prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'),phase_lattice(ph))
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), & prm%forestProjection = lattice_forestProjection_edge(N_sl,phase_lattice(ph),phase_cOverA(ph))
phase%get_asString('lattice'))
prm%forestProjection = lattice_forestProjection_edge(N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%forestProjection = transpose(prm%forestProjection) prm%forestProjection = transpose(prm%forestProjection)
prm%n0_sl = lattice_slip_normal(N_sl,phase%get_asString('lattice'),& prm%n0_sl = lattice_slip_normal(N_sl,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. (N_sl(1) == 12) prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. (N_sl(1) == 12)
if(prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR if (prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR
rho_mob_0 = pl%get_as1dFloat('rho_mob_0', requiredSize=size(N_sl)) rho_mob_0 = pl%get_as1dFloat('rho_mob_0', requiredSize=size(N_sl))
rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl)) rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl))
@ -265,11 +261,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray) N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(abs(N_tw)) prm%sum_N_tw = sum(abs(N_tw))
twinActive: if (prm%sum_N_tw > 0) then twinActive: if (prm%sum_N_tw > 0) then
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),& prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'), &
prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,& phase_lattice(ph))
pl%get_as1dFloat('h_tw-tw'), &
phase%get_asString('lattice'))
prm%b_tw = pl%get_as1dFloat('b_tw', requiredSize=size(N_tw)) prm%b_tw = pl%get_as1dFloat('b_tw', requiredSize=size(N_tw))
prm%t_tw = pl%get_as1dFloat('t_tw', requiredSize=size(N_tw)) prm%t_tw = pl%get_as1dFloat('t_tw', requiredSize=size(N_tw))
@ -279,11 +273,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%L_tw = pl%get_asFloat('L_tw') prm%L_tw = pl%get_asFloat('L_tw')
prm%i_tw = pl%get_asFloat('i_tw') prm%i_tw = pl%get_asFloat('i_tw')
prm%gamma_char= lattice_characteristicShear_Twin(N_tw,phase%get_asString('lattice'),& prm%gamma_char= lattice_characteristicShear_Twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%C66_tw = lattice_C66_twin(N_tw,prm%C66,phase%get_asString('lattice'),& prm%C66_tw = lattice_C66_twin(N_tw,prm%C66,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if (.not. prm%fccTwinTransNucleation) then if (.not. prm%fccTwinTransNucleation) then
prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw') prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw')
@ -324,8 +316,8 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%x_c_tr = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%x_c_tr = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%L_tr = pl%get_asFloat('L_tr') prm%L_tr = pl%get_asFloat('L_tr')
prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_as1dFloat('h_tr-tr'), & prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_as1dFloat('h_tr-tr'),&
phase%get_asString('lattice')) phase_lattice(ph))
prm%C66_tr = lattice_C66_trans(N_tr,prm%C66,pl%get_asString('lattice_tr'), & prm%C66_tr = lattice_C66_trans(N_tr,prm%C66,pl%get_asString('lattice_tr'), &
0.0_pReal, & 0.0_pReal, &
@ -391,17 +383,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
endif endif
slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,& prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), &
pl%get_as1dFloat('h_sl-tw'), & phase_lattice(ph))
phase%get_asString('lattice')) if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' interaction_sliptwin'
endif slipAndTwinActive endif slipAndTwinActive
slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then
prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,& prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,pl%get_as1dFloat('h_sl-tr'), &
pl%get_as1dFloat('h_sl-tr'), & phase_lattice(ph))
phase%get_asString('lattice')) if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' interaction_sliptrans'
endif slipAndTransActive endif slipAndTransActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -491,8 +481,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
real(pReal) :: f_unrotated real(pReal) :: f_unrotated
associate(prm => param(ph),& associate(prm => param(ph), stt => state(ph))
stt => state(ph))
f_unrotated = 1.0_pReal & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,en)) & - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
@ -531,7 +520,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
ddot_gamma_dtau, & ddot_gamma_dtau, &
tau tau
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_sl,ddot_gamma_dtau_slip dot_gamma_sl,ddot_gamma_dtau_sl
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
dot_gamma_tw,ddot_gamma_dtau_tw dot_gamma_tw,ddot_gamma_dtau_tw
real(pReal), dimension(param(ph)%sum_N_tr) :: & real(pReal), dimension(param(ph)%sum_N_tr) :: &
@ -568,15 +557,15 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
call kinetics_slip(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_slip) call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl)
slipContribution: do i = 1, prm%sum_N_sl slipContribution: do i = 1, prm%sum_N_sl
Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i) Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
enddo slipContribution enddo slipContribution
call kinetics_twin(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
twinContibution: do i = 1, prm%sum_N_tw twinContibution: do i = 1, prm%sum_N_tw
Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i) Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -584,7 +573,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
+ ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinContibution enddo twinContibution
call kinetics_trans(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
transContibution: do i = 1, prm%sum_N_tr transContibution: do i = 1, prm%sum_N_tr
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -664,7 +653,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
- sum(stt%f_tw(1:prm%sum_N_tw,en)) & - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,en)) - sum(stt%f_tr(1:prm%sum_N_tr,en))
call kinetics_slip(Mp,T,ph,en,dot_gamma_sl) call kinetics_sl(Mp,T,ph,en,dot_gamma_sl)
dot%gamma_sl(:,en) = abs(dot_gamma_sl) dot%gamma_sl(:,en) = abs(dot_gamma_sl)
rho_dip_distance_min = prm%D_a*prm%b_sl rho_dip_distance_min = prm%D_a*prm%b_sl
@ -709,10 +698,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) &
- dot_rho_dip_climb - dot_rho_dip_climb
call kinetics_twin(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw)
dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char
call kinetics_trans(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr)
dot%f_tr(:,en) = f_unrotated*dot_gamma_tr dot%f_tr(:,en) = f_unrotated*dot_gamma_tr
end associate end associate
@ -745,9 +734,7 @@ module subroutine dislotwin_dependentState(T,ph,en)
x0 x0
associate(prm => param(ph),& associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
stt => state(ph),&
dst => dependentState(ph))
sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en)) sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en))
sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en)) sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
@ -857,8 +844,8 @@ end subroutine plastic_dislotwin_results
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end ! have the optional arguments at the end
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,T,ph,en, & pure subroutine kinetics_sl(Mp,T,ph,en, &
dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip) dot_gamma_sl,ddot_gamma_dtau_sl,tau_sl)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -871,8 +858,8 @@ pure subroutine kinetics_slip(Mp,T,ph,en, &
real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: & real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: &
dot_gamma_sl dot_gamma_sl
real(pReal), dimension(param(ph)%sum_N_sl), optional, intent(out) :: & real(pReal), dimension(param(ph)%sum_N_sl), optional, intent(out) :: &
ddot_gamma_dtau_slip, & ddot_gamma_dtau_sl, &
tau_slip tau_sl
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
ddot_gamma_dtau ddot_gamma_dtau
@ -891,9 +878,7 @@ pure subroutine kinetics_slip(Mp,T,ph,en, &
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
do i = 1, prm%sum_N_sl tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)]
tau(i) = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
enddo
tau_eff = abs(tau)-dst%tau_pass(:,en) tau_eff = abs(tau)-dst%tau_pass(:,en)
@ -921,10 +906,10 @@ pure subroutine kinetics_slip(Mp,T,ph,en, &
end associate end associate
if(present(ddot_gamma_dtau_slip)) ddot_gamma_dtau_slip = ddot_gamma_dtau if(present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau
if(present(tau_slip)) tau_slip = tau if(present(tau_sl)) tau_sl = tau
end subroutine kinetics_slip end subroutine kinetics_sl
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -934,7 +919,7 @@ end subroutine kinetics_slip
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,en,& pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
dot_gamma_tw,ddot_gamma_dtau_tw) dot_gamma_tw,ddot_gamma_dtau_tw)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
@ -993,7 +978,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,en,&
if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau
end subroutine kinetics_twin end subroutine kinetics_tw
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1003,7 +988,7 @@ end subroutine kinetics_twin
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,en,& pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
dot_gamma_tr,ddot_gamma_dtau_tr) dot_gamma_tr,ddot_gamma_dtau_tr)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
@ -1061,6 +1046,6 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,en,&
if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau
end subroutine kinetics_trans end subroutine kinetics_tr
end submodule dislotwin end submodule dislotwin

2
src/phase_mechanical_plastic_isotropic.f90
View File

@ -170,6 +170,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
integer :: & integer :: &
k, l, m, n k, l, m, n
associate(prm => param(ph), stt => state(ph)) associate(prm => param(ph), stt => state(ph))
Mp_dev = math_deviatoric33(Mp) Mp_dev = math_deviatoric33(Mp)
@ -218,6 +219,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
integer :: & integer :: &
k, l, m, n k, l, m, n
associate(prm => param(ph), stt => state(ph)) associate(prm => param(ph), stt => state(ph))
tr=math_trace33(math_spherical33(Mi)) tr=math_trace33(math_spherical33(Mi))

199
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -19,11 +19,11 @@ submodule(phase:plastic) kinehardening
xi_inf_f, & xi_inf_f, &
xi_inf_b xi_inf_b
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
interaction_slipslip !< slip resistance from slip activity h_sl_sl !< slip resistance from slip activity
real(pReal), allocatable, dimension(:,:,:) :: & real(pReal), allocatable, dimension(:,:,:) :: &
P, & P, &
nonSchmid_pos, & P_nS_pos, &
nonSchmid_neg P_nS_neg
integer :: & integer :: &
sum_N_sl sum_N_sl
logical :: & logical :: &
@ -33,13 +33,14 @@ submodule(phase:plastic) kinehardening
end type tParameters end type tParameters
type :: tKinehardeningState type :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance real(pReal), pointer, dimension(:,:) :: &
crss, & !< critical resolved stress xi, & !< resistance against plastic slip
crss_back, & !< critical resolved back stress chi, & !< back stress
sense, & !< sense of acting shear stress (-1 or +1) chi_0, & !< back stress at last switch of stress sense
chi0, & !< backstress at last switch of stress sense gamma, & !< accumulated (absolute) shear
gamma0, & !< accumulated shear at last switch of stress sense gamma_0, & !< accumulated shear at last switch of stress sense
accshear !< accumulated (absolute) shear sgn_gamma !< sense of acting shear stress (-1 or +1)
end type tKinehardeningState end type tKinehardeningState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -112,21 +113,19 @@ module function plastic_kinehardening_init() result(myPlasticity)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%P = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),& prm%P = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(trim(phase%get_asString('lattice')) == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray) a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray)
if(size(a) > 0) prm%nonSchmidActive = .true. if(size(a) > 0) prm%nonSchmidActive = .true.
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
else else
prm%nonSchmid_pos = prm%P prm%P_nS_pos = prm%P
prm%nonSchmid_neg = prm%P prm%P_nS_neg = prm%P
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(N_sl, & prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), &
pl%get_as1dFloat('h_sl-sl'), & phase_lattice(ph))
phase%get_asString('lattice'))
xi_0 = pl%get_as1dFloat('xi_0', requiredSize=size(N_sl)) xi_0 = pl%get_as1dFloat('xi_0', requiredSize=size(N_sl))
prm%xi_inf_f = pl%get_as1dFloat('xi_inf_f', requiredSize=size(N_sl)) prm%xi_inf_f = pl%get_as1dFloat('xi_inf_f', requiredSize=size(N_sl))
@ -156,18 +155,17 @@ module function plastic_kinehardening_init() result(myPlasticity)
if (any(prm%xi_inf_f <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_f' if (any(prm%xi_inf_f <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_f'
if (any(prm%xi_inf_b <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_b' if (any(prm%xi_inf_b <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_b'
!ToDo: Any sensible checks for theta?
else slipActive else slipActive
xi_0 = emptyRealArray xi_0 = emptyRealArray
allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray) allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray)
allocate(prm%interaction_SlipSlip(0,0)) allocate(prm%h_sl_sl(0,0))
endif slipActive endif slipActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseID == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeDotState = size(['xi ','chi ', 'gamma']) * prm%sum_N_sl
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeDeltaState = size(['sgn_gamma', 'chi_0 ', 'gamma_0 ']) * prm%sum_N_sl
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState) call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
@ -176,40 +174,40 @@ module function plastic_kinehardening_init() result(myPlasticity)
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%crss => plasticState(ph)%state (startIndex:endIndex,:) stt%xi => plasticState(ph)%state (startIndex:endIndex,:)
stt%crss = spread(xi_0, 2, Nmembers) stt%xi = spread(xi_0, 2, Nmembers)
dot%crss => plasticState(ph)%dotState(startIndex:endIndex,:) dot%xi => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%crss_back => plasticState(ph)%state (startIndex:endIndex,:) stt%chi => plasticState(ph)%state (startIndex:endIndex,:)
dot%crss_back => plasticState(ph)%dotState(startIndex:endIndex,:) dot%chi => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%accshear => plasticState(ph)%state (startIndex:endIndex,:) stt%gamma => plasticState(ph)%state (startIndex:endIndex,:)
dot%accshear => plasticState(ph)%dotState(startIndex:endIndex,:) dot%gamma => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
o = plasticState(ph)%offsetDeltaState o = plasticState(ph)%offsetDeltaState
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%sense => plasticState(ph)%state (startIndex :endIndex ,:) stt%sgn_gamma => plasticState(ph)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) dlt%sgn_gamma => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%chi0 => plasticState(ph)%state (startIndex :endIndex ,:) stt%chi_0 => plasticState(ph)%state (startIndex :endIndex ,:)
dlt%chi0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) dlt%chi_0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%gamma0 => plasticState(ph)%state (startIndex :endIndex ,:) stt%gamma_0 => plasticState(ph)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:) dlt%gamma_0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
end associate end associate
@ -242,22 +240,22 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos,gdot_neg, & dot_gamma_pos,dot_gamma_neg, &
dgdot_dtau_pos,dgdot_dtau_neg ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
associate(prm => param(ph)) associate(prm => param(ph))
call kinetics(Mp,ph,en,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
do i = 1, prm%sum_N_sl do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i) Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtau_pos(i) * prm%P(k,l,i) * prm%nonSchmid_pos(m,n,i) & + ddot_gamma_dtau_pos(i) * prm%P(k,l,i) * prm%P_nS_pos(m,n,i) &
+ dgdot_dtau_neg(i) * prm%P(k,l,i) * prm%nonSchmid_neg(m,n,i) + ddot_gamma_dtau_neg(i) * prm%P(k,l,i) * prm%P_nS_neg(m,n,i)
enddo enddo
end associate end associate
@ -279,28 +277,27 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en)
real(pReal) :: & real(pReal) :: &
sumGamma sumGamma
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos,gdot_neg dot_gamma_pos,dot_gamma_neg
associate(prm => param(ph), stt => state(ph),& associate(prm => param(ph), stt => state(ph),dot => dotState(ph))
dot => dotState(ph))
call kinetics(Mp,ph,en,gdot_pos,gdot_neg) call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg)
dot%accshear(:,en) = abs(gdot_pos+gdot_neg) dot%gamma(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
sumGamma = sum(stt%accshear(:,en)) sumGamma = sum(stt%gamma(:,en))
dot%crss(:,en) = matmul(prm%interaction_SlipSlip,dot%accshear(:,en)) & dot%xi(:,en) = matmul(prm%h_sl_sl,dot%gamma(:,en)) &
* ( prm%h_inf_f & * ( prm%h_inf_f &
+ (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) & + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) &
* exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) & * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) &
) )
dot%crss_back(:,en) = stt%sense(:,en)*dot%accshear(:,en) * & dot%chi(:,en) = stt%sgn_gamma(:,en)*dot%gamma(:,en) * &
( prm%h_inf_b + & ( prm%h_inf_b + &
(prm%h_0_b - prm%h_inf_b & (prm%h_0_b - prm%h_inf_b &
+ prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,en))*(stt%accshear(:,en)-stt%gamma0(:,en))& + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi_0(:,en))*(stt%gamma(:,en)-stt%gamma_0(:,en))&
) *exp(-(stt%accshear(:,en)-stt%gamma0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,en))) & ) *exp(-(stt%gamma(:,en)-stt%gamma_0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi_0(:,en))) &
) )
end associate end associate
@ -320,27 +317,25 @@ module subroutine plastic_kinehardening_deltaState(Mp,ph,en)
en en
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos,gdot_neg, & dot_gamma_pos,dot_gamma_neg, &
sense sgn_gamma
associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph)) associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph))
call kinetics(Mp,ph,en,gdot_pos,gdot_neg) call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg)
sense = merge(state(ph)%sense(:,en), & ! keep existing... sgn_gamma = merge(state(ph)%sgn_gamma(:,en), &
sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined sign(1.0_pReal,dot_gamma_pos+dot_gamma_neg), &
dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction dEq0(dot_gamma_pos+dot_gamma_neg,1e-10_pReal))
where(dNeq(sgn_gamma,stt%sgn_gamma(:,en),0.1_pReal)) ! ToDo sgn_gamma*stt%sgn_gamma(:,en)<0
!-------------------------------------------------------------------------------------------------- dlt%sgn_gamma (:,en) = sgn_gamma - stt%sgn_gamma(:,en)
! switch in sense of shear? dlt%chi_0 (:,en) = abs(stt%chi(:,en)) - stt%chi_0(:,en)
where(dNeq(sense,stt%sense(:,en),0.1_pReal)) dlt%gamma_0(:,en) = stt%gamma(:,en) - stt%gamma_0(:,en)
dlt%sense (:,en) = sense - stt%sense(:,en) ! switch sense
dlt%chi0 (:,en) = abs(stt%crss_back(:,en)) - stt%chi0(:,en) ! remember current backstress magnitude
dlt%gamma0(:,en) = stt%accshear(:,en) - stt%gamma0(:,en) ! remember current accumulated shear
else where else where
dlt%sense (:,en) = 0.0_pReal dlt%sgn_gamma (:,en) = 0.0_pReal
dlt%chi0 (:,en) = 0.0_pReal dlt%chi_0 (:,en) = 0.0_pReal
dlt%gamma0(:,en) = 0.0_pReal dlt%gamma_0(:,en) = 0.0_pReal
end where end where
end associate end associate
@ -362,22 +357,22 @@ module subroutine plastic_kinehardening_results(ph,group)
outputsLoop: do o = 1,size(prm%output) outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o))) select case(trim(prm%output(o)))
case ('xi') case ('xi')
if(prm%sum_N_sl>0) call results_writeDataset(stt%crss,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%xi,group,trim(prm%output(o)), &
'resistance against plastic slip','Pa') 'resistance against plastic slip','Pa')
case ('tau_b') case ('tau_b') !ToDo: chi
if(prm%sum_N_sl>0) call results_writeDataset(stt%crss_back,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%chi,group,trim(prm%output(o)), &
'back stress against plastic slip','Pa') 'back stress','Pa')
case ('sgn(gamma)') case ('sgn(gamma)')
if(prm%sum_N_sl>0) call results_writeDataset(stt%sense,group,trim(prm%output(o)), & ! ToDo: could be int if(prm%sum_N_sl>0) call results_writeDataset(stt%sgn_gamma,group,trim(prm%output(o)), & ! ToDo: could be int
'sense of shear','1') 'sense of shear','1')
case ('chi_0') case ('chi_0')
if(prm%sum_N_sl>0) call results_writeDataset(stt%chi0,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%chi_0,group,trim(prm%output(o)), &
'tbd','Pa') 'back stress at last switch of stress sense','Pa')
case ('gamma_0') case ('gamma_0')
if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma0,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_0,group,trim(prm%output(o)), &
'tbd','1') 'plastic shear at last switch of stress sense','1')
case ('gamma') case ('gamma')
if(prm%sum_N_sl>0) call results_writeDataset(stt%accshear,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma,group,trim(prm%output(o)), &
'plastic shear','1') 'plastic shear','1')
end select end select
enddo outputsLoop enddo outputsLoop
@ -394,7 +389,7 @@ end subroutine plastic_kinehardening_results
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics(Mp,ph,en, & pure subroutine kinetics(Mp,ph,en, &
gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -403,53 +398,55 @@ pure subroutine kinetics(Mp,ph,en, &
en en
real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: &
gdot_pos, & dot_gamma_pos, &
gdot_neg dot_gamma_neg
real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: & real(pReal), intent(out), dimension(param(ph)%sum_N_sl), optional :: &
dgdot_dtau_pos, & ddot_gamma_dtau_pos, &
dgdot_dtau_neg ddot_gamma_dtau_neg
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau_pos, & tau_pos, &
tau_neg tau_neg
integer :: i integer :: i
associate(prm => param(ph), stt => state(ph)) associate(prm => param(ph), stt => state(ph))
do i = 1, prm%sum_N_sl do i = 1, prm%sum_N_sl
tau_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,en) tau_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) - stt%chi(i,en)
tau_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,en), & tau_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)) - stt%chi(i,en), &
0.0_pReal, prm%nonSchmidActive) 0.0_pReal, prm%nonSchmidActive)
enddo enddo
where(dNeq0(tau_pos)) where(dNeq0(tau_pos))
gdot_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active dot_gamma_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_pos/stt%crss(:,en))**prm%n, tau_pos) * sign(abs(tau_pos/stt%xi(:,en))**prm%n, tau_pos)
else where else where
gdot_pos = 0.0_pReal dot_gamma_pos = 0.0_pReal
end where end where
where(dNeq0(tau_neg)) where(dNeq0(tau_neg))
gdot_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 dot_gamma_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_neg/stt%crss(:,en))**prm%n, tau_neg) * sign(abs(tau_neg/stt%xi(:,en))**prm%n, tau_neg)
else where else where
gdot_neg = 0.0_pReal dot_gamma_neg = 0.0_pReal
end where end where
if (present(dgdot_dtau_pos)) then if (present(ddot_gamma_dtau_pos)) then
where(dNeq0(gdot_pos)) where(dNeq0(dot_gamma_pos))
dgdot_dtau_pos = gdot_pos*prm%n/tau_pos ddot_gamma_dtau_pos = dot_gamma_pos*prm%n/tau_pos
else where else where
dgdot_dtau_pos = 0.0_pReal ddot_gamma_dtau_pos = 0.0_pReal
end where end where
endif endif
if (present(dgdot_dtau_neg)) then if (present(ddot_gamma_dtau_neg)) then
where(dNeq0(gdot_neg)) where(dNeq0(dot_gamma_neg))
dgdot_dtau_neg = gdot_neg*prm%n/tau_neg ddot_gamma_dtau_neg = dot_gamma_neg*prm%n/tau_neg
else where else where
dgdot_dtau_neg = 0.0_pReal ddot_gamma_dtau_neg = 0.0_pReal
end where end where
endif endif
end associate end associate
end subroutine kinetics end subroutine kinetics

167
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -108,9 +108,9 @@ submodule(phase:plastic) nonlocal
forestProjection_Edge, & !< matrix of forest projections of edge dislocations forestProjection_Edge, & !< matrix of forest projections of edge dislocations
forestProjection_Screw !< matrix of forest projections of screw dislocations forestProjection_Screw !< matrix of forest projections of screw dislocations
real(pReal), dimension(:,:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
Schmid, & !< Schmid contribution P_sl, & !< Schmid contribution
nonSchmid_pos, & P_nS_pos, &
nonSchmid_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws) P_nS_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
integer :: & integer :: &
sum_N_sl = 0 sum_N_sl = 0
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -240,41 +240,35 @@ module function plastic_nonlocal_init() result(myPlasticity)
prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
! This data is read in already in lattice prm%mu = elastic_mu(ph)
prm%mu = lattice_mu(ph) prm%nu = elastic_nu(ph)
prm%nu = lattice_nu(ph)
ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(ini%N_sl)) prm%sum_N_sl = sum(abs(ini%N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%Schmid = lattice_SchmidMatrix_slip(ini%N_sl,phase%get_asString('lattice'),& prm%P_sl = lattice_SchmidMatrix_slip(ini%N_sl,phase_lattice(ph), phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(trim(phase%get_asString('lattice')) == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray) a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray)
if(size(a) > 0) prm%nonSchmidActive = .true. if(size(a) > 0) prm%nonSchmidActive = .true.
prm%nonSchmid_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1) prm%P_nS_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1) prm%P_nS_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1)
else else
prm%nonSchmid_pos = prm%Schmid prm%P_nS_pos = prm%P_sl
prm%nonSchmid_neg = prm%Schmid prm%P_nS_neg = prm%P_sl
endif endif
prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl, & prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl,pl%get_as1dFloat('h_sl-sl'), &
pl%get_as1dFloat('h_sl-sl'), & phase_lattice(ph))
phase%get_asString('lattice'))
prm%forestProjection_edge = lattice_forestProjection_edge (ini%N_sl,phase%get_asString('lattice'),& prm%forestProjection_edge = lattice_forestProjection_edge (ini%N_sl,phase_lattice(ph),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase_cOverA(ph))
prm%forestProjection_screw = lattice_forestProjection_screw(ini%N_sl,phase%get_asString('lattice'),& prm%forestProjection_screw = lattice_forestProjection_screw(ini%N_sl,phase_lattice(ph),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase_cOverA(ph))
prm%slip_direction = lattice_slip_direction (ini%N_sl,phase%get_asString('lattice'),& prm%slip_direction = lattice_slip_direction (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) prm%slip_transverse = lattice_slip_transverse(ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%slip_transverse = lattice_slip_transverse(ini%N_sl,phase%get_asString('lattice'),& prm%slip_normal = lattice_slip_normal (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%slip_normal = lattice_slip_normal (ini%N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
! collinear systems (only for octahedral slip systems in fcc) ! collinear systems (only for octahedral slip systems in fcc)
allocate(prm%colinearSystem(prm%sum_N_sl), source = -1) allocate(prm%colinearSystem(prm%sum_N_sl), source = -1)
@ -761,11 +755,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp Mp
!< derivative of Lp with respect to Mp !< derivative of Lp with respect to Mp
integer :: & integer :: &
ns, & !< short notation for the total number of active slip systems i, j, k, l, &
i, &
j, &
k, &
l, &
t, & !< dislocation type t, & !< dislocation type
s !< index of my current slip system s !< index of my current slip system
real(pReal), dimension(param(ph)%sum_N_sl,8) :: & real(pReal), dimension(param(ph)%sum_N_sl,8) :: &
@ -779,26 +769,24 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
dv_dtauNS !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau, & !< resolved shear stress including backstress terms tau, & !< resolved shear stress including backstress terms
gdotTotal !< shear rate dot_gamma !< shear rate
associate(prm => param(ph),dst=>microstructure(ph),& associate(prm => param(ph),dst=>microstructure(ph),stt=>state(ph))
stt=>state(ph))
ns = prm%sum_N_sl
!*** shortcut to state variables !*** shortcut to state variables
rho = getRho(ph,en) rho = getRho(ph,en)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
do s = 1,ns do s = 1,prm%sum_N_sl
tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s))
tauNS(s,1) = tau(s) tauNS(s,1) = tau(s)
tauNS(s,2) = tau(s) tauNS(s,2) = tau(s)
if (tau(s) > 0.0_pReal) then if (tau(s) > 0.0_pReal) then
tauNS(s,3) = math_tensordot(Mp, +prm%nonSchmid_pos(1:3,1:3,s)) tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s))
tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_neg(1:3,1:3,s)) tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s))
else else
tauNS(s,3) = math_tensordot(Mp, +prm%nonSchmid_neg(1:3,1:3,s)) tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s))
tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_pos(1:3,1:3,s)) tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s))
endif endif
enddo enddo
tauNS = tauNS + spread(dst%tau_back(:,en),2,4) tauNS = tauNS + spread(dst%tau_back(:,en),2,4)
@ -826,22 +814,22 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
stt%v(:,en) = pack(v,.true.) stt%v(:,en) = pack(v,.true.)
!*** Bauschinger effect !*** Bauschinger effect
forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
gdotTotal = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
do s = 1,ns do s = 1,prm%sum_N_sl
Lp = Lp + gdotTotal(s) * prm%Schmid(1:3,1:3,s) Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s)
forall (i=1:3,j=1:3,k=1:3,l=1:3) & forall (i=1:3,j=1:3,k=1:3,l=1:3) &
dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
+ prm%Schmid(i,j,s) * prm%Schmid(k,l,s) & + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) &
* sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) & * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) &
+ prm%Schmid(i,j,s) & + prm%P_sl(i,j,s) &
* (+ prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
- prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s)
enddo enddo
end associate end associate
@ -861,7 +849,6 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
en en
integer :: & integer :: &
ns, & ! short notation for the total number of active slip systems
c, & ! character of dislocation c, & ! character of dislocation
t, & ! type of dislocation t, & ! type of dislocation
s ! index of my current slip system s ! index of my current slip system
@ -881,11 +868,10 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
deltaDUpper ! change in maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws
associate(prm => param(ph),dst => microstructure(ph),del => deltaState(ph)) associate(prm => param(ph),dst => microstructure(ph),del => deltaState(ph))
ns = prm%sum_N_sl
!*** shortcut to state variables !*** shortcut to state variables
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en) forall (s = 1:prm%sum_N_sl, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en)
rho = getRho(ph,en) rho = getRho(ph,en)
rhoDip = rho(:,dip) rhoDip = rho(:,dip)
@ -908,7 +894,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
!*** calculate limits for stable dipole height !*** calculate limits for stable dipole height
do s = 1,prm%sum_N_sl do s = 1,prm%sum_N_sl
tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,en) tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) +dst%tau_back(s,en)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo enddo
@ -926,16 +912,16 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
!*** dissociation by stress increase !*** dissociation by stress increase
deltaRhoDipole2SingleStress = 0.0_pReal deltaRhoDipole2SingleStress = 0.0_pReal
forall (c=1:2, s=1:ns, deltaDUpper(s,c) < 0.0_pReal .and. & forall (c=1:2, s=1:prm%sum_N_sl, deltaDUpper(s,c) < 0.0_pReal .and. &
dNeq0(dUpperOld(s,c) - prm%minDipoleHeight(s,c))) & dNeq0(dUpperOld(s,c) - prm%minDipoleHeight(s,c))) &
deltaRhoDipole2SingleStress(s,8+c) = rhoDip(s,c) * deltaDUpper(s,c) & deltaRhoDipole2SingleStress(s,8+c) = rhoDip(s,c) * deltaDUpper(s,c) &
/ (dUpperOld(s,c) - prm%minDipoleHeight(s,c)) / (dUpperOld(s,c) - prm%minDipoleHeight(s,c))
forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9) forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9)
forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ph),en) = dUpper(s,c) forall (s = 1:prm%sum_N_sl, c = 1:2) plasticState(ph)%state(iD(s,c,ph),en) = dUpper(s,c)
plasticState(ph)%deltaState(:,en) = 0.0_pReal plasticState(ph)%deltaState(:,en) = 0.0_pReal
del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*prm%sum_N_sl])
end associate end associate
@ -960,7 +946,6 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
el !< current element number el !< current element number
integer :: & integer :: &
ns, & !< short notation for the total number of active slip systems
c, & !< character of dislocation c, & !< character of dislocation
t, & !< type of dislocation t, & !< type of dislocation
s !< index of my current slip system s !< index of my current slip system
@ -978,7 +963,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
real(pReal), dimension(param(ph)%sum_N_sl,4) :: & real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
v0, & v0, &
gdot !< shear rates dot_gamma !< shear rates
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau, & !< current resolved shear stress tau, & !< current resolved shear stress
v_climb !< climb velocity of edge dipoles v_climb !< climb velocity of edge dipoles
@ -994,14 +979,10 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
return return
endif endif
associate(prm => param(ph), & associate(prm => param(ph), dst => microstructure(ph), dot => dotState(ph), stt => state(ph))
dst => microstructure(ph), &
dot => dotState(ph), &
stt => state(ph))
ns = prm%sum_N_sl
tau = 0.0_pReal tau = 0.0_pReal
gdot = 0.0_pReal dot_gamma = 0.0_pReal
rho = getRho(ph,en) rho = getRho(ph,en)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
@ -1009,14 +990,14 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
rho0 = getRho0(ph,en) rho0 = getRho0(ph,en)
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
! limits for stable dipole height ! limits for stable dipole height
do s = 1,ns do s = 1,prm%sum_N_sl
tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,en) tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + dst%tau_back(s,en)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo enddo
@ -1035,22 +1016,22 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
! dislocation multiplication ! dislocation multiplication
rhoDotMultiplication = 0.0_pReal rhoDotMultiplication = 0.0_pReal
isBCC: if (phase_lattice(ph) == 'cI') then isBCC: if (phase_lattice(ph) == 'cI') then
forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) forall (s = 1:prm%sum_N_sl, sum(abs(v(s,1:4))) > 0.0_pReal)
rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication rhoDotMultiplication(s,1:2) = sum(abs(dot_gamma(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation
rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication rhoDotMultiplication(s,3:4) = sum(abs(dot_gamma(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation
endforall endforall
else isBCC else isBCC
rhoDotMultiplication(:,1:4) = spread( & rhoDotMultiplication(:,1:4) = spread( &
(sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) & (sum(abs(dot_gamma(:,1:2)),2) * prm%f_ed_mult + sum(abs(dot_gamma(:,3:4)),2)) &
* sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26 * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26
endif isBCC endif isBCC
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) forall (s = 1:prm%sum_N_sl, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
!**************************************************************************** !****************************************************************************
@ -1059,20 +1040,20 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
! formation by glide ! formation by glide
do c = 1,2 do c = 1,2
rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile * ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile * ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile + abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c))) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile * rhoSgl(:,2*c+3) * abs(dot_gamma(:,2*c)) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile * rhoSgl(:,2*c+4) * abs(dot_gamma(:,2*c-1)) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) &
+ abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) &
@ -1085,13 +1066,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
rhoDotAthermalAnnihilation = 0.0_pReal rhoDotAthermalAnnihilation = 0.0_pReal
forall (c=1:2) & forall (c=1:2) &
rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl & rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl &
* ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single * ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1))) & ! was single hitting single
+ 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single + 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single
+ rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent + rhoDip(:,c) * (abs(dot_gamma(:,2*c-1)) + abs(dot_gamma(:,2*c)))) ! single knocks dipole constituent
! annihilated screw dipoles leave edge jogs behind on the colinear system ! annihilated screw dipoles leave edge jogs behind on the colinear system
if (phase_lattice(ph) == 'cF') & if (phase_lattice(ph) == 'cF') &
forall (s = 1:ns, prm%colinearSystem(s) > 0) & forall (s = 1:prm%sum_N_sl, prm%colinearSystem(s) > 0) &
rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) &
* 0.25_pReal * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed * 0.25_pReal * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed
@ -1101,7 +1082,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
D_SD = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) ! eq. 3.53 D_SD = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) ! eq. 3.53
v_climb = D_SD * prm%mu * prm%V_at & v_climb = D_SD * prm%mu * prm%V_at &
/ (PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1)) * kB * Temperature) ! eq. 3.54 / (PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1)) * kB * Temperature) ! eq. 3.54
forall (s = 1:ns, dUpper(s,1) > dLower(s,1)) & forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) &
rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), & rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
@ -1124,7 +1105,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
else else
dot%rho(:,en) = pack(rhoDot,.true.) dot%rho(:,en) = pack(rhoDot,.true.)
dot%gamma(:,en) = sum(gdot,2) dot%gamma(:,en) = sum(dot_gamma,2)
endif endif
end associate end associate
@ -1174,7 +1155,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
v0, & v0, &
neighbor_v0, & !< dislocation glide velocity of enighboring ip neighbor_v0, & !< dislocation glide velocity of enighboring ip
gdot !< shear rates dot_gamma !< shear rates
real(pReal), dimension(3,param(ph)%sum_N_sl,4) :: & real(pReal), dimension(3,param(ph)%sum_N_sl,4) :: &
m !< direction of dislocation motion m !< direction of dislocation motion
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
@ -1200,7 +1181,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
stt => state(ph)) stt => state(ph))
ns = prm%sum_N_sl ns = prm%sum_N_sl
gdot = 0.0_pReal dot_gamma = 0.0_pReal
rho = getRho(ph,en) rho = getRho(ph,en)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
@ -1208,7 +1189,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
@ -1218,14 +1199,14 @@ function rhoDotFlux(timestep,ph,en,ip,el)
if (plasticState(ph)%nonlocal) then if (plasticState(ph)%nonlocal) then
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ...
.and. prm%C_CFL * abs(v0) * timestep & .and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG #ifdef DEBUG
if (debugConstitutive%extensive) then if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', & print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', &
maxval(abs(v0), abs(gdot) > 0.0_pReal & maxval(abs(v0), abs(dot_gamma) > 0.0_pReal &
.and. prm%C_CFL * abs(v0) * timestep & .and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el))), & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
' at a timestep of ',timestep ' at a timestep of ',timestep

238
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -33,8 +33,8 @@ submodule(phase:plastic) phenopowerlaw
real(pReal), allocatable, dimension(:,:,:) :: & real(pReal), allocatable, dimension(:,:,:) :: &
P_sl, & P_sl, &
P_tw, & P_tw, &
nonSchmid_pos, & P_nS_pos, &
nonSchmid_neg P_nS_neg
integer :: & integer :: &
sum_N_sl, & !< total number of active slip system sum_N_sl, & !< total number of active slip system
sum_N_tw !< total number of active twin systems sum_N_tw !< total number of active twin systems
@ -46,10 +46,10 @@ submodule(phase:plastic) phenopowerlaw
type :: tPhenopowerlawState type :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: & real(pReal), pointer, dimension(:,:) :: &
xi_slip, & xi_sl, &
xi_twin, & xi_tw, &
gamma_slip, & gamma_sl, &
gamma_twin gamma_tw
end type tPhenopowerlawState end type tPhenopowerlawState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -115,21 +115,19 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),& prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(phase%get_asString('lattice') == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dFloat('a_nonSchmid',defaultVal=emptyRealArray) a = pl%get_as1dFloat('a_nonSchmid',defaultVal=emptyRealArray)
if(size(a) > 0) prm%nonSchmidActive = .true. if(size(a) > 0) prm%nonSchmidActive = .true.
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
else else
prm%nonSchmid_pos = prm%P_sl prm%P_nS_pos = prm%P_sl
prm%nonSchmid_neg = prm%P_sl prm%P_nS_neg = prm%P_sl
endif endif
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl, & prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), &
pl%get_as1dFloat('h_sl-sl'), & phase_lattice(ph))
phase%get_asString('lattice'))
xi_0_sl = pl%get_as1dFloat('xi_0_sl', requiredSize=size(N_sl)) xi_0_sl = pl%get_as1dFloat('xi_0_sl', requiredSize=size(N_sl))
prm%xi_inf_sl = pl%get_as1dFloat('xi_inf_sl', requiredSize=size(N_sl)) prm%xi_inf_sl = pl%get_as1dFloat('xi_inf_sl', requiredSize=size(N_sl))
@ -164,13 +162,11 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray) N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(abs(N_tw)) prm%sum_N_tw = sum(abs(N_tw))
twinActive: if (prm%sum_N_tw > 0) then twinActive: if (prm%sum_N_tw > 0) then
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),& prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'), &
prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,& phase_lattice(ph))
pl%get_as1dFloat('h_tw-tw'), & prm%gamma_char = lattice_characteristicShear_twin(N_tw,phase_lattice(ph),&
phase%get_asString('lattice')) phase_cOverA(ph))
prm%gamma_char = lattice_characteristicShear_twin(N_tw,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
xi_0_tw = pl%get_as1dFloat('xi_0_tw',requiredSize=size(N_tw)) xi_0_tw = pl%get_as1dFloat('xi_0_tw',requiredSize=size(N_tw))
@ -200,12 +196,10 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
! slip-twin related parameters ! slip-twin related parameters
slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
prm%h_0_tw_sl = pl%get_asFloat('h_0_tw-sl') prm%h_0_tw_sl = pl%get_asFloat('h_0_tw-sl')
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,& prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), &
pl%get_as1dFloat('h_sl-tw'), & phase_lattice(ph))
phase%get_asString('lattice')) prm%h_tw_sl = lattice_interaction_TwinBySlip(N_tw,N_sl,pl%get_as1dFloat('h_tw-sl'), &
prm%h_tw_sl = lattice_interaction_TwinBySlip(N_tw,N_sl,& phase_lattice(ph))
pl%get_as1dFloat('h_tw-sl'), &
phase%get_asString('lattice'))
else slipAndTwinActive else slipAndTwinActive
allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0 allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0
allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0 allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0
@ -235,30 +229,30 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%xi_slip => plasticState(ph)%state (startIndex:endIndex,:) stt%xi_sl => plasticState(ph)%state (startIndex:endIndex,:)
stt%xi_slip = spread(xi_0_sl, 2, Nmembers) stt%xi_sl = spread(xi_0_sl, 2, Nmembers)
dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:) dot%xi_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
stt%xi_twin => plasticState(ph)%state (startIndex:endIndex,:) stt%xi_tw => plasticState(ph)%state (startIndex:endIndex,:)
stt%xi_twin = spread(xi_0_tw, 2, Nmembers) stt%xi_tw = spread(xi_0_tw, 2, Nmembers)
dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:) dot%xi_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%gamma_slip => plasticState(ph)%state (startIndex:endIndex,:) stt%gamma_sl => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_slip => plasticState(ph)%dotState(startIndex:endIndex,:) dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
stt%gamma_twin => plasticState(ph)%state (startIndex:endIndex,:) stt%gamma_tw => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(ph)%dotState(startIndex:endIndex,:) dot%gamma_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
end associate end associate
@ -293,31 +287,31 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_slip_pos,gdot_slip_neg, & dot_gamma_sl_pos,dot_gamma_sl_neg, &
dgdot_dtauslip_pos,dgdot_dtauslip_neg ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
gdot_twin,dgdot_dtautwin dot_gamma_tw,ddot_gamma_dtau_tw
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
associate(prm => param(ph)) associate(prm => param(ph))
call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg,ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg)
slipSystems: do i = 1, prm%sum_N_sl slipSystems: do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i) Lp = Lp + (dot_gamma_sl_pos(i)+dot_gamma_sl_neg(i))*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtauslip_pos(i) * prm%P_sl(k,l,i) * prm%nonSchmid_pos(m,n,i) & + ddot_gamma_dtau_sl_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) &
+ dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i) + ddot_gamma_dtau_sl_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i)
enddo slipSystems enddo slipSystems
call kinetics_twin(Mp,ph,en,gdot_twin,dgdot_dtautwin) call kinetics_tw(Mp,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
twinSystems: do i = 1, prm%sum_N_tw twinSystems: do i = 1, prm%sum_N_tw
Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i) Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtautwin(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + ddot_gamma_dtau_tw(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinSystems enddo twinSystems
end associate end associate
@ -337,46 +331,32 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en)
en en
real(pReal) :: & real(pReal) :: &
c_SlipSlip,c_TwinSlip,c_TwinTwin, & xi_sl_sat_offset,&
xi_slip_sat_offset,& sumF
sumGamma,sumF
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
left_SlipSlip,right_SlipSlip, & dot_gamma_sl_pos,dot_gamma_sl_neg, &
gdot_slip_pos,gdot_slip_neg right_SlipSlip
associate(prm => param(ph), stt => state(ph), &
dot => dotState(ph))
sumGamma = sum(stt%gamma_slip(:,en)) associate(prm => param(ph), stt => state(ph), dot => dotState(ph))
sumF = sum(stt%gamma_twin(:,en)/prm%gamma_char)
!-------------------------------------------------------------------------------------------------- call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg)
! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices dot%gamma_sl(:,en) = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
c_SlipSlip = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en))
c_TwinSlip = prm%h_0_tw_sl * sumGamma**prm%c_3
c_TwinTwin = prm%h_0_tw_tw * sumF**prm%c_4
!-------------------------------------------------------------------------------------------------- sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
! calculate left and right vectors xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
left_SlipSlip = 1.0_pReal + prm%h_int right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, &
xi_slip_sat_offset = prm%f_sat_sl_tw*sqrt(sumF) 1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))
right_SlipSlip = sign(abs(1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl, &
1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset))
!-------------------------------------------------------------------------------------------------- dot%xi_sl(:,en) = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) &
! shear rates * matmul(prm%h_sl_sl,dot%gamma_sl(:,en)*right_SlipSlip) &
call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg) + matmul(prm%h_sl_tw,dot%gamma_tw(:,en))
dot%gamma_slip(:,en) = abs(gdot_slip_pos+gdot_slip_neg)
call kinetics_twin(Mp,ph,en,dot%gamma_twin(:,en))
!-------------------------------------------------------------------------------------------------- dot%xi_tw(:,en) = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
! hardening * matmul(prm%h_tw_sl,dot%gamma_sl(:,en)) &
dot%xi_slip(:,en) = c_SlipSlip * left_SlipSlip * & + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot%gamma_tw(:,en))
matmul(prm%h_sl_sl,dot%gamma_slip(:,en)*right_SlipSlip) &
+ matmul(prm%h_sl_tw,dot%gamma_twin(:,en))
dot%xi_twin(:,en) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,en)) &
+ c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,en))
end associate end associate
end subroutine phenopowerlaw_dotState end subroutine phenopowerlaw_dotState
@ -397,17 +377,17 @@ module subroutine plastic_phenopowerlaw_results(ph,group)
select case(trim(prm%output(o))) select case(trim(prm%output(o)))
case('xi_sl') case('xi_sl')
if(prm%sum_N_sl>0) call results_writeDataset(stt%xi_slip,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%xi_sl,group,trim(prm%output(o)), &
'resistance against plastic slip','Pa') 'resistance against plastic slip','Pa')
case('gamma_sl') case('gamma_sl')
if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_slip,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_sl,group,trim(prm%output(o)), &
'plastic shear','1') 'plastic shear','1')
case('xi_tw') case('xi_tw')
if(prm%sum_N_tw>0) call results_writeDataset(stt%xi_twin,group,trim(prm%output(o)), & if(prm%sum_N_tw>0) call results_writeDataset(stt%xi_tw,group,trim(prm%output(o)), &
'resistance against twinning','Pa') 'resistance against twinning','Pa')
case('gamma_tw') case('gamma_tw')
if(prm%sum_N_tw>0) call results_writeDataset(stt%gamma_twin,group,trim(prm%output(o)), & if(prm%sum_N_tw>0) call results_writeDataset(stt%gamma_tw,group,trim(prm%output(o)), &
'twinning shear','1') 'twinning shear','1')
end select end select
@ -424,8 +404,8 @@ end subroutine plastic_phenopowerlaw_results
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,ph,en, & pure subroutine kinetics_sl(Mp,ph,en, &
gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) dot_gamma_sl_pos,dot_gamma_sl_neg,ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -434,56 +414,57 @@ pure subroutine kinetics_slip(Mp,ph,en, &
en en
real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: &
gdot_slip_pos, & dot_gamma_sl_pos, &
gdot_slip_neg dot_gamma_sl_neg
real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: & real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: &
dgdot_dtau_slip_pos, & ddot_gamma_dtau_sl_pos, &
dgdot_dtau_slip_neg ddot_gamma_dtau_sl_neg
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau_slip_pos, & tau_sl_pos, &
tau_slip_neg tau_sl_neg
integer :: i integer :: i
associate(prm => param(ph), stt => state(ph)) associate(prm => param(ph), stt => state(ph))
do i = 1, prm%sum_N_sl do i = 1, prm%sum_N_sl
tau_slip_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) tau_sl_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i))
tau_slip_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & tau_sl_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)), &
0.0_pReal, prm%nonSchmidActive) 0.0_pReal, prm%nonSchmidActive)
enddo enddo
where(dNeq0(tau_slip_pos)) where(dNeq0(tau_sl_pos))
gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active dot_gamma_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_slip_pos/stt%xi_slip(:,en))**prm%n_sl, tau_slip_pos) * sign(abs(tau_sl_pos/stt%xi_sl(:,en))**prm%n_sl, tau_sl_pos)
else where else where
gdot_slip_pos = 0.0_pReal dot_gamma_sl_pos = 0.0_pReal
end where end where
where(dNeq0(tau_slip_neg)) where(dNeq0(tau_sl_neg))
gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2 dot_gamma_sl_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_slip_neg/stt%xi_slip(:,en))**prm%n_sl, tau_slip_neg) * sign(abs(tau_sl_neg/stt%xi_sl(:,en))**prm%n_sl, tau_sl_neg)
else where else where
gdot_slip_neg = 0.0_pReal dot_gamma_sl_neg = 0.0_pReal
end where end where
if (present(dgdot_dtau_slip_pos)) then if (present(ddot_gamma_dtau_sl_pos)) then
where(dNeq0(gdot_slip_pos)) where(dNeq0(dot_gamma_sl_pos))
dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_sl/tau_slip_pos ddot_gamma_dtau_sl_pos = dot_gamma_sl_pos*prm%n_sl/tau_sl_pos
else where else where
dgdot_dtau_slip_pos = 0.0_pReal ddot_gamma_dtau_sl_pos = 0.0_pReal
end where end where
endif endif
if (present(dgdot_dtau_slip_neg)) then if (present(ddot_gamma_dtau_sl_neg)) then
where(dNeq0(gdot_slip_neg)) where(dNeq0(dot_gamma_sl_neg))
dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_sl/tau_slip_neg ddot_gamma_dtau_sl_neg = dot_gamma_sl_neg*prm%n_sl/tau_sl_neg
else where else where
dgdot_dtau_slip_neg = 0.0_pReal ddot_gamma_dtau_sl_neg = 0.0_pReal
end where end where
endif endif
end associate end associate
end subroutine kinetics_slip end subroutine kinetics_sl
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -493,8 +474,8 @@ end subroutine kinetics_slip
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,ph,en,& pure subroutine kinetics_tw(Mp,ph,en,&
gdot_twin,dgdot_dtau_twin) dot_gamma_tw,ddot_gamma_dtau_tw)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -503,37 +484,36 @@ pure subroutine kinetics_twin(Mp,ph,en,&
en en
real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
gdot_twin dot_gamma_tw
real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: & real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: &
dgdot_dtau_twin ddot_gamma_dtau_tw
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
tau_twin tau_tw
integer :: i integer :: i
associate(prm => param(ph), stt => state(ph)) associate(prm => param(ph), stt => state(ph))
do i = 1, prm%sum_N_tw tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)]
tau_twin(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
enddo
where(tau_twin > 0.0_pReal) where(tau_tw > 0.0_pReal)
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction dot_gamma_tw = (1.0_pReal-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction
* prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,en))**prm%n_tw * prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw
else where else where
gdot_twin = 0.0_pReal dot_gamma_tw = 0.0_pReal
end where end where
if (present(dgdot_dtau_twin)) then if (present(ddot_gamma_dtau_tw)) then
where(dNeq0(gdot_twin)) where(dNeq0(dot_gamma_tw))
dgdot_dtau_twin = gdot_twin*prm%n_tw/tau_twin ddot_gamma_dtau_tw = dot_gamma_tw*prm%n_tw/tau_tw
else where else where
dgdot_dtau_twin = 0.0_pReal ddot_gamma_dtau_tw = 0.0_pReal
end where end where
endif endif
end associate end associate
end subroutine kinetics_twin end subroutine kinetics_tw
end submodule phenopowerlaw end submodule phenopowerlaw

2
src/phase_thermal.f90
View File

@ -103,7 +103,7 @@ module subroutine thermal_init(phases)
param(ph)%C_p = thermal%get_asFloat('C_p',defaultVal=0.0_pReal) ! ToDo: make mandatory? param(ph)%C_p = thermal%get_asFloat('C_p',defaultVal=0.0_pReal) ! ToDo: make mandatory?
param(ph)%K(1,1) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal) ! ToDo: make mandatory? param(ph)%K(1,1) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal) ! ToDo: make mandatory?
param(ph)%K(3,3) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal) ! ToDo: depends on symmtery param(ph)%K(3,3) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal) ! ToDo: depends on symmtery
param(ph)%K = lattice_applyLatticeSymmetry33(param(ph)%K,phase_lattice(ph)) param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph))
sources => thermal%get('source',defaultVal=emptyList) sources => thermal%get('source',defaultVal=emptyList)
thermal_Nsources(ph) = sources%length thermal_Nsources(ph) = sources%length