using central functionality

This commit is contained in:
Martin Diehl 2020-12-28 11:03:29 +01:00
parent d59cb81ca8
commit 1ac5465d65
2 changed files with 38 additions and 53 deletions

View File

@ -8,6 +8,7 @@
!--------------------------------------------------------------------------------------------------
submodule(homogenization:homogenization_mech) homogenization_mech_RGC
use rotations
use lattice
type :: tParameters
integer, dimension(:), allocatable :: &
@ -524,8 +525,10 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
real(pReal), dimension (3) :: nVect,surfCorr
real(pReal), dimension (2) :: Gmoduli
integer :: iGrain,iGNghb,iFace,i,j,k,l
real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
real(pReal), parameter :: nDefToler = 1.0e-10_pReal
real(pReal) :: muGrain,muGNghb,nDefNorm
real(pReal), parameter :: &
nDefToler = 1.0e-10_pReal, &
b = 2.5e-10_pReal ! Length of Burgers vector
nGDim = param(instance)%N_constituents
rPen = 0.0_pReal
@ -543,9 +546,7 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
!-----------------------------------------------------------------------------------------------
! computing the mismatch and penalty stress tensor of all grains
grainLoop: do iGrain = 1,product(prm%N_constituents)
Gmoduli = equivalentModuli(iGrain,ip,el)
muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain
bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector
muGrain = equivalentMu(iGrain,ip,el)
iGrain3 = grain1to3(iGrain,prm%N_constituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position
interfaceLoop: do iFace = 1,6
@ -557,9 +558,7 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
where(iGNghb3 < 1) iGNghb3 = nGDim
where(iGNghb3 >nGDim) iGNghb3 = 1
iGNghb = grain3to1(iGNghb3,prm%N_constituents) ! get the ID of the neighboring grain
Gmoduli = equivalentModuli(iGNghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor
muGNghb = Gmoduli(1)
bgGNghb = Gmoduli(2)
muGNghb = equivalentMu(iGNghb,ip,el)
gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! difference/jump in deformation gradeint across the neighbor
!-------------------------------------------------------------------------------------------
@ -579,7 +578,7 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
!-------------------------------------------------------------------------------------------
! compute the stress penalty of all interfaces
do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xi_alpha &
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*b + muGNghb*b)*prm%xi_alpha &
*surfCorr(abs(intFace(1)))/prm%D_alpha(abs(intFace(1))) &
*cosh(prm%c_alpha*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
@ -666,44 +665,26 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
end function surfaceCorrection
!--------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
!> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
!--------------------------------------------------------------------------------------------------
function equivalentModuli(grainID,ip,el)
real(pReal), dimension(2) :: equivalentModuli
!-------------------------------------------------------------------------------------------------
real(pReal) function equivalentMu(grainID,ip,el)
integer, intent(in) :: &
grainID,&
ip, & !< integration point number
el !< element number
real(pReal), dimension(6,6) :: elasTens
real(pReal) :: &
cEquiv_11, &
cEquiv_12, &
cEquiv_44
elasTens = constitutive_homogenizedC(grainID,ip,el)
!----------------------------------------------------------------------------------------------
! compute the equivalent shear modulus after Turterltaub and Suiker, JMPS (2005)
cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal
cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + &
elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal
cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal
equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44
!----------------------------------------------------------------------------------------------
! obtain the length of Burgers vector (could be model dependend)
equivalentModuli(2) = 2.5e-10_pReal
end function equivalentModuli
!--------------------------------------------------------------------------------------------------
equivalentMu = lattice_equivalent_mu(constitutive_homogenizedC(grainID,ip,el),'voigt')
end function equivalentMu
!-------------------------------------------------------------------------------------------------
!> @brief calculating the grain deformation gradient (the same with
! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
!--------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
subroutine grainDeformation(F, avgF, instance, of)
real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain
@ -718,7 +699,7 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
integer, dimension(3) :: iGrain3
integer :: iGrain,iFace,i,j
!-------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
associate(prm => param(instance))

View File

@ -421,6 +421,8 @@ module lattice
lattice_BCT_ID, &
lattice_HEX_ID, &
lattice_ORT_ID, &
lattice_equivalent_nu, &
lattice_equivalent_mu, &
lattice_applyLatticeSymmetry33, &
lattice_SchmidMatrix_slip, &
lattice_SchmidMatrix_twin, &
@ -508,8 +510,8 @@ subroutine lattice_init
lattice_C66(1:6,1:6,p) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,p),phase%get_asString('lattice'))
lattice_mu(p) = equivalent_mu(lattice_C66(1:6,1:6,p),'voigt')
lattice_nu(p) = equivalent_nu(lattice_C66(1:6,1:6,p),'voigt')
lattice_nu(p) = lattice_equivalent_nu(lattice_C66(1:6,1:6,p),'voigt')
lattice_mu(p) = lattice_equivalent_mu(lattice_C66(1:6,1:6,p),'voigt')
lattice_C66(1:6,1:6,p) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,p))) ! Literature data is in Voigt notation
do i = 1, 6
@ -2188,15 +2190,16 @@ end function getlabels
!> @brief Equivalent Poisson's ratio (ν)
!> @details https://doi.org/10.1143/JPSJ.20.635
!--------------------------------------------------------------------------------------------------
function equivalent_nu(C,assumption) result(nu)
function lattice_equivalent_nu(C,assumption) result(nu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
real(pReal) :: K, mu, nu
logical :: error
real(pReal), dimension(6,6) :: S
if (IO_lc(assumption) == 'voigt') then
K = (C(1,1)+C(2,2)+C(3,3) +2.0_pReal*(C(1,2)+C(2,3)+C(1,3))) &
/ 9.0_pReal
@ -2210,25 +2213,26 @@ function equivalent_nu(C,assumption) result(nu)
K = 0.0_pReal
endif
mu = equivalent_mu(C,assumption)
mu = lattice_equivalent_mu(C,assumption)
nu = (1.5_pReal*K -mu)/(3.0_pReal*K+mu)
end function equivalent_nu
end function lattice_equivalent_nu
!--------------------------------------------------------------------------------------------------
!> @brief Equivalent shear modulus (μ)
!> @details https://doi.org/10.1143/JPSJ.20.635
!--------------------------------------------------------------------------------------------------
function equivalent_mu(C,assumption) result(mu)
function lattice_equivalent_mu(C,assumption) result(mu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
real(pReal) :: mu
logical :: error
real(pReal), dimension(6,6) :: S
if (IO_lc(assumption) == 'voigt') then
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))) &
/ 15.0_pReal
@ -2242,7 +2246,7 @@ function equivalent_mu(C,assumption) result(mu)
mu = 0.0_pReal
endif
end function equivalent_mu
end function lattice_equivalent_mu
!--------------------------------------------------------------------------------------------------
@ -2266,14 +2270,14 @@ subroutine selfTest
call random_number(C)
C(1,1) = C(1,1) + 1.0_pReal
C = applyLatticeSymmetryC66(C,'aP')
if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
if(dNeq(C(6,6),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
if(dNeq(C(6,6),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
lambda = C(1,2)
if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'voigt')),equivalent_nu(C,'voigt'),1.0e-12_pReal)) &
error stop 'equivalent_nu/voigt'
if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'reuss')),equivalent_nu(C,'reuss'),1.0e-12_pReal)) &
error stop 'equivalent_nu/reuss'
if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), &
lattice_equivalent_nu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_nu/voigt'
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'
end subroutine selfTest