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

Merge remote-tracking branch 'origin/select-mu' into development

This commit is contained in:
Martin Diehl 2022-12-20 00:41:26 +01:00
parent 902f9d2e97 dc54f5935c
commit 38e6f4b49b
7 changed files with 165 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/homogenization_mechanical_RGC.f90
View File

@ -652,9 +652,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
C = phase_homogenizedC66(material_phaseID(co,ce),material_phaseEntry(co,ce)) ! damage not included! C = phase_homogenizedC66(material_phaseID(co,ce),material_phaseEntry(co,ce)) ! damage not included!
equivalentMu = lattice_equivalent_mu(C,'voigt')
equivalentMu = lattice_isotropic_mu(C,'isostrain')
end function equivalentMu end function equivalentMu

163
src/lattice.f90
View File

@ -376,8 +376,8 @@ module lattice
public :: & public :: &
lattice_init, & lattice_init, &
lattice_equivalent_nu, & lattice_isotropic_nu, &
lattice_equivalent_mu, & lattice_isotropic_mu, &
lattice_symmetrize_33, & lattice_symmetrize_33, &
lattice_symmetrize_C66, & lattice_symmetrize_C66, &
lattice_SchmidMatrix_slip, & lattice_SchmidMatrix_slip, &
@ -422,7 +422,7 @@ end subroutine lattice_init
function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear) function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(sum(Ntwin)) :: characteristicShear real(pReal), dimension(sum(Ntwin)) :: characteristicShear
@ -496,7 +496,7 @@ end function lattice_characteristicShear_Twin
function lattice_C66_twin(Ntwin,C66,lattice,CoverA) function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), dimension(6,6), intent(in) :: C66 !< unrotated parent stiffness matrix real(pReal), dimension(6,6), intent(in) :: C66 !< unrotated parent stiffness matrix
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin
@ -535,7 +535,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
cOverA_trans,a_cF,a_cI) cOverA_trans,a_cF,a_cI)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol)
real(pReal), dimension(6,6), intent(in) :: C_parent66 real(pReal), dimension(6,6), intent(in) :: C_parent66
real(pReal), optional, intent(in) :: cOverA_trans, a_cF, a_cI real(pReal), optional, intent(in) :: cOverA_trans, a_cF, a_cI
real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans
@ -647,7 +647,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), dimension(sum(Nslip),sum(Nslip)) :: interactionMatrix real(pReal), dimension(sum(Nslip),sum(Nslip)) :: interactionMatrix
integer, dimension(:), allocatable :: NslipMax integer, dimension(:), allocatable :: NslipMax
@ -965,7 +965,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: interactionMatrix real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: interactionMatrix
integer, dimension(:), allocatable :: NtwinMax integer, dimension(:), allocatable :: NtwinMax
@ -1064,7 +1064,7 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) resu
integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction
character(len=2), intent(in) :: lattice !<Bravais lattice (Pearson symbol) (parent crystal) character(len=*), intent(in) :: lattice !<Bravais lattice (Pearson symbol) (parent crystal)
real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: interactionMatrix real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: interactionMatrix
integer, dimension(:), allocatable :: NtransMax integer, dimension(:), allocatable :: NtransMax
@ -1107,7 +1107,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntwin !< number of active twin systems per family Ntwin !< number of active twin systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-twin interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-twin interaction
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: interactionMatrix real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: interactionMatrix
integer, dimension(:), allocatable :: NslipMax, & integer, dimension(:), allocatable :: NslipMax, &
@ -1267,7 +1267,7 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice)
integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntrans !< number of active trans systems per family Ntrans !< number of active trans systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-trans interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-trans interaction
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) (parent crystal) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) (parent crystal)
real(pReal), dimension(sum(Nslip),sum(Ntrans)) :: interactionMatrix real(pReal), dimension(sum(Nslip),sum(Ntrans)) :: interactionMatrix
integer, dimension(:), allocatable :: NslipMax, & integer, dimension(:), allocatable :: NslipMax, &
@ -1320,7 +1320,7 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) r
integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family
Nslip !< number of active slip systems per family Nslip !< number of active slip systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: interactionMatrix real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: interactionMatrix
integer, dimension(:), allocatable :: NtwinMax, & integer, dimension(:), allocatable :: NtwinMax, &
@ -1394,7 +1394,7 @@ end function lattice_interaction_TwinBySlip
function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix) function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA real(pReal), intent(in) :: cOverA
real(pReal), dimension(3,3,sum(Nslip)) :: SchmidMatrix real(pReal), dimension(3,3,sum(Nslip)) :: SchmidMatrix
@ -1444,7 +1444,7 @@ end function lattice_SchmidMatrix_slip
function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix) function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,3,sum(Ntwin)) :: SchmidMatrix real(pReal), dimension(3,3,sum(Ntwin)) :: SchmidMatrix
@ -1491,7 +1491,7 @@ end function lattice_SchmidMatrix_twin
function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_cF,a_cI) result(SchmidMatrix) function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_cF,a_cI) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol)
real(pReal), optional, intent(in) :: cOverA, a_cI, a_cF real(pReal), optional, intent(in) :: cOverA, a_cI, a_cF
real(pReal), dimension(3,3,sum(Ntrans)) :: SchmidMatrix real(pReal), dimension(3,3,sum(Ntrans)) :: SchmidMatrix
@ -1520,7 +1520,7 @@ end function lattice_SchmidMatrix_trans
function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMatrix) function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,3,3,sum(Ncleavage)) :: SchmidMatrix real(pReal), dimension(3,3,3,sum(Ncleavage)) :: SchmidMatrix
@ -1563,7 +1563,7 @@ end function lattice_SchmidMatrix_cleavage
function lattice_slip_direction(Nslip,lattice,cOverA) result(d) function lattice_slip_direction(Nslip,lattice,cOverA) result(d)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,sum(Nslip)) :: d real(pReal), dimension(3,sum(Nslip)) :: d
@ -1581,7 +1581,7 @@ end function lattice_slip_direction
function lattice_slip_normal(Nslip,lattice,cOverA) result(n) function lattice_slip_normal(Nslip,lattice,cOverA) result(n)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,sum(Nslip)) :: n real(pReal), dimension(3,sum(Nslip)) :: n
@ -1599,7 +1599,7 @@ end function lattice_slip_normal
function lattice_slip_transverse(Nslip,lattice,cOverA) result(t) function lattice_slip_transverse(Nslip,lattice,cOverA) result(t)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,sum(Nslip)) :: t real(pReal), dimension(3,sum(Nslip)) :: t
@ -1618,7 +1618,7 @@ end function lattice_slip_transverse
function lattice_labels_slip(Nslip,lattice) result(labels) function lattice_labels_slip(Nslip,lattice) result(labels)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
character(len=:), dimension(:), allocatable :: labels character(len=:), dimension(:), allocatable :: labels
@ -1660,7 +1660,7 @@ pure function lattice_symmetrize_33(T,lattice) result(T_sym)
real(pReal), dimension(3,3) :: T_sym real(pReal), dimension(3,3) :: T_sym
real(pReal), dimension(3,3), intent(in) :: T real(pReal), dimension(3,3), intent(in) :: T
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
T_sym = 0.0_pReal T_sym = 0.0_pReal
@ -1688,7 +1688,7 @@ pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym)
real(pReal), dimension(6,6) :: C66_sym real(pReal), dimension(6,6) :: C66_sym
real(pReal), dimension(6,6), intent(in) :: C66 real(pReal), dimension(6,6), intent(in) :: C66
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
integer :: i,j integer :: i,j
@ -1732,7 +1732,7 @@ end function lattice_symmetrize_C66
function lattice_labels_twin(Ntwin,lattice) result(labels) function lattice_labels_twin(Ntwin,lattice) result(labels)
integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
character(len=:), dimension(:), allocatable :: labels character(len=:), dimension(:), allocatable :: labels
@ -1770,7 +1770,7 @@ end function lattice_labels_twin
function slipProjection_transverse(Nslip,lattice,cOverA) result(projection) function slipProjection_transverse(Nslip,lattice,cOverA) result(projection)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection
@ -1794,7 +1794,7 @@ end function slipProjection_transverse
function slipProjection_direction(Nslip,lattice,cOverA) result(projection) function slipProjection_direction(Nslip,lattice,cOverA) result(projection)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection
@ -1818,7 +1818,7 @@ end function slipProjection_direction
function coordinateSystem_slip(Nslip,lattice,cOverA) result(coordinateSystem) function coordinateSystem_slip(Nslip,lattice,cOverA) result(coordinateSystem)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=2), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem
@ -1907,7 +1907,7 @@ function buildCoordinateSystem(active,potential,system,lattice,cOverA)
potential !< # of potential systems per family potential !< # of potential systems per family
real(pReal), dimension(:,:), intent(in) :: & real(pReal), dimension(:,:), intent(in) :: &
system system
character(len=2), intent(in) :: & character(len=*), intent(in) :: &
lattice !< Bravais lattice (Pearson symbol) lattice !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
cOverA cOverA
@ -2149,62 +2149,85 @@ end function getlabels
!> @brief Equivalent Poisson's ratio (ν) !> @brief Equivalent Poisson's ratio (ν)
!> @details https://doi.org/10.1143/JPSJ.20.635 !> @details https://doi.org/10.1143/JPSJ.20.635
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_equivalent_nu(C,assumption) result(nu) pure function lattice_isotropic_nu(C,assumption,lattice) result(nu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss')
character(len=*), optional, intent(in) :: lattice
real(pReal) :: nu real(pReal) :: nu
real(pReal) :: K, mu real(pReal) :: K, mu
logical :: error logical :: error
real(pReal), dimension(6,6) :: S real(pReal), dimension(6,6) :: S
character(len=:), allocatable :: lattice_
if (IO_lc(assumption) == 'voigt') then lattice_ = IO_WHITESPACE
K = (C(1,1)+C(2,2)+C(3,3) +2.0_pReal*(C(1,2)+C(2,3)+C(1,3))) & if (present(lattice)) lattice_ = lattice
/ 9.0_pReal
elseif (IO_lc(assumption) == 'reuss') then if (IO_lc(assumption) == 'isostrain') then
K = sum(C(1:3,1:3)) / 9.0_pReal
elseif (IO_lc(assumption) == 'isostress') then
call math_invert(S,error,C) call math_invert(S,error,C)
if (error) error stop 'matrix inversion failed' if (error) error stop 'matrix inversion failed'
K = 1.0_pReal & K = 1.0_pReal / sum(S(1:3,1:3))
/ (S(1,1)+S(2,2)+S(3,3) +2.0_pReal*(S(1,2)+S(2,3)+S(1,3)))
else else
error stop 'invalid assumption' error stop 'invalid assumption'
end if end if
mu = lattice_equivalent_mu(C,assumption) mu = lattice_isotropic_mu(C,assumption,lattice_)
nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu) nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu)
end function lattice_equivalent_nu end function lattice_isotropic_nu
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Equivalent shear modulus (μ) !> @brief Equivalent shear modulus (μ)
!> @details https://doi.org/10.1143/JPSJ.20.635 !> @details https://doi.org/10.1143/JPSJ.20.635
!> @details Nonlinear Mechanics of Crystals 10.1007/978-94-007-0350-6, pp 563
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_equivalent_mu(C,assumption) result(mu) pure function lattice_isotropic_mu(C,assumption,lattice) result(mu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss')
character(len=*), optional, intent(in) :: lattice
real(pReal) :: mu real(pReal) :: mu
logical :: error logical :: error
real(pReal), dimension(6,6) :: S real(pReal), dimension(6,6) :: S
character(len=:), allocatable :: lattice_
if (IO_lc(assumption) == 'voigt') then lattice_ = IO_WHITESPACE
mu = (1.0_pReal*(C(1,1)+C(2,2)+C(3,3)) -1.0_pReal*(C(1,2)+C(2,3)+C(1,3)) +3.0_pReal*(C(4,4)+C(5,5)+C(6,6))) & if (present(lattice)) lattice_ = lattice
/ 15.0_pReal
elseif (IO_lc(assumption) == 'reuss') then if (IO_lc(assumption) == 'isostrain') then
select case(lattice_)
case('cF','cI')
mu = ( C(1,1) - C(1,2) + C(4,4)*3.0_pReal) / 5.0_pReal
case default
mu = ( C(1,1)+C(2,2)+C(3,3) &
- C(1,2)-C(2,3)-C(1,3) &
+(C(4,4)+C(5,5)+C(6,6)) * 3.0_pReal &
) / 15.0_pReal
end select
elseif (IO_lc(assumption) == 'isostress') then
select case(lattice_)
case('cF','cI')
mu = 5.0_pReal &
/ (4.0_pReal/(C(1,1)-C(1,2)) + 3.0_pReal/C(4,4))
case default
call math_invert(S,error,C) call math_invert(S,error,C)
if (error) error stop 'matrix inversion failed' if (error) error stop 'matrix inversion failed'
mu = 15.0_pReal & mu = 15.0_pReal &
/ (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6))) / (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)-S(1,2)-S(2,3)-S(1,3)) + 3.0_pReal*(S(4,4)+S(5,5)+S(6,6)))
end select
else else
error stop 'invalid assumption' error stop 'invalid assumption'
end if end if
end function lattice_equivalent_mu end function lattice_isotropic_mu
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -2270,16 +2293,52 @@ subroutine selfTest
call random_number(C) call random_number(C)
C(1,1) = C(1,1) + C(1,2) + 0.1_pReal C(1,1) = C(1,1) + C(1,2) + 0.1_pReal
C(1,3) = C(1,2)
C(3,3) = C(1,1)
C(4,4) = 0.5_pReal * (C(1,1) - C(1,2)) C(4,4) = 0.5_pReal * (C(1,1) - C(1,2))
C = lattice_symmetrize_C66(C,'cI') C(6,6) = C(4,4)
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'
lambda = C(1,2) C_cI = lattice_symmetrize_C66(C,'cI')
if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), & if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'isostrain','cI'),1.0e-12_pReal)) error stop 'isotropic_mu/isostrain/cI'
lattice_equivalent_nu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_nu/voigt' if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'isostress','cI'),1.0e-12_pReal)) error stop 'isotropic_mu/isostress/cI'
if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), &
lattice_equivalent_nu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_nu/reuss' lambda = C_cI(1,2)
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_cI,'isostrain','cI')), &
lattice_isotropic_nu(C_cI,'isostrain','cI'),1.0e-12_pReal)) error stop 'isotropic_nu/isostrain/cI'
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_cI,'isostress','cI')), &
lattice_isotropic_nu(C_cI,'isostress','cI'),1.0e-12_pReal)) error stop 'isotropic_nu/isostress/cI'
C_hP = lattice_symmetrize_C66(C,'hP')
if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'isostrain','hP'),1.0e-12_pReal)) error stop 'isotropic_mu/isostrain/hP'
if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'isostress','hP'),1.0e-12_pReal)) error stop 'isotropic_mu/isostress/hP'
lambda = C_hP(1,2)
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_hP,'isostrain','hP')), &
lattice_isotropic_nu(C_hP,'isostrain','hP'),1.0e-12_pReal)) error stop 'isotropic_nu/isostrain/hP'
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_hP,'isostress','hP')), &
lattice_isotropic_nu(C_hP,'isostress','hP'),1.0e-12_pReal)) error stop 'isotropic_nu/isostress/hP'
C_tI = lattice_symmetrize_C66(C,'tI')
if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'isostrain','tI'),1.0e-12_pReal)) error stop 'isotropic_mu/isostrain/tI'
if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'isostress','tI'),1.0e-12_pReal)) error stop 'isotropic_mu/isostress/tI'
lambda = C_tI(1,2)
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_tI,'isostrain','tI')), &
lattice_isotropic_nu(C_tI,'isostrain','tI'),1.0e-12_pReal)) error stop 'isotropic_nu/isostrain/tI'
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_tI,'isostress','tI')), &
lattice_isotropic_nu(C_tI,'isostress','tI'),1.0e-12_pReal)) error stop 'isotropic_nu/isostress/tI'
call random_number(C)
C = lattice_symmetrize_C66(C,'cI')
if (dNeq(lattice_isotropic_mu(C,'isostrain','cI'), lattice_isotropic_mu(C,'isostrain','hP'), 1.0e-9_pReal)) &
error stop 'isotropic_mu/isostrain/cI-hP'
if (dNeq(lattice_isotropic_nu(C,'isostrain','cF'), lattice_isotropic_nu(C,'isostrain','cI'), 1.0e-9_pReal)) &
error stop 'isotropic_nu/isostrain/cF-tI'
if (dNeq(lattice_isotropic_mu(C,'isostress','cI'), lattice_isotropic_mu(C,'isostress'), 1.0e-9_pReal)) &
error stop 'isotropic_mu/isostress/cI-hP'
if (dNeq(lattice_isotropic_nu(C,'isostress','cF'), lattice_isotropic_nu(C,'isostress'), 1.0e-9_pReal)) &
error stop 'isotropic_nu/isostress/cF-tI'
end subroutine selfTest end subroutine selfTest

6
src/phase_mechanical.f90
View File

@ -169,14 +169,16 @@ submodule(phase) mechanical
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
end function elastic_C66 end function elastic_C66
pure module function elastic_mu(ph,en) result(mu) pure module function elastic_mu(ph,en,isotropic_bound) result(mu)
real(pReal) :: mu real(pReal) :: mu
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
character(len=*), intent(in) :: isotropic_bound
end function elastic_mu end function elastic_mu
pure module function elastic_nu(ph,en) result(nu) pure module function elastic_nu(ph,en,isotropic_bound) result(nu)
real(pReal) :: nu real(pReal) :: nu
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
character(len=*), intent(in) :: isotropic_bound
end function elastic_nu end function elastic_nu
end interface end interface

18
src/phase_mechanical_elastic.f90
View File

@ -102,16 +102,21 @@ end function elastic_C66
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return shear modulus !> @brief return shear modulus
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure module function elastic_mu(ph,en) result(mu) pure module function elastic_mu(ph,en,isotropic_bound) result(mu)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
character(len=*), intent(in) :: isotropic_bound
real(pReal) :: & real(pReal) :: &
mu mu
mu = lattice_equivalent_mu(elastic_C66(ph,en),'voigt') associate(prm => param(ph))
mu = lattice_isotropic_mu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph))
end associate
end function elastic_mu end function elastic_mu
@ -119,16 +124,21 @@ end function elastic_mu
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return Poisson ratio !> @brief return Poisson ratio
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure module function elastic_nu(ph,en) result(nu) pure module function elastic_nu(ph,en,isotropic_bound) result(nu)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
character(len=*), intent(in) :: isotropic_bound
real(pReal) :: & real(pReal) :: &
nu nu
nu = lattice_equivalent_nu(elastic_C66(ph,en),'voigt') associate(prm => param(ph))
nu = lattice_isotropic_nu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph))
end associate
end function elastic_nu end function elastic_nu

8
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -35,6 +35,8 @@ submodule(phase:plastic) dislotungsten
P_nS_neg P_nS_neg
integer :: & integer :: &
sum_N_sl !< total number of active slip system sum_N_sl !< total number of active slip system
character(len=:), allocatable :: &
isotropic_bound
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
logical :: & logical :: &
@ -131,6 +133,8 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray) prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
#endif #endif
prm%isotropic_bound = pl%get_asString('isotropic_bound',defaultVal='isostrain')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
@ -333,7 +337,7 @@ module function dislotungsten_dotState(Mp,ph,en) result(dotState)
dot_rho_dip => dotState(indexDotState(ph)%rho_dip(1):indexDotState(ph)%rho_dip(2)), & dot_rho_dip => dotState(indexDotState(ph)%rho_dip(1):indexDotState(ph)%rho_dip(2)), &
dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2))) dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)))
mu = elastic_mu(ph,en) mu = elastic_mu(ph,en,prm%isotropic_bound)
T = thermal_T(ph,en) T = thermal_T(ph,en)
call kinetics(Mp,T,ph,en,& call kinetics(Mp,T,ph,en,&
@ -384,7 +388,7 @@ module subroutine dislotungsten_dependentState(ph,en)
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
dst%tau_pass(:,en) = elastic_mu(ph,en)*prm%b_sl & dst%tau_pass(:,en) = elastic_mu(ph,en,prm%isotropic_bound)*prm%b_sl &
* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
Lambda_sl_inv = 1.0_pReal/prm%D & Lambda_sl_inv = 1.0_pReal/prm%D &

19
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -73,7 +73,8 @@ submodule(phase:plastic) dislotwin
integer, allocatable, dimension(:,:) :: & integer, allocatable, dimension(:,:) :: &
fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans
character(len=:), allocatable :: & character(len=:), allocatable :: &
lattice_tr lattice_tr, &
isotropic_bound
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
logical :: & logical :: &
@ -186,6 +187,8 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray) prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
#endif #endif
prm%isotropic_bound = pl%get_asString('isotropic_bound',defaultVal='isostrain')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
@ -644,8 +647,8 @@ module function dislotwin_dotState(Mp,ph,en) result(dotState)
dot_f_tw => dotState(indexDotState(ph)%f_tw(1):indexDotState(ph)%f_tw(2)), & dot_f_tw => dotState(indexDotState(ph)%f_tw(1):indexDotState(ph)%f_tw(2)), &
dot_f_tr => dotState(indexDotState(ph)%f_tr(1):indexDotState(ph)%f_tr(2))) dot_f_tr => dotState(indexDotState(ph)%f_tr(1):indexDotState(ph)%f_tr(2)))
mu = elastic_mu(ph,en) mu = elastic_mu(ph,en,prm%isotropic_bound)
nu = elastic_nu(ph,en) nu = elastic_nu(ph,en,prm%isotropic_bound)
T = thermal_T(ph,en) T = thermal_T(ph,en)
f_matrix = 1.0_pReal & f_matrix = 1.0_pReal &
@ -732,7 +735,7 @@ module subroutine dislotwin_dependentState(ph,en)
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
mu = elastic_mu(ph,en) mu = elastic_mu(ph,en,prm%isotropic_bound)
sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en)) sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en))
sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en)) sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
@ -930,8 +933,8 @@ pure subroutine kinetics_tw(Mp,T,abs_dot_gamma_sl,ph,en,&
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
mu = elastic_mu(ph,en) mu = elastic_mu(ph,en,prm%isotropic_bound)
nu = elastic_nu(ph,en) nu = elastic_nu(ph,en,prm%isotropic_bound)
Gamma_sf = prm%Gamma_sf%at(T) Gamma_sf = prm%Gamma_sf%at(T)
tau_hat = 3.0_pReal*prm%b_tw(1)*mu/prm%L_tw & tau_hat = 3.0_pReal*prm%b_tw(1)*mu/prm%L_tw &
@ -1006,8 +1009,8 @@ pure subroutine kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,&
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
mu = elastic_mu(ph,en) mu = elastic_mu(ph,en,prm%isotropic_bound)
nu = elastic_nu(ph,en) nu = elastic_nu(ph,en,prm%isotropic_bound)
Gamma_sf = prm%Gamma_sf%at(T) Gamma_sf = prm%Gamma_sf%at(T)
tau_hat = 3.0_pReal*prm%b_tr(1)*mu/prm%L_tr & tau_hat = 3.0_pReal*prm%b_tr(1)*mu/prm%L_tr &

15
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -115,6 +115,8 @@ submodule(phase:plastic) nonlocal
sum_N_sl = 0 sum_N_sl = 0
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
colinearSystem !< colinear system to the active slip system (only valid for fcc!) colinearSystem !< colinear system to the active slip system (only valid for fcc!)
character(len=:), allocatable :: &
isotropic_bound
character(len=pStringLen), dimension(:), allocatable :: & character(len=pStringLen), dimension(:), allocatable :: &
output output
logical :: & logical :: &
@ -241,6 +243,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray) prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
#endif #endif
prm%isotropic_bound = pl%get_asString('isotropic_bound',defaultVal='isostrain')
prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
@ -609,8 +612,8 @@ module subroutine nonlocal_dependentState(ph, en)
associate(prm => param(ph),dst => dependentState(ph), stt => state(ph)) associate(prm => param(ph),dst => dependentState(ph), stt => state(ph))
mu = elastic_mu(ph,en) mu = elastic_mu(ph,en,prm%isotropic_bound)
nu = elastic_nu(ph,en) nu = elastic_nu(ph,en,prm%isotropic_bound)
rho = getRho(ph,en) rho = getRho(ph,en)
stt%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & stt%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
@ -880,8 +883,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph)) associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph))
mu = elastic_mu(ph,en) mu = elastic_mu(ph,en,prm%isotropic_bound)
nu = elastic_nu(ph,en) nu = elastic_nu(ph,en,prm%isotropic_bound)
!*** shortcut to state variables !*** shortcut to state variables
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, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
@ -994,8 +997,8 @@ module subroutine nonlocal_dotState(Mp,timestep, &
associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph)) associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph))
mu = elastic_mu(ph,en) mu = elastic_mu(ph,en,prm%isotropic_bound)
nu = elastic_nu(ph,en) nu = elastic_nu(ph,en,prm%isotropic_bound)
Temperature = thermal_T(ph,en) Temperature = thermal_T(ph,en)
tau = 0.0_pReal tau = 0.0_pReal