have isostrain/stress replace voigt/reuss

This commit is contained in:
Philip Eisenlohr 2022-11-29 14:56:13 -05:00 committed by Sharan Roongta
parent 599e4472e8
commit 3a3be7cdc8
7 changed files with 48 additions and 50 deletions

View File

@ -655,7 +655,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
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_isotropic_mu(C,'voigt',phase_lattice_structure(co,ce)) equivalentMu = lattice_isotropic_mu(C,phase_lattice_structure(co,ce),'isostrain')
end function equivalentMu end function equivalentMu

View File

@ -2149,11 +2149,11 @@ 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_isotropic_nu(C,assumption,lattice) result(nu) pure function lattice_isotropic_nu(C,lattice,assumption) result(nu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
character(len=2), intent(in) :: lattice character(len=2), intent(in) :: lattice
character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss')
real(pReal) :: nu real(pReal) :: nu
real(pReal) :: K, mu real(pReal) :: K, mu
@ -2161,19 +2161,17 @@ pure function lattice_isotropic_nu(C,assumption,lattice) result(nu)
real(pReal), dimension(6,6) :: S real(pReal), dimension(6,6) :: S
if (IO_lc(assumption) == 'voigt') then if (IO_lc(assumption) == 'isostrain') then
K = (C(1,1)+C(2,2)+C(3,3) +2.0_pReal*(C(1,2)+C(2,3)+C(1,3))) & K = sum(C(1:3,1:3)) / 9.0_pReal
/ 9.0_pReal elseif (IO_lc(assumption) == 'isostress') then
elseif (IO_lc(assumption) == 'reuss') 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_isotropic_mu(C,assumption,lattice) mu = lattice_isotropic_mu(C,lattice,assumption)
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_isotropic_nu end function lattice_isotropic_nu
@ -2184,18 +2182,18 @@ end function lattice_isotropic_nu
!> @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 !> @details Nonlinear Mechanics of Crystals 10.1007/978-94-007-0350-6, pp 563
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_isotropic_mu(C,assumption,lattice) result(mu) pure function lattice_isotropic_mu(C,lattice,assumption) result(mu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
character(len=2), intent(in) :: lattice character(len=2), intent(in) :: lattice
character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss')
real(pReal) :: mu real(pReal) :: mu
logical :: error logical :: error
real(pReal), dimension(6,6) :: S real(pReal), dimension(6,6) :: S
if (IO_lc(assumption) == 'voigt') then if (IO_lc(assumption) == 'isostrain') then
select case(lattice) select case(lattice)
case('cF','cI') case('cF','cI')
mu = ( C(1,1) - C(1,2) + C(4,4)*3.0_pReal) / 5.0_pReal mu = ( C(1,1) - C(1,2) + C(4,4)*3.0_pReal) / 5.0_pReal
@ -2206,16 +2204,16 @@ pure function lattice_isotropic_mu(C,assumption,lattice) result(mu)
) / 15.0_pReal ) / 15.0_pReal
end select end select
elseif (IO_lc(assumption) == 'reuss') then elseif (IO_lc(assumption) == 'isostress') then
select case(lattice) select case(lattice)
case('cF','cI') case('cF','cI')
mu = 1.0_pReal & mu = 1.0_pReal &
/ ((4.0_pReal/(5.0_pReal * (C(1,1)-C(1,2)))) + (3.0_pReal/(5.0_pReal*C(4,4)))) / (4.0_pReal/(5.0_pReal*(C(1,1)-C(1,2))) + 3.0_pReal/(5.0_pReal*C(4,4)))
case default 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 end select
else else
error stop 'invalid assumption' error stop 'invalid assumption'
@ -2293,35 +2291,35 @@ subroutine selfTest
C(6,6) = C(4,4) C(6,6) = C(4,4)
C_cI = lattice_symmetrize_C66(C,'cI') C_cI = lattice_symmetrize_C66(C,'cI')
if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'voigt','cI'),1.0e-12_pReal)) error stop 'isotropic_mu/cI/voigt' if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'cI','isostrain'),1.0e-12_pReal)) error stop 'isotropic_mu/cI/isostrain'
if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'reuss','cI'),1.0e-12_pReal)) error stop 'isotropic_mu/cI/reuss' if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'cI','isostress'),1.0e-12_pReal)) error stop 'isotropic_mu/cI/isostress'
lambda = C_cI(1,2) lambda = C_cI(1,2)
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_cI,'voigt','cI')), & if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_cI,'cI','isostrain')), &
lattice_isotropic_nu(C_cI,'voigt','cI'),1.0e-12_pReal)) error stop 'isotropic_nu/cI/voigt' lattice_isotropic_nu(C_cI,'cI','isostrain'),1.0e-12_pReal)) error stop 'isotropic_nu/cI/isostrain'
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_cI,'reuss','cI')), & if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_cI,'cI','isostress')), &
lattice_isotropic_nu(C_cI,'reuss','cI'),1.0e-12_pReal)) error stop 'isotropic_nu/cI/reuss' lattice_isotropic_nu(C_cI,'cI','isostress'),1.0e-12_pReal)) error stop 'isotropic_nu/cI/isostress'
C_hP = lattice_symmetrize_C66(C,'hP') C_hP = lattice_symmetrize_C66(C,'hP')
if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'voigt','hP'),1.0e-12_pReal)) error stop 'isotropic_mu/hP/voigt' if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'hP','isostrain'),1.0e-12_pReal)) error stop 'isotropic_mu/hP/isostrain'
if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'reuss','hP'),1.0e-12_pReal)) error stop 'isotropic_mu/hP/reuss' if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'hP','isostress'),1.0e-12_pReal)) error stop 'isotropic_mu/hP/isostress'
lambda = C_hP(1,2) lambda = C_hP(1,2)
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_hP,'voigt','hP')), & if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_hP,'hP','isostrain')), &
lattice_isotropic_nu(C_hP,'voigt','hP'),1.0e-12_pReal)) error stop 'isotropic_nu/hP/voigt' lattice_isotropic_nu(C_hP,'hP','isostrain'),1.0e-12_pReal)) error stop 'isotropic_nu/hP/isostrain'
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_hP,'reuss','hP')), & if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_hP,'hP','isostress')), &
lattice_isotropic_nu(C_hP,'reuss','hP'),1.0e-12_pReal)) error stop 'isotropic_nu/hP/reuss' lattice_isotropic_nu(C_hP,'hP','isostress'),1.0e-12_pReal)) error stop 'isotropic_nu/hP/isostress'
C_tI = lattice_symmetrize_C66(C,'tI') C_tI = lattice_symmetrize_C66(C,'tI')
if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'voigt','tI'),1.0e-12_pReal)) error stop 'isotropic_mu/tI/voigt' if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'tI','isostrain'),1.0e-12_pReal)) error stop 'isotropic_mu/tI/isostrain'
if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'reuss','tI'),1.0e-12_pReal)) error stop 'isotropic_mu/tI/reuss' if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'tI','isostress'),1.0e-12_pReal)) error stop 'isotropic_mu/tI/isostress'
lambda = C_tI(1,2) lambda = C_tI(1,2)
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_tI,'voigt','tI')), & if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_tI,'tI','isostrain')), &
lattice_isotropic_nu(C_tI,'voigt','tI'),1.0e-12_pReal)) error stop 'isotropic_nu/tI/voigt' lattice_isotropic_nu(C_tI,'tI','isostrain'),1.0e-12_pReal)) error stop 'isotropic_nu/tI/isostrain'
if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_tI,'reuss','tI')), & if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_tI,'tI','isostress')), &
lattice_isotropic_nu(C_tI,'reuss','tI'),1.0e-12_pReal)) error stop 'isotropic_nu/tI/reuss' lattice_isotropic_nu(C_tI,'tI','isostress'),1.0e-12_pReal)) error stop 'isotropic_nu/tI/isostress'
end subroutine selfTest end subroutine selfTest

View File

@ -172,13 +172,13 @@ submodule(phase) mechanical
pure module function elastic_mu(ph,en,isotropic_bound) 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=5), intent(in) :: isotropic_bound character(len=*), intent(in) :: isotropic_bound
end function elastic_mu end function elastic_mu
pure module function elastic_nu(ph,en,isotropic_bound) 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=5), intent(in) :: isotropic_bound character(len=*), intent(in) :: isotropic_bound
end function elastic_nu end function elastic_nu
end interface end interface

View File

@ -107,13 +107,13 @@ pure module function elastic_mu(ph,en,isotropic_bound) result(mu)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
character(len=5), intent(in) :: isotropic_bound character(len=*), intent(in) :: isotropic_bound
real(pReal) :: & real(pReal) :: &
mu mu
associate(prm => param(ph)) associate(prm => param(ph))
mu = lattice_isotropic_mu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph)) mu = lattice_isotropic_mu(elastic_C66(ph,en),phase_lattice(ph),isotropic_bound)
end associate end associate
@ -128,13 +128,13 @@ pure module function elastic_nu(ph,en,isotropic_bound) result(nu)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
character(len=5), intent(in) :: isotropic_bound character(len=*), intent(in) :: isotropic_bound
real(pReal) :: & real(pReal) :: &
nu nu
associate(prm => param(ph)) associate(prm => param(ph))
nu = lattice_isotropic_nu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph)) nu = lattice_isotropic_nu(elastic_C66(ph,en),phase_lattice(ph),isotropic_bound)
end associate end associate

View File

@ -35,7 +35,7 @@ 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=5) :: & character(len=9) :: &
isotropic_bound isotropic_bound
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
@ -133,7 +133,7 @@ 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='Voigt') prm%isotropic_bound = pl%get_asString('isotropic_bound',defaultVal='isostrain')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters

View File

@ -74,7 +74,7 @@ submodule(phase:plastic) dislotwin
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
character(len=5) :: & character(len=9) :: &
isotropic_bound isotropic_bound
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
@ -188,7 +188,7 @@ 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='Voigt') prm%isotropic_bound = pl%get_asString('isotropic_bound',defaultVal='isostrain')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters

View File

@ -115,7 +115,7 @@ 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=5) :: & character(len=9) :: &
isotropic_bound isotropic_bound
character(len=pStringLen), dimension(:), allocatable :: & character(len=pStringLen), dimension(:), allocatable :: &
output output
@ -243,7 +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='Voigt') 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)