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

Merge branch 'polishing' into 'development' 

whitespace adjustments and consistent names

See merge request damask/DAMASK!416
This commit is contained in:
Franz Roters 2021-07-23 14:15:32 +00:00
parent 2fff191421 d216e77ca6
commit b69b857542
16 changed files with 714 additions and 768 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
src/IO.f90
View File

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

344
src/lattice.f90
View File

@ -3,7 +3,7 @@
!> @author Philip Eisenlohr, 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
!> @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
!--------------------------------------------------------------------------------------------------
module lattice
@ -365,12 +365,6 @@ module lattice
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)
! 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
module procedure slipProjection_transverse
@ -384,7 +378,8 @@ module lattice
lattice_init, &
lattice_equivalent_nu, &
lattice_equivalent_mu, &
lattice_applyLatticeSymmetry33, &
lattice_symmetrize_33, &
lattice_symmetrize_C66, &
lattice_SchmidMatrix_slip, &
lattice_SchmidMatrix_twin, &
lattice_SchmidMatrix_trans, &
@ -414,53 +409,9 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine lattice_init
integer :: Nphases, ph,i
class(tNode), pointer :: &
phases, &
phase, &
mech, &
elasticity
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
enddo
call selfTest
end subroutine lattice_init
@ -468,10 +419,10 @@ end subroutine lattice_init
!--------------------------------------------------------------------------------------------------
!> @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
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), dimension(sum(Ntwin)) :: characteristicShear
@ -513,7 +464,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
myFamilies: do f = 1,size(Ntwin,1)
mySystems: do s = 1,Ntwin(f)
a = a + 1
select case(structure)
select case(lattice)
case('cF','cI')
characteristicShear(a) = 0.5_pReal*sqrt(2.0_pReal)
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
end select
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
enddo mySystems
enddo myFamilies
@ -542,10 +493,10 @@ end function lattice_characteristicShear_Twin
!--------------------------------------------------------------------------------------------------
!> @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
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), intent(in) :: cOverA !< c/a ratio
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
integer :: i
select case(structure)
select case(lattice)
case('cF')
coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,&
structure,0.0_pReal)
lattice,0.0_pReal)
case('cI')
coordinateSystem = buildCoordinateSystem(Ntwin,BCC_NSLIPSYSTEM,BCC_SYSTEMTWIN,&
structure,0.0_pReal)
lattice,0.0_pReal)
case('hP')
coordinateSystem = buildCoordinateSystem(Ntwin,HEX_NSLIPSYSTEM,HEX_SYSTEMTWIN,&
structure,cOverA)
lattice,cOverA)
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
do i = 1, sum(Ntwin)
@ -579,11 +530,11 @@ end function lattice_C66_twin
!--------------------------------------------------------------------------------------------------
!> @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)
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,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
if (structure_target == 'hP') then
if (lattice_target == 'hP') then
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,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
@ -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(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 = applyLatticeSymmetryC66(C_target_unrotated66,'hP')
elseif (structure_target == 'cI') then
C_target_unrotated66 = lattice_symmetrize_C66(C_target_unrotated66,'hP')
elseif (lattice_target == 'cI') then
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
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
do i = 1, 6
@ -686,11 +637,11 @@ end function lattice_nonSchmidMatrix
!> @brief Slip-slip interaction matrix
!> 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
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
integer, dimension(:), allocatable :: NslipMax
@ -908,7 +859,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
],shape(BCT_INTERACTIONSLIPSLIP))
select case(structure)
select case(lattice)
case('cF')
interactionTypes = FCC_INTERACTIONSLIPSLIP
NslipMax = FCC_NSLIPSYSTEM
@ -922,7 +873,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
interactionTypes = BCT_INTERACTIONSLIPSLIP
NslipMax = BCT_NSLIPSYSTEM
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
interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes)
@ -934,11 +885,11 @@ end function lattice_interaction_SlipBySlip
!> @brief Twin-twin interaction matrix
!> 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
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
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 &
],shape(HEX_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for hex
select case(structure)
select case(lattice)
case('cF')
interactionTypes = FCC_INTERACTIONTWINTWIN
NtwinMax = FCC_NTWINSYSTEM
@ -1020,7 +971,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul
interactionTypes = HEX_INTERACTIONTWINTWIN
NtwinMax = HEX_NTWINSYSTEM
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
interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes)
@ -1032,11 +983,11 @@ end function lattice_interaction_TwinByTwin
!> @brief Trans-trans interaction matrix
!> 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
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
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 &
],shape(FCC_INTERACTIONTRANSTRANS)) !< Trans-trans interaction types for fcc
if(structure == 'cF') then
if(lattice == 'cF') then
interactionTypes = FCC_INTERACTIONTRANSTRANS
NtransMax = FCC_NTRANSSYSTEM
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
interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes)
@ -1074,12 +1025,12 @@ end function lattice_interaction_TransByTrans
!> @brief Slip-twin interaction matrix
!> 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
Ntwin !< number of active twin systems per family
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
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 &
],shape(HEX_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for hex
select case(structure)
select case(lattice)
case('cF')
interactionTypes = FCC_INTERACTIONSLIPTWIN
NslipMax = FCC_NSLIPSYSTEM
@ -1225,7 +1176,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure)
NslipMax = HEX_NSLIPSYSTEM
NtwinMax = HEX_NTWINSYSTEM
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
interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes)
@ -1237,12 +1188,12 @@ end function lattice_interaction_SlipByTwin
!> @brief Slip-trans interaction matrix
!> 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
Ntrans !< number of active trans systems per family
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
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 &
],shape(FCC_INTERACTIONSLIPTRANS)) !< Slip-trans interaction types for fcc
select case(structure)
select case(lattice)
case('cF')
interactionTypes = FCC_INTERACTIONSLIPTRANS
NslipMax = FCC_NSLIPSYSTEM
NtransMax = FCC_NTRANSSYSTEM
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
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
!> 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
Nslip !< number of active slip systems per family
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
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 &
],shape(HEX_INTERACTIONTWINSLIP)) !< Twin-slip interaction types for hex
select case(structure)
select case(lattice)
case('cF')
interactionTypes = FCC_INTERACTIONTWINSLIP
NtwinMax = FCC_NTWINSYSTEM
@ -1353,7 +1304,7 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure)
NtwinMax = HEX_NTWINSYSTEM
NslipMax = HEX_NSLIPSYSTEM
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
interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes)
@ -1365,10 +1316,10 @@ end function lattice_interaction_TwinBySlip
!> @brief Schmid matrix for slip
!> 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
character(len=*), intent(in) :: structure !< lattice structure
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA
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 :: i
select case(structure)
select case(lattice)
case('cF')
NslipMax = FCC_NSLIPSYSTEM
slipSystems = FCC_SYSTEMSLIP
@ -1392,15 +1343,15 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix)
slipSystems = BCT_SYSTEMSLIP
case default
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
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)) &
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)
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
!> 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
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), dimension(3,3,sum(Ntwin)) :: SchmidMatrix
@ -1427,7 +1378,7 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
integer, dimension(:), allocatable :: NtwinMax
integer :: i
select case(structure)
select case(lattice)
case('cF')
NtwinMax = FCC_NTWINSYSTEM
twinSystems = FCC_SYSTEMTWIN
@ -1439,15 +1390,15 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
twinSystems = HEX_SYSTEMTWIN
case default
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
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)) &
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)
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
!> 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
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), dimension(3,3,sum(Ntrans)) :: SchmidMatrix
real(pReal), dimension(3,3,sum(Ntrans)) :: devNull
real(pReal) :: a_bcc, a_fcc
if (structure_target /= 'cI' .and. structure_target /= 'hP') &
call IO_error(137,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target))
if (lattice_target /= 'cI' .and. lattice_target /= 'hP') &
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)) &
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target))
if (lattice_target == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) &
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)) &
call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target))
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(lattice_target))
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA,a_fcc,a_bcc)
@ -1490,10 +1441,10 @@ end function lattice_SchmidMatrix_trans
!> @brief Schmid matrix for cleavage
!> 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
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), 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 :: i
select case(structure)
select case(lattice)
case('cF')
NcleavageMax = FCC_NCLEAVAGESYSTEM
cleavageSystems = FCC_SYSTEMCLEAVAGE
@ -1511,15 +1462,15 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid
cleavageSystems = BCC_SYSTEMCLEAVAGE
case default
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
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)) &
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)
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)
!--------------------------------------------------------------------------------------------------
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
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), dimension(3,sum(Nslip)) :: d
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))
end function lattice_slip_direction
@ -1551,16 +1502,16 @@ end function lattice_slip_direction
!--------------------------------------------------------------------------------------------------
!> @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
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), dimension(3,sum(Nslip)) :: n
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))
end function lattice_slip_normal
@ -1569,16 +1520,16 @@ end function lattice_slip_normal
!--------------------------------------------------------------------------------------------------
!> @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
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), dimension(3,sum(Nslip)) :: t
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))
end function lattice_slip_transverse
@ -1588,17 +1539,17 @@ end function lattice_slip_transverse
!> @brief Labels for slip systems
!> 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
character(len=*), intent(in) :: structure !< lattice structure
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
character(len=:), dimension(:), allocatable :: labels
real(pReal), dimension(:,:), allocatable :: slipSystems
integer, dimension(:), allocatable :: NslipMax
select case(structure)
select case(lattice)
case('cF')
NslipMax = FCC_NSLIPSYSTEM
slipSystems = FCC_SYSTEMSLIP
@ -1612,13 +1563,13 @@ function lattice_labels_slip(Nslip,structure) result(labels)
NslipMax = BCT_NSLIPSYSTEM
slipSystems = BCT_SYSTEMSLIP
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
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)) &
call IO_error(144,ext_msg='Nslip '//trim(structure))
call IO_error(144,ext_msg='Nslip '//trim(lattice))
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), intent(in) :: T
character(len=*), intent(in) :: structure
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
T_sym = 0.0_pReal
select case(structure)
select case(lattice)
case('cF','cI')
T_sym(1,1) = 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)
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
!--------------------------------------------------------------------------------------------------
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), intent(in) :: C66
character(len=*), intent(in) :: structure
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
integer :: i,j
C66_sym = 0.0_pReal
select case(structure)
select case(lattice)
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,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
end function applyLatticeSymmetryC66
end function lattice_symmetrize_C66
!--------------------------------------------------------------------------------------------------
!> @brief Labels for twin systems
!> 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
character(len=*), intent(in) :: structure !< lattice structure
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
character(len=:), dimension(:), allocatable :: labels
real(pReal), dimension(:,:), allocatable :: twinSystems
integer, dimension(:), allocatable :: NtwinMax
select case(structure)
select case(lattice)
case('cF')
NtwinMax = FCC_NTWINSYSTEM
twinSystems = FCC_SYSTEMTWIN
@ -1723,13 +1674,13 @@ function lattice_labels_twin(Ntwin,structure) result(labels)
NtwinMax = HEX_NTWINSYSTEM
twinSystems = HEX_SYSTEMTWIN
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
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)) &
call IO_error(144,ext_msg='Ntwin '//trim(structure))
call IO_error(144,ext_msg='Ntwin '//trim(lattice))
labels = getLabels(Ntwin,NtwinMax,twinSystems)
@ -1740,18 +1691,18 @@ end function lattice_labels_twin
!> @brief Projection of the transverse direction onto the slip plane
!> @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
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), dimension(sum(Nslip),sum(Nslip)) :: projection
real(pReal), dimension(3,sum(Nslip)) :: n, t
integer :: i, j
n = lattice_slip_normal (Nslip,structure,cOverA)
t = lattice_slip_transverse(Nslip,structure,cOverA)
n = lattice_slip_normal (Nslip,lattice,cOverA)
t = lattice_slip_transverse(Nslip,lattice,cOverA)
do i=1, sum(Nslip); do j=1, sum(Nslip)
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
!> @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
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), dimension(sum(Nslip),sum(Nslip)) :: projection
real(pReal), dimension(3,sum(Nslip)) :: n, d
integer :: i, j
n = lattice_slip_normal (Nslip,structure,cOverA)
d = lattice_slip_direction(Nslip,structure,cOverA)
n = lattice_slip_normal (Nslip,lattice,cOverA)
d = lattice_slip_direction(Nslip,lattice,cOverA)
do i=1, sum(Nslip); do j=1, sum(Nslip)
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
!> @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
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), dimension(3,3,sum(Nslip)) :: coordinateSystem
real(pReal), dimension(:,:), allocatable :: slipSystems
integer, dimension(:), allocatable :: NslipMax
select case(structure)
select case(lattice)
case('cF')
NslipMax = FCC_NSLIPSYSTEM
slipSystems = FCC_SYSTEMSLIP
@ -1813,15 +1764,15 @@ function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem)
slipSystems = BCT_SYSTEMSLIP
case default
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
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)) &
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
@ -1873,15 +1824,15 @@ end function buildInteraction
!> @brief Build a local coordinate system on slip, twin, trans, cleavage systems
!> @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) :: &
active, & !< # of active systems per family
potential !< # of potential systems per family
real(pReal), dimension(:,:), intent(in) :: &
system
character(len=*), intent(in) :: &
structure !< lattice structure
character(len=2), intent(in) :: &
lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: &
cOverA
real(pReal), dimension(3,3,sum(active)) :: &
@ -1895,10 +1846,10 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA)
f, & !< index of my family
s !< index of my system in current family
if (structure == 'tI' .and. cOverA > 2.0_pReal) &
call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure))
if (structure == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) &
call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure))
if (lattice == 'tI' .and. cOverA > 2.0_pReal) &
call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(lattice))
if (lattice == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) &
call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(lattice))
a = 0
activeFamilies: do f = 1,size(active,1)
@ -1906,7 +1857,7 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA)
a = a + 1
p = sum(potential(1:f-1))+s
select case(structure)
select case(lattice)
case ('cF','cI','tI')
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))
case default
call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(structure))
call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(lattice))
end select
@ -1950,9 +1901,9 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
Q, & !< Total rotation: Q = R*B
S !< Eigendeformation tensor for phase transformation
real(pReal), intent(in) :: &
cOverA, & !< c/a for target hex structure
a_bcc, & !< lattice parameter a for target bcc structure
a_fcc !< lattice parameter a for parent fcc structure
cOverA, & !< c/a for target hex lattice
a_bcc, & !< lattice parameter a for bcc target lattice
a_fcc !< lattice parameter a for fcc parent lattice
type(rotation) :: &
R, & !< Pitsch rotation
@ -2126,9 +2077,10 @@ end function getlabels
function lattice_equivalent_nu(C,assumption) result(nu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
real(pReal) :: K, mu, nu
character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
real(pReal) :: nu
real(pReal) :: K, mu
logical :: error
real(pReal), dimension(6,6) :: S
@ -2158,8 +2110,8 @@ end function lattice_equivalent_nu
function lattice_equivalent_mu(C,assumption) result(mu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
real(pReal) :: mu
character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
real(pReal) :: mu
logical :: error
real(pReal), dimension(6,6) :: S
@ -2203,10 +2155,10 @@ subroutine selfTest
do i = 1, 10
call random_number(C)
C_cF = applyLatticeSymmetryC66(C,'cI')
C_cI = applyLatticeSymmetryC66(C,'cF')
C_hP = applyLatticeSymmetryC66(C,'hP')
C_tI = applyLatticeSymmetryC66(C,'tI')
C_cF = lattice_symmetrize_C66(C,'cI')
C_cI = lattice_symmetrize_C66(C,'cF')
C_hP = lattice_symmetrize_C66(C,'hP')
C_tI = lattice_symmetrize_C66(C,'tI')
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'
@ -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'
call random_number(T)
T_cF = lattice_applyLatticeSymmetry33(T,'cI')
T_cI = lattice_applyLatticeSymmetry33(T,'cF')
T_hP = lattice_applyLatticeSymmetry33(T,'hP')
T_tI = lattice_applyLatticeSymmetry33(T,'tI')
T_cF = lattice_symmetrize_33(T,'cI')
T_cI = lattice_symmetrize_33(T,'cF')
T_hP = lattice_symmetrize_33(T,'hP')
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_hP) .and. math_I3<1.0_pReal)) error stop 'Symmetry33/hP'
@ -2244,7 +2196,7 @@ subroutine selfTest
call random_number(C)
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 = 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,'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
math_factorial = product(math_range(n))
end function math_factorial
@ -1081,6 +1082,7 @@ integer pure function math_binomial(n,k)
integer, intent(in) :: n, k
integer :: i, k_, n_
k_ = min(k,n-k)
n_ = 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
!--------------------------------------------------------------------------------------------------
integer pure function math_multinomial(alpha)
integer pure function math_multinomial(k)
integer, intent(in), dimension(:) :: alpha
integer, intent(in), dimension(:) :: k
integer :: i
math_multinomial = 1
do i = 1, size(alpha)
math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i))
enddo
math_multinomial = product([(math_binomial(sum(k(1:i)),k(i)),i=1,size(k))])
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,3) :: m
m(1:3,1) = v1-v2
m(1:3,2) = v1-v3
m(1:3,3) = v1-v4
@ -1289,6 +1290,9 @@ subroutine selfTest
if(math_binomial(49,6) /= 13983816) &
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))
if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
error stop 'math_LeviCivita(even)'

2
src/phase_damage.f90
View File

@ -107,7 +107,7 @@ module subroutine damage_init
param(ph)%mu = source%get_asFloat('mu')
param(ph)%D(1,1) = source%get_asFloat('D_11')
if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%D(3,3) = source%get_asFloat('D_33')
param(ph)%D = lattice_applyLatticeSymmetry33(param(ph)%D,phase_lattice(ph))
param(ph)%D = lattice_symmetrize_33(param(ph)%D,phase_lattice(ph))
endif
enddo

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%g_crit = src%get_as1dFloat('g_crit', requiredSize=size(N_cl))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase_lattice(ph),phase_cOverA(ph))
! expand: family => system
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
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
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_type, &
mechanics
integer :: p,k
integer :: ph,k
phases => config_material%get('phase')
allocate(active_kinematics(kinematics_length,phases%length), source = .false. )
do p = 1, phases%length
phase => phases%get(p)
do ph = 1, phases%length
phase => phases%get(ph)
mechanics => phase%get('mechanical')
kinematics => mechanics%get('eigen',defaultVal=emptyList)
do k = 1, kinematics%length
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
@ -118,17 +118,17 @@ function kinematics_active2(kinematics_label) result(active_kinematics)
phase, &
kinematics, &
kinematics_type
integer :: p
integer :: ph
phases => config_material%get('phase')
allocate(active_kinematics(phases%length), source = .false. )
do p = 1, phases%length
phase => phases%get(p)
allocate(active_kinematics(phases%length), source = .false.)
do ph = 1, phases%length
phase => phases%get(ph)
kinematics => phase%get('damage',defaultVal=emptyList)
if(kinematics%length < 1) return
kinematics_type => kinematics%get(1)
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

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)
endif
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
end associate
endif

44
src/phase_mechanical_elastic.f90
View File

@ -2,7 +2,10 @@ submodule(phase:mechanical) elastic
type :: tParameters
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
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'))
associate(prm => param(ph))
prm%C66(1,1) = elastic%get_asFloat('C_11')
prm%C66(1,2) = elastic%get_asFloat('C_12')
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')
endif
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
enddo
@ -104,10 +116,38 @@ module function phase_homogenizedC(ph,en) result(C)
case (PLASTICITY_DISLOTWIN_ID) plasticType
C = plastic_dislotwin_homogenizedC(ph,en)
case default plasticType
C = lattice_C66(1:6,1:6,ph)
C = param(ph)%C66
end select plasticType
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

95
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -36,8 +36,8 @@ submodule(phase:plastic) dislotungsten
forestProjection
real(pReal), allocatable, dimension(:,:,:) :: &
P_sl, &
nonSchmid_pos, &
nonSchmid_neg
P_nS_pos, &
P_nS_neg
integer :: &
sum_N_sl !< total number of active slip system
character(len=pStringLen), allocatable, dimension(:) :: &
@ -129,30 +129,28 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
#endif
! This data is read in already in lattice
prm%mu = lattice_mu(ph)
prm%mu = elastic_mu(ph)
!--------------------------------------------------------------------------------------------------
! slip related parameters
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
if(trim(phase%get_asString('lattice')) == 'cI') then
if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
else
prm%nonSchmid_pos = prm%P_sl
prm%nonSchmid_neg = prm%P_sl
prm%P_nS_pos = prm%P_sl
prm%P_nS_neg = prm%P_sl
endif
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), &
phase%get_asString('lattice'))
prm%forestProjection = lattice_forestProjection_edge(N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
phase_lattice(ph))
prm%forestProjection = lattice_forestProjection_edge(N_sl,phase_lattice(ph),&
phase_cOverA(ph))
prm%forestProjection = transpose(prm%forestProjection)
rho_mob_0 = pl%get_as1dFloat('rho_mob_0', requiredSize=size(N_sl))
@ -163,39 +161,37 @@ module function plastic_dislotungsten_init() result(myPlasticity)
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%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), &
defaultVal=[(1.0_pReal,i=1,size(N_sl))])
prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl))
prm%q = pl%get_as1dFloat('q_sl', 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%omega = pl%get_as1dFloat('omega', requiredSize=size(N_sl))
prm%B = pl%get_as1dFloat('B', requiredSize=size(N_sl))
prm%D = pl%get_asFloat('D')
prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl')
prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal
prm%D_a = pl%get_asFloat('D_a') * prm%b_sl
prm%D = pl%get_asFloat('D')
prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl')
prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal
prm%D_a = pl%get_asFloat('D_a') * prm%b_sl
prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.)
! expand: family => system
rho_mob_0 = math_expand(rho_mob_0, N_sl)
rho_dip_0 = math_expand(rho_dip_0, N_sl)
prm%q = math_expand(prm%q, N_sl)
prm%p = math_expand(prm%p, N_sl)
prm%Q_s = math_expand(prm%Q_s, N_sl)
prm%b_sl = math_expand(prm%b_sl, N_sl)
prm%h = math_expand(prm%h, N_sl)
prm%w = math_expand(prm%w, N_sl)
prm%omega = math_expand(prm%omega, N_sl)
prm%tau_Peierls = math_expand(prm%tau_Peierls, N_sl)
prm%v_0 = math_expand(prm%v_0, N_sl)
prm%B = math_expand(prm%B, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%f_at = math_expand(prm%f_at, N_sl)
prm%D_a = math_expand(prm%D_a, N_sl)
rho_mob_0 = math_expand(rho_mob_0, N_sl)
rho_dip_0 = math_expand(rho_dip_0, N_sl)
prm%q = math_expand(prm%q, N_sl)
prm%p = math_expand(prm%p, N_sl)
prm%Q_s = math_expand(prm%Q_s, N_sl)
prm%b_sl = math_expand(prm%b_sl, N_sl)
prm%h = math_expand(prm%h, N_sl)
prm%w = math_expand(prm%w, N_sl)
prm%omega = math_expand(prm%omega, N_sl)
prm%tau_Peierls = math_expand(prm%tau_Peierls, N_sl)
prm%v_0 = math_expand(prm%v_0, N_sl)
prm%B = math_expand(prm%B, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%f_at = math_expand(prm%f_at, N_sl)
prm%D_a = math_expand(prm%D_a, N_sl)
! sanity checks
if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0'
@ -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)
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) &
+ ddot_gamma_dtau_pos(i) * prm%P_sl(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ ddot_gamma_dtau_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(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%P_nS_neg(m,n,i)
enddo
end associate
@ -323,7 +319,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
real(pReal) :: &
VacancyDiffusion
real(pReal), dimension(param(ph)%sum_N_sl) :: &
gdot_pos, gdot_neg,&
dot_gamma_pos, dot_gamma_neg,&
tau_pos,&
tau_neg, &
v_cl, &
@ -335,10 +331,10 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
dot => dotState(ph), dst => dependentState(ph))
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)
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))
where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg
@ -380,13 +376,14 @@ module subroutine dislotungsten_dependentState(ph,en)
real(pReal), dimension(param(ph)%sum_N_sl) :: &
dislocationSpacing
associate(prm => param(ph), stt => state(ph),dst => dependentState(ph))
dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
dst%threshold_stress(:,en) = prm%mu*prm%b_sl &
* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
dst%threshold_stress(:,en) = prm%mu*prm%b_sl &
* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
end associate
@ -466,8 +463,8 @@ pure subroutine kinetics(Mp,T,ph,en, &
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
do j = 1, prm%sum_N_sl
tau_pos(j) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,j))
tau_neg(j) = math_tensordot(Mp,prm%nonSchmid_neg(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%P_nS_neg(1:3,1:3,j))
enddo

197
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -184,27 +184,23 @@ module function plastic_dislotwin_init() result(myPlasticity)
#endif
! This data is read in already in lattice
prm%mu = lattice_mu(ph)
prm%nu = lattice_nu(ph)
prm%C66 = lattice_C66(1:6,1:6,ph)
prm%mu = elastic_mu(ph)
prm%nu = elastic_nu(ph)
prm%C66 = elastic_C66(ph)
!--------------------------------------------------------------------------------------------------
! slip related parameters
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),&
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%get_asString('lattice'))
prm%forestProjection = lattice_forestProjection_edge(N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'),phase_lattice(ph))
prm%forestProjection = lattice_forestProjection_edge(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%forestProjection = transpose(prm%forestProjection)
prm%n0_sl = lattice_slip_normal(N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%n0_sl = lattice_slip_normal(N_sl,phase_lattice(ph),phase_cOverA(ph))
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_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)
prm%sum_N_tw = sum(abs(N_tw))
twinActive: if (prm%sum_N_tw > 0) then
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,&
pl%get_as1dFloat('h_tw-tw'), &
phase%get_asString('lattice'))
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'), &
phase_lattice(ph))
prm%b_tw = pl%get_as1dFloat('b_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%i_tw = pl%get_asFloat('i_tw')
prm%gamma_char= lattice_characteristicShear_Twin(N_tw,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%gamma_char= lattice_characteristicShear_Twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%C66_tw = lattice_C66_twin(N_tw,prm%C66,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%C66_tw = lattice_C66_twin(N_tw,prm%C66,phase_lattice(ph),phase_cOverA(ph))
if (.not. prm%fccTwinTransNucleation) then
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%L_tr = pl%get_asFloat('L_tr')
prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_as1dFloat('h_tr-tr'), &
phase%get_asString('lattice'))
prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_as1dFloat('h_tr-tr'),&
phase_lattice(ph))
prm%C66_tr = lattice_C66_trans(N_tr,prm%C66,pl%get_asString('lattice_tr'), &
0.0_pReal, &
@ -391,17 +383,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
endif
slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,&
pl%get_as1dFloat('h_sl-tw'), &
phase%get_asString('lattice'))
if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' interaction_sliptwin'
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), &
phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
endif slipAndTwinActive
slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then
prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,&
pl%get_as1dFloat('h_sl-tr'), &
phase%get_asString('lattice'))
if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' interaction_sliptrans'
prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,pl%get_as1dFloat('h_sl-tr'), &
phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
endif slipAndTransActive
!--------------------------------------------------------------------------------------------------
@ -491,22 +481,21 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
real(pReal) :: f_unrotated
associate(prm => param(ph),&
stt => state(ph))
associate(prm => param(ph), stt => state(ph))
f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,en))
f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,en))
homogenizedC = f_unrotated * prm%C66
do i=1,prm%sum_N_tw
homogenizedC = homogenizedC &
+ stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i)
enddo
do i=1,prm%sum_N_tr
homogenizedC = homogenizedC &
+ stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i)
enddo
homogenizedC = f_unrotated * prm%C66
do i=1,prm%sum_N_tw
homogenizedC = homogenizedC &
+ stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i)
enddo
do i=1,prm%sum_N_tr
homogenizedC = homogenizedC &
+ stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i)
enddo
end associate
@ -531,7 +520,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
ddot_gamma_dtau, &
tau
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) :: &
dot_gamma_tw,ddot_gamma_dtau_tw
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
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
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) &
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
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
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) &
@ -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)
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
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) &
@ -609,8 +598,8 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau)
ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb &
* (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
* (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
* (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
* (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
Lp = Lp + dot_gamma_sb * P_sb
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_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)
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) &
- 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
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
end associate
@ -745,54 +734,52 @@ module subroutine dislotwin_dependentState(T,ph,en)
x0
associate(prm => param(ph),&
stt => state(ph),&
dst => dependentState(ph))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
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_tw = sum(stt%f_tw(1:prm%sum_N_tw,en))
sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * T
Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * T
!* rescaled volume fraction for topology
f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw ! this is per system ...
f_over_t_tr = sumf_tr/prm%t_tr ! but this not
!* rescaled volume fraction for topology
f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw ! this is per system ...
f_over_t_tr = sumf_tr/prm%t_tr ! but this not
! ToDo ...Physically correct, but naming could be adjusted
inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl
if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
dst%Lambda_sl(:,en) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl)
inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl
if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
dst%Lambda_sl(:,en) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl)
inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
dst%Lambda_tw(:,en) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
dst%Lambda_tw(:,en) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
!* threshold stress for dislocation motion
dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
!* threshold stress for dislocation motion
dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
!* threshold stress for growing twin/martensite
if(prm%sum_N_tw == prm%sum_N_sl) &
dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) &
+ 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct?
if(prm%sum_N_tr == prm%sum_N_sl) &
dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
+ 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct?
+ prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr)
!* threshold stress for growing twin/martensite
if(prm%sum_N_tw == prm%sum_N_sl) &
dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) &
+ 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct?
if(prm%sum_N_tr == prm%sum_N_sl) &
dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
+ 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct?
+ prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr)
dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal
dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal
dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal
dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal
x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans
dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans
dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
end associate
@ -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
! have the optional arguments at the end
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,T,ph,en, &
dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip)
pure subroutine kinetics_sl(Mp,T,ph,en, &
dot_gamma_sl,ddot_gamma_dtau_sl,tau_sl)
real(pReal), dimension(3,3), intent(in) :: &
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) :: &
dot_gamma_sl
real(pReal), dimension(param(ph)%sum_N_sl), optional, intent(out) :: &
ddot_gamma_dtau_slip, &
tau_slip
ddot_gamma_dtau_sl, &
tau_sl
real(pReal), dimension(param(ph)%sum_N_sl) :: &
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))
do i = 1, prm%sum_N_sl
tau(i) = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
enddo
tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)]
tau_eff = abs(tau)-dst%tau_pass(:,en)
@ -921,10 +906,10 @@ pure subroutine kinetics_slip(Mp,T,ph,en, &
end associate
if(present(ddot_gamma_dtau_slip)) ddot_gamma_dtau_slip = ddot_gamma_dtau
if(present(tau_slip)) tau_slip = tau
if(present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau
if(present(tau_sl)) tau_sl = tau
end subroutine kinetics_slip
end subroutine kinetics_sl
!--------------------------------------------------------------------------------------------------
@ -934,8 +919,8 @@ end subroutine kinetics_slip
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,en,&
dot_gamma_tw,ddot_gamma_dtau_tw)
pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
dot_gamma_tw,ddot_gamma_dtau_tw)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -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
end subroutine kinetics_twin
end subroutine kinetics_tw
!--------------------------------------------------------------------------------------------------
@ -1003,8 +988,8 @@ end subroutine kinetics_twin
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,en,&
dot_gamma_tr,ddot_gamma_dtau_tr)
pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
dot_gamma_tr,ddot_gamma_dtau_tr)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -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
end subroutine kinetics_trans
end subroutine kinetics_tr
end submodule dislotwin

56
src/phase_mechanical_plastic_isotropic.f90
View File

@ -170,27 +170,28 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
integer :: &
k, l, m, n
associate(prm => param(ph), stt => state(ph))
Mp_dev = math_deviatoric33(Mp)
squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev)
norm_Mp_dev = sqrt(squarenorm_Mp_dev)
Mp_dev = math_deviatoric33(Mp)
squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev)
norm_Mp_dev = sqrt(squarenorm_Mp_dev)
if (norm_Mp_dev > 0.0_pReal) then
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n
if (norm_Mp_dev > 0.0_pReal) then
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n
Lp = dot_gamma * Mp_dev/norm_Mp_dev
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev
forall (k=1:3,l=1:3) &
dLp_dMp(k,l,k,l) = dLp_dMp(k,l,k,l) + 1.0_pReal
forall (k=1:3,m=1:3) &
dLp_dMp(k,k,m,m) = dLp_dMp(k,k,m,m) - 1.0_pReal/3.0_pReal
dLp_dMp = dot_gamma * dLp_dMp / norm_Mp_dev
else
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
end if
Lp = dot_gamma * Mp_dev/norm_Mp_dev
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev
forall (k=1:3,l=1:3) &
dLp_dMp(k,l,k,l) = dLp_dMp(k,l,k,l) + 1.0_pReal
forall (k=1:3,m=1:3) &
dLp_dMp(k,k,m,m) = dLp_dMp(k,k,m,m) - 1.0_pReal/3.0_pReal
dLp_dMp = dot_gamma * dLp_dMp / norm_Mp_dev
else
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
end if
end associate
@ -218,19 +219,20 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
integer :: &
k, l, m, n
associate(prm => param(ph), stt => state(ph))
tr=math_trace33(math_spherical33(Mi))
tr=math_trace33(math_spherical33(Mi))
if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero
Li = math_I3 &
* prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(en))**(-prm%n) &
* tr * abs(tr)**(prm%n-1.0_pReal)
forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n)
else
Li = 0.0_pReal
dLi_dMi = 0.0_pReal
endif
if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero
Li = math_I3 &
* prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(en))**(-prm%n) &
* tr * abs(tr)**(prm%n-1.0_pReal)
forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n)
else
Li = 0.0_pReal
dLi_dMi = 0.0_pReal
endif
end associate

245
src/phase_mechanical_plastic_kinehardening.f90
View File

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

169
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_Screw !< matrix of forest projections of screw dislocations
real(pReal), dimension(:,:,:), allocatable :: &
Schmid, & !< Schmid contribution
nonSchmid_pos, &
nonSchmid_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
P_sl, & !< Schmid contribution
P_nS_pos, &
P_nS_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
integer :: &
sum_N_sl = 0
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)
! This data is read in already in lattice
prm%mu = lattice_mu(ph)
prm%nu = lattice_nu(ph)
prm%mu = elastic_mu(ph)
prm%nu = elastic_nu(ph)
ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(ini%N_sl))
slipActive: if (prm%sum_N_sl > 0) then
prm%Schmid = lattice_SchmidMatrix_slip(ini%N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%P_sl = lattice_SchmidMatrix_slip(ini%N_sl,phase_lattice(ph), phase_cOverA(ph))
if(trim(phase%get_asString('lattice')) == 'cI') then
if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray)
if(size(a) > 0) prm%nonSchmidActive = .true.
prm%nonSchmid_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1)
prm%P_nS_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1)
prm%P_nS_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1)
else
prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid
prm%P_nS_pos = prm%P_sl
prm%P_nS_neg = prm%P_sl
endif
prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl, &
pl%get_as1dFloat('h_sl-sl'), &
phase%get_asString('lattice'))
prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl,pl%get_as1dFloat('h_sl-sl'), &
phase_lattice(ph))
prm%forestProjection_edge = lattice_forestProjection_edge (ini%N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%forestProjection_screw = lattice_forestProjection_screw(ini%N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%forestProjection_edge = lattice_forestProjection_edge (ini%N_sl,phase_lattice(ph),&
phase_cOverA(ph))
prm%forestProjection_screw = lattice_forestProjection_screw(ini%N_sl,phase_lattice(ph),&
phase_cOverA(ph))
prm%slip_direction = lattice_slip_direction (ini%N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%slip_transverse = lattice_slip_transverse(ini%N_sl,phase%get_asString('lattice'),&
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))
prm%slip_direction = lattice_slip_direction (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%slip_transverse = lattice_slip_transverse(ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%slip_normal = lattice_slip_normal (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
! collinear systems (only for octahedral slip systems in fcc)
allocate(prm%colinearSystem(prm%sum_N_sl), source = -1)
@ -761,11 +755,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp
!< derivative of Lp with respect to Mp
integer :: &
ns, & !< short notation for the total number of active slip systems
i, &
j, &
k, &
l, &
i, j, k, l, &
t, & !< dislocation type
s !< index of my current slip system
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
real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau, & !< resolved shear stress including backstress terms
gdotTotal !< shear rate
dot_gamma !< shear rate
associate(prm => param(ph),dst=>microstructure(ph),&
stt=>state(ph))
ns = prm%sum_N_sl
associate(prm => param(ph),dst=>microstructure(ph),stt=>state(ph))
!*** shortcut to state variables
rho = getRho(ph,en)
rhoSgl = rho(:,sgl)
do s = 1,ns
tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s))
do s = 1,prm%sum_N_sl
tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s))
tauNS(s,1) = tau(s)
tauNS(s,2) = tau(s)
if (tau(s) > 0.0_pReal) then
tauNS(s,3) = math_tensordot(Mp, +prm%nonSchmid_pos(1:3,1:3,s))
tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_neg(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%P_nS_neg(1:3,1:3,s))
else
tauNS(s,3) = math_tensordot(Mp, +prm%nonSchmid_neg(1:3,1:3,s))
tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_pos(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%P_nS_pos(1:3,1:3,s))
endif
enddo
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.)
!*** 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))
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
dLp_dMp = 0.0_pReal
do s = 1,ns
Lp = Lp + gdotTotal(s) * prm%Schmid(1:3,1:3,s)
do s = 1,prm%sum_N_sl
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) &
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) &
+ prm%Schmid(i,j,s) &
* (+ prm%nonSchmid_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_sl(i,j,s) &
* (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
- prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s)
enddo
end associate
@ -861,7 +849,6 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
en
integer :: &
ns, & ! short notation for the total number of active slip systems
c, & ! character of dislocation
t, & ! type of dislocation
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
associate(prm => param(ph),dst => microstructure(ph),del => deltaState(ph))
ns = prm%sum_N_sl
!*** 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:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,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:prm%sum_N_sl, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en)
rho = getRho(ph,en)
rhoDip = rho(:,dip)
@ -908,7 +894,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
!*** calculate limits for stable dipole height
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
enddo
@ -926,16 +912,16 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
!*** dissociation by stress increase
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))) &
deltaRhoDipole2SingleStress(s,8+c) = rhoDip(s,c) * deltaDUpper(s,c) &
/ (dUpperOld(s,c) - prm%minDipoleHeight(s,c))
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
del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns])
del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*prm%sum_N_sl])
end associate
@ -946,7 +932,7 @@ end subroutine plastic_nonlocal_deltaState
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
ph,en,ip,el)
ph,en,ip,el)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
@ -960,7 +946,6 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
el !< current element number
integer :: &
ns, & !< short notation for the total number of active slip systems
c, & !< character of dislocation
t, & !< type of dislocation
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) :: &
v, & !< current dislocation glide velocity
v0, &
gdot !< shear rates
dot_gamma !< shear rates
real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau, & !< current resolved shear stress
v_climb !< climb velocity of edge dipoles
@ -994,14 +979,10 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
return
endif
associate(prm => param(ph), &
dst => microstructure(ph), &
dot => dotState(ph), &
stt => state(ph))
ns = prm%sum_N_sl
associate(prm => param(ph), dst => microstructure(ph), dot => dotState(ph), stt => state(ph))
tau = 0.0_pReal
gdot = 0.0_pReal
dot_gamma = 0.0_pReal
rho = getRho(ph,en)
rhoSgl = rho(:,sgl)
@ -1009,14 +990,14 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
rho0 = getRho0(ph,en)
my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, 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)
forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
! limits for stable dipole height
do s = 1,ns
tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,en)
do s = 1,prm%sum_N_sl
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
enddo
@ -1035,22 +1016,22 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
! dislocation multiplication
rhoDotMultiplication = 0.0_pReal
isBCC: if (phase_lattice(ph) == 'cI') then
forall (s = 1:ns, 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
forall (s = 1:prm%sum_N_sl, sum(abs(v(s,1:4))) > 0.0_pReal)
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
! * 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
! * 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
else isBCC
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
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
do c = 1,2
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) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile
* ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile
+ 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 &
* ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile
* ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile
+ 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 &
* 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 &
* 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)) &
+ abs(rhoDotSingle2DipoleGlide(:,2*c+4)) &
@ -1085,13 +1066,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
rhoDotAthermalAnnihilation = 0.0_pReal
forall (c=1:2) &
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 * (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
+ rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent
* ( 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(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(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
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) &
* 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
v_climb = D_SD * prm%mu * prm%V_at &
/ (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)), &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- 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)
else
dot%rho(:,en) = pack(rhoDot,.true.)
dot%gamma(:,en) = sum(gdot,2)
dot%gamma(:,en) = sum(dot_gamma,2)
endif
end associate
@ -1174,7 +1155,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
v, & !< current dislocation glide velocity
v0, &
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) :: &
m !< direction of dislocation motion
real(pReal), dimension(3,3) :: &
@ -1200,7 +1181,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
stt => state(ph))
ns = prm%sum_N_sl
gdot = 0.0_pReal
dot_gamma = 0.0_pReal
rho = getRho(ph,en)
rhoSgl = rho(:,sgl)
@ -1208,7 +1189,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
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
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)
@ -1218,14 +1199,14 @@ function rhoDotFlux(timestep,ph,en,ip,el)
if (plasticState(ph)%nonlocal) then
!*** 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 &
> 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
if (debugConstitutive%extensive) then
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 ', &
maxval(abs(v0), abs(gdot) > 0.0_pReal &
maxval(abs(v0), abs(dot_gamma) > 0.0_pReal &
.and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
' at a timestep of ',timestep

266
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

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