functions only used within updatestate

This commit is contained in:
Martin Diehl 2018-11-03 21:11:43 +01:00
parent c16fdec749
commit 0aa21e507a
1 changed files with 322 additions and 323 deletions

View File

@ -74,11 +74,6 @@ module homogenization_RGC
homogenization_RGC_updateState, & homogenization_RGC_updateState, &
homogenization_RGC_postResults homogenization_RGC_postResults
private :: & private :: &
stressPenalty, &
volumePenalty, &
grainDeformation, &
surfaceCorrection, &
equivalentModuli, &
relaxationVector, & relaxationVector, &
interfaceNormal, & interfaceNormal, &
getInterface, & getInterface, &
@ -798,6 +793,328 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
contains
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to deformation mismatch
!--------------------------------------------------------------------------------------------------
subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance)
use debug, only: &
debug_level, &
debug_homogenization,&
debug_levelExtensive, &
debug_e, &
debug_i
use mesh, only: &
mesh_homogenizationAt
use constitutive, only: &
constitutive_homogenizedC
use math, only: &
math_civita
use material, only: &
homogenization_maxNgrains,&
homogenization_Ngrains
use numerics, only: &
xSmoo_RGC
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen !< stress-like penalty
real(pReal), dimension (3,homogenization_maxNgrains), intent(out) :: nMis !< total amount of mismatch
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef !< deformation gradients
real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor
integer(pInt), intent(in) :: ip,el,instance
integer(pInt), dimension (4) :: intFace
integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim
real(pReal), dimension (3,3) :: gDef,nDef
real(pReal), dimension (3) :: nVect,surfCorr
real(pReal), dimension (2) :: Gmoduli
integer(pInt) :: iGrain,iGNghb,iFace,i,j,k,l
real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
type(tParameters) :: prm
integer(pInt), parameter :: nFace = 6_pInt
real(pReal), parameter :: nDefToler = 1.0e-10_pReal
nGDim = param(instance)%Nconstituents
rPen = 0.0_pReal
nMis = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! get the correction factor the modulus of penalty stress representing the evolution of area of
! the interfaces due to deformations
surfCorr = surfaceCorrection(avgF,ip,el)
associate(prm => param(instance))
!--------------------------------------------------------------------------------------------------
! debugging the surface correction factor
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el
write(6,'(1x,3(e11.4,1x))')(surfCorr(i), i = 1,3)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! computing the mismatch and penalty stress tensor of all grains
do iGrain = 1_pInt,homogenization_Ngrains(mesh_homogenizationAt(el))
Gmoduli = equivalentModuli(iGrain,ip,el)
muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain
bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector
iGrain3 = grain1to3(iGrain,instance) ! get the grain ID in local 3-dimensional index (x,y,z)-position
!* Looping over all six interfaces of each grain
do iFace = 1_pInt,nFace
intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain
nVect = interfaceNormal(intFace,ip,el) ! get the interface normal
iGNghb3 = iGrain3 ! identify the neighboring grain across the interface
iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal),pInt)
if (iGNghb3(1) < 1) iGNghb3(1) = nGDim(1) ! with periodicity along e1 direction
if (iGNghb3(1) > nGDim(1)) iGNghb3(1) = 1_pInt
if (iGNghb3(2) < 1) iGNghb3(2) = nGDim(2) ! with periodicity along e2 direction
if (iGNghb3(2) > nGDim(2)) iGNghb3(2) = 1_pInt
if (iGNghb3(3) < 1) iGNghb3(3) = nGDim(3) ! with periodicity along e3 direction
if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1_pInt
iGNghb = grain3to1(iGNghb3,instance) ! get the ID of the neighboring grain
Gmoduli = equivalentModuli(iGNghb,ip,el) ! collecting the shear modulus and Burgers vector of the neighbor
muGNghb = Gmoduli(1)
bgGNghb = Gmoduli(2)
gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! compute the difference/jump in deformation gradeint across the neighbor
!--------------------------------------------------------------------------------------------------
! compute the mismatch tensor of all interfaces
nDefNorm = 0.0_pReal
nDef = 0.0_pReal
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt
nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient
enddo; enddo
nDefNorm = nDefNorm + nDef(i,j)*nDef(i,j) ! compute the norm of the mismatch tensor
enddo; enddo
nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity)
nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces)
!--------------------------------------------------------------------------------------------------
! debuggin the mismatch tensor
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a20,i2,1x,a20,1x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb
do i = 1,3
write(6,'(1x,3(e11.4,1x))')(nDef(i,j), j = 1,3)
enddo
write(6,'(1x,a20,e11.4)')'with magnitude: ',nDefNorm
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! compute the stress penalty of all interfaces
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha &
*surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) &
*cosh(prm%ciAlpha*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) &
*tanh(nDefNorm/xSmoo_RGC)
enddo; enddo
enddo; enddo
enddo
!--------------------------------------------------------------------------------------------------
! debugging the stress-like penalty
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain
do i = 1,3
write(6,'(1x,3(e11.4,1x))')(rPen(i,j,iGrain), j = 1,3)
enddo
!$OMP END CRITICAL (write2out)
endif
enddo
end associate
end subroutine stressPenalty
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to volume discrepancy
!--------------------------------------------------------------------------------------------------
subroutine volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el)
use debug, only: &
debug_level, &
debug_homogenization,&
debug_levelExtensive, &
debug_e, &
debug_i
use mesh, only: &
mesh_homogenizationAt
use math, only: &
math_det33, &
math_inv33
use material, only: &
homogenization_maxNgrains,&
homogenization_Ngrains
use numerics, only: &
maxVolDiscr_RGC,&
volDiscrMod_RGC,&
volDiscrPow_RGC
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: vPen ! stress-like penalty due to volume
real(pReal), intent(out) :: vDiscrep ! total volume discrepancy
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef ! deformation gradients
real(pReal), dimension (3,3), intent(in) :: fAvg ! overall deformation gradient
integer(pInt), intent(in) :: ip,& ! integration point
el
real(pReal), dimension (homogenization_maxNgrains) :: gVol
integer(pInt) :: iGrain,nGrain,i,j
nGrain = homogenization_Ngrains(mesh_homogenizationAt(el))
!--------------------------------------------------------------------------------------------------
! compute the volumes of grains and of cluster
vDiscrep = math_det33(fAvg) ! compute the volume of the cluster
do iGrain = 1_pInt,nGrain
gVol(iGrain) = math_det33(fDef(1:3,1:3,iGrain)) ! compute the volume of individual grains
vDiscrep = vDiscrep - gVol(iGrain)/real(nGrain,pReal) ! calculate the difference/dicrepancy between
! the volume of the cluster and the the total volume of grains
enddo
!--------------------------------------------------------------------------------------------------
! calculate the stress and penalty due to volume discrepancy
vPen = 0.0_pReal
do iGrain = 1_pInt,nGrain
vPen(:,:,iGrain) = -1.0_pReal/real(nGrain,pReal)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* &
sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* &
gVol(iGrain)*transpose(math_inv33(fDef(:,:,iGrain)))
!--------------------------------------------------------------------------------------------------
! debugging the stress-like penalty
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30,i2)')'Volume penalty of grain: ',iGrain
do i = 1,3
write(6,'(1x,3(e11.4,1x))')(vPen(i,j,iGrain), j = 1,3)
enddo
!$OMP END CRITICAL (write2out)
endif
enddo
end subroutine volumePenalty
!--------------------------------------------------------------------------------------------------
!> @brief compute the correction factor accouted for surface evolution (area change) due to
! deformation
!--------------------------------------------------------------------------------------------------
function surfaceCorrection(avgF,ip,el)
use math, only: &
math_invert33, &
math_mul33x33
implicit none
real(pReal), dimension(3) :: surfaceCorrection
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
integer(pInt), intent(in) :: ip,& !< integration point number
el !< element number
real(pReal), dimension(3,3) :: invC
real(pReal), dimension(3) :: nVect
real(pReal) :: detF
integer(pInt) :: i,j,iBase
logical :: error
call math_invert33(math_mul33x33(transpose(avgF),avgF),invC,detF,error)
surfaceCorrection = 0.0_pReal
do iBase = 1_pInt,3_pInt
nVect = interfaceNormal([iBase,1_pInt,1_pInt,1_pInt],ip,el) ! get the normal of the interface
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal
enddo; enddo
surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF ! get the surface correction factor (area contraction/enlargement)
enddo
end function surfaceCorrection
!--------------------------------------------------------------------------------------------------
!> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
!--------------------------------------------------------------------------------------------------
function equivalentModuli(grainID,ip,el)
use constitutive, only: &
constitutive_homogenizedC
implicit none
integer(pInt), intent(in) :: &
grainID,&
ip, & !< integration point number
el !< element number
real(pReal), dimension (6,6) :: elasTens
real(pReal), dimension(2) :: equivalentModuli
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
!--------------------------------------------------------------------------------------------------
!> @brief calculating the grain deformation gradient (the same with
! homogenization_RGC_partionDeformation, but used only for perturbation scheme)
!--------------------------------------------------------------------------------------------------
subroutine grainDeformation(F, avgF, ip, el)
use mesh, only: &
mesh_homogenizationAt
use material, only: &
homogenization_maxNgrains,&
homogenization_typeInstance
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain
real(pReal), dimension (3,3), intent(in) :: avgF !<
integer(pInt), intent(in) :: &
el, & !< element number
ip !< integration point number
real(pReal), dimension (3) :: aVect,nVect
integer(pInt), dimension (4) :: intFace
integer(pInt), dimension (3) :: iGrain3
integer(pInt) :: instance, iGrain,iFace,i,j
!--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
instance = homogenization_typeInstance(mesh_homogenizationAt(el))
F = 0.0_pReal
do iGrain = 1_pInt,product(param(instance)%Nconstituents)
iGrain3 = grain1to3(iGrain,instance)
do iFace = 1_pInt,6_pInt
intFace = getInterface(iFace,iGrain3)
aVect = relaxationVector(intFace,instance, ip, el)
nVect = interfaceNormal(intFace,ip,el)
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations
enddo
F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! relaxed deformation gradient
enddo
end subroutine grainDeformation
end function homogenization_RGC_updateState end function homogenization_RGC_updateState
@ -887,285 +1204,7 @@ pure function homogenization_RGC_postResults(ip,el) result(postResults)
end function homogenization_RGC_postResults end function homogenization_RGC_postResults
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to deformation mismatch
!--------------------------------------------------------------------------------------------------
subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance)
use debug, only: &
debug_level, &
debug_homogenization,&
debug_levelExtensive, &
debug_e, &
debug_i
use mesh, only: &
mesh_homogenizationAt
use constitutive, only: &
constitutive_homogenizedC
use math, only: &
math_civita
use material, only: &
homogenization_maxNgrains,&
homogenization_Ngrains
use numerics, only: &
xSmoo_RGC
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen !< stress-like penalty
real(pReal), dimension (3,homogenization_maxNgrains), intent(out) :: nMis !< total amount of mismatch
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef !< deformation gradients
real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor
integer(pInt), intent(in) :: ip,el,instance
integer(pInt), dimension (4) :: intFace
integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim
real(pReal), dimension (3,3) :: gDef,nDef
real(pReal), dimension (3) :: nVect,surfCorr
real(pReal), dimension (2) :: Gmoduli
integer(pInt) :: iGrain,iGNghb,iFace,i,j,k,l
real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
type(tParameters) :: prm
integer(pInt), parameter :: nFace = 6_pInt
real(pReal), parameter :: nDefToler = 1.0e-10_pReal
nGDim = param(instance)%Nconstituents
rPen = 0.0_pReal
nMis = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! get the correction factor the modulus of penalty stress representing the evolution of area of
! the interfaces due to deformations
surfCorr = surfaceCorrection(avgF,ip,el)
associate(prm => param(instance))
!--------------------------------------------------------------------------------------------------
! debugging the surface correction factor
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el
write(6,'(1x,3(e11.4,1x))')(surfCorr(i), i = 1,3)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! computing the mismatch and penalty stress tensor of all grains
do iGrain = 1_pInt,homogenization_Ngrains(mesh_homogenizationAt(el))
Gmoduli = equivalentModuli(iGrain,ip,el)
muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain
bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector
iGrain3 = grain1to3(iGrain,instance) ! get the grain ID in local 3-dimensional index (x,y,z)-position
!* Looping over all six interfaces of each grain
do iFace = 1_pInt,nFace
intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain
nVect = interfaceNormal(intFace,ip,el) ! get the interface normal
iGNghb3 = iGrain3 ! identify the neighboring grain across the interface
iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal),pInt)
if (iGNghb3(1) < 1) iGNghb3(1) = nGDim(1) ! with periodicity along e1 direction
if (iGNghb3(1) > nGDim(1)) iGNghb3(1) = 1_pInt
if (iGNghb3(2) < 1) iGNghb3(2) = nGDim(2) ! with periodicity along e2 direction
if (iGNghb3(2) > nGDim(2)) iGNghb3(2) = 1_pInt
if (iGNghb3(3) < 1) iGNghb3(3) = nGDim(3) ! with periodicity along e3 direction
if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1_pInt
iGNghb = grain3to1(iGNghb3,instance) ! get the ID of the neighboring grain
Gmoduli = equivalentModuli(iGNghb,ip,el) ! collecting the shear modulus and Burgers vector of the neighbor
muGNghb = Gmoduli(1)
bgGNghb = Gmoduli(2)
gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! compute the difference/jump in deformation gradeint across the neighbor
!--------------------------------------------------------------------------------------------------
! compute the mismatch tensor of all interfaces
nDefNorm = 0.0_pReal
nDef = 0.0_pReal
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt
nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient
enddo; enddo
nDefNorm = nDefNorm + nDef(i,j)*nDef(i,j) ! compute the norm of the mismatch tensor
enddo; enddo
nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity)
nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces)
!--------------------------------------------------------------------------------------------------
! debuggin the mismatch tensor
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a20,i2,1x,a20,1x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb
do i = 1,3
write(6,'(1x,3(e11.4,1x))')(nDef(i,j), j = 1,3)
enddo
write(6,'(1x,a20,e11.4)')'with magnitude: ',nDefNorm
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! compute the stress penalty of all interfaces
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha &
*surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) &
*cosh(prm%ciAlpha*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) &
*tanh(nDefNorm/xSmoo_RGC)
enddo; enddo
enddo; enddo
enddo
!--------------------------------------------------------------------------------------------------
! debugging the stress-like penalty
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain
do i = 1,3
write(6,'(1x,3(e11.4,1x))')(rPen(i,j,iGrain), j = 1,3)
enddo
!$OMP END CRITICAL (write2out)
endif
enddo
end associate
end subroutine stressPenalty
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to volume discrepancy
!--------------------------------------------------------------------------------------------------
subroutine volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el)
use debug, only: &
debug_level, &
debug_homogenization,&
debug_levelExtensive, &
debug_e, &
debug_i
use mesh, only: &
mesh_homogenizationAt
use math, only: &
math_det33, &
math_inv33
use material, only: &
homogenization_maxNgrains,&
homogenization_Ngrains
use numerics, only: &
maxVolDiscr_RGC,&
volDiscrMod_RGC,&
volDiscrPow_RGC
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: vPen ! stress-like penalty due to volume
real(pReal), intent(out) :: vDiscrep ! total volume discrepancy
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef ! deformation gradients
real(pReal), dimension (3,3), intent(in) :: fAvg ! overall deformation gradient
integer(pInt), intent(in) :: ip,& ! integration point
el
real(pReal), dimension (homogenization_maxNgrains) :: gVol
integer(pInt) :: iGrain,nGrain,i,j
nGrain = homogenization_Ngrains(mesh_homogenizationAt(el))
!--------------------------------------------------------------------------------------------------
! compute the volumes of grains and of cluster
vDiscrep = math_det33(fAvg) ! compute the volume of the cluster
do iGrain = 1_pInt,nGrain
gVol(iGrain) = math_det33(fDef(1:3,1:3,iGrain)) ! compute the volume of individual grains
vDiscrep = vDiscrep - gVol(iGrain)/real(nGrain,pReal) ! calculate the difference/dicrepancy between
! the volume of the cluster and the the total volume of grains
enddo
!--------------------------------------------------------------------------------------------------
! calculate the stress and penalty due to volume discrepancy
vPen = 0.0_pReal
do iGrain = 1_pInt,nGrain
vPen(:,:,iGrain) = -1.0_pReal/real(nGrain,pReal)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* &
sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* &
gVol(iGrain)*transpose(math_inv33(fDef(:,:,iGrain)))
!--------------------------------------------------------------------------------------------------
! debugging the stress-like penalty
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30,i2)')'Volume penalty of grain: ',iGrain
do i = 1,3
write(6,'(1x,3(e11.4,1x))')(vPen(i,j,iGrain), j = 1,3)
enddo
!$OMP END CRITICAL (write2out)
endif
enddo
end subroutine volumePenalty
!--------------------------------------------------------------------------------------------------
!> @brief compute the correction factor accouted for surface evolution (area change) due to
! deformation
!--------------------------------------------------------------------------------------------------
function surfaceCorrection(avgF,ip,el)
use math, only: &
math_invert33, &
math_mul33x33
implicit none
real(pReal), dimension(3) :: surfaceCorrection
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
integer(pInt), intent(in) :: ip,& !< integration point number
el !< element number
real(pReal), dimension(3,3) :: invC
real(pReal), dimension(3) :: nVect
real(pReal) :: detF
integer(pInt) :: i,j,iBase
logical :: error
call math_invert33(math_mul33x33(transpose(avgF),avgF),invC,detF,error)
surfaceCorrection = 0.0_pReal
do iBase = 1_pInt,3_pInt
nVect = interfaceNormal([iBase,1_pInt,1_pInt,1_pInt],ip,el) ! get the normal of the interface
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal
enddo; enddo
surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF ! get the surface correction factor (area contraction/enlargement)
enddo
end function surfaceCorrection
!--------------------------------------------------------------------------------------------------
!> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
!--------------------------------------------------------------------------------------------------
function equivalentModuli(grainID,ip,el)
use constitutive, only: &
constitutive_homogenizedC
implicit none
integer(pInt), intent(in) :: &
grainID,&
ip, & !< integration point number
el !< element number
real(pReal), dimension (6,6) :: elasTens
real(pReal), dimension(2) :: equivalentModuli
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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1394,44 +1433,4 @@ pure function interface1to4(iFace1D, instance)
end function interface1to4 end function interface1to4
!--------------------------------------------------------------------------------------------------
!> @brief calculating the grain deformation gradient (the same with
! homogenization_RGC_partionDeformation, but used only for perturbation scheme)
!--------------------------------------------------------------------------------------------------
subroutine grainDeformation(F, avgF, ip, el)
use mesh, only: &
mesh_homogenizationAt
use material, only: &
homogenization_maxNgrains,&
homogenization_typeInstance
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain
real(pReal), dimension (3,3), intent(in) :: avgF !<
integer(pInt), intent(in) :: &
el, & !< element number
ip !< integration point number
real(pReal), dimension (3) :: aVect,nVect
integer(pInt), dimension (4) :: intFace
integer(pInt), dimension (3) :: iGrain3
integer(pInt) :: instance, iGrain,iFace,i,j
!--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
instance = homogenization_typeInstance(mesh_homogenizationAt(el))
F = 0.0_pReal
do iGrain = 1_pInt,product(param(instance)%Nconstituents)
iGrain3 = grain1to3(iGrain,instance)
do iFace = 1_pInt,6_pInt
intFace = getInterface(iFace,iGrain3)
aVect = relaxationVector(intFace,instance, ip, el)
nVect = interfaceNormal(intFace,ip,el)
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations
enddo
F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! relaxed deformation gradient
enddo
end subroutine grainDeformation
end module homogenization_RGC end module homogenization_RGC