exploit symmetry for stress calculation

This commit is contained in:
Martin Diehl 2021-11-20 22:36:01 +01:00
parent a857285e34
commit 28ff027b38
2 changed files with 37 additions and 6 deletions

View File

@ -841,14 +841,14 @@ pure function math_66toSym3333(m66,weighted)
w = merge(INVNRMMANDEL,1.0_pReal,weighted)
else
w = INVNRMMANDEL
endif
end if
do i=1,6; do j=1,6
math_66toSym3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) = w(i)*w(j)*m66(i,j)
math_66toSym3333(MAPNYE(2,i),MAPNYE(1,i),MAPNYE(1,j),MAPNYE(2,j)) = w(i)*w(j)*m66(i,j)
math_66toSym3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(2,j),MAPNYE(1,j)) = w(i)*w(j)*m66(i,j)
math_66toSym3333(MAPNYE(2,i),MAPNYE(1,i),MAPNYE(2,j),MAPNYE(1,j)) = w(i)*w(j)*m66(i,j)
enddo; enddo
end do; end do
end function math_66toSym3333
@ -866,7 +866,6 @@ pure function math_Voigt6to33_stress(sigma_tilde) result(sigma)
sigma_tilde(6), sigma_tilde(2), sigma_tilde(4), &
sigma_tilde(5), sigma_tilde(4), sigma_tilde(3)],[3,3])
end function math_Voigt6to33_stress
@ -883,10 +882,40 @@ pure function math_Voigt6to33_strain(epsilon_tilde) result(epsilon)
0.5_pReal*epsilon_tilde(6), epsilon_tilde(2), 0.5_pReal*epsilon_tilde(4), &
0.5_pReal*epsilon_tilde(5), 0.5_pReal*epsilon_tilde(4), epsilon_tilde(3)],[3,3])
end function math_Voigt6to33_strain
!--------------------------------------------------------------------------------------------------
!> @brief Convert 3x3 tensor into 6 Voigt stress vector.
!--------------------------------------------------------------------------------------------------
pure function math_33toVoigt6_stress(sigma) result(sigma_tilde)
real(pReal), dimension(6) :: sigma_tilde
real(pReal), dimension(3,3), intent(in) :: sigma
sigma_tilde = [sigma(1,1), sigma(2,2), sigma(3,3), &
sigma(3,2), sigma(3,1), sigma(1,2)]
end function math_33toVoigt6_stress
!--------------------------------------------------------------------------------------------------
!> @brief Convert 3x3 tensor into 6 Voigt strain vector.
!--------------------------------------------------------------------------------------------------
pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde)
real(pReal), dimension(6) :: epsilon_tilde
real(pReal), dimension(3,3), intent(in) :: epsilon
epsilon_tilde = [ epsilon(1,1), epsilon(2,2), epsilon(3,3), &
2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)]
end function math_33toVoigt6_strain
!--------------------------------------------------------------------------------------------------
!> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix.
!--------------------------------------------------------------------------------------------------

View File

@ -148,15 +148,17 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
real(pReal), dimension(3,3) :: E
real(pReal), dimension(6,6) :: C66
real(pReal), dimension(3,3,3,3) :: C
integer :: &
i, j
C = math_Voigt66to3333(phase_damage_C66(phase_homogenizedC66(ph,en),ph,en))
C66 = phase_damage_C66(phase_homogenizedC66(ph,en),ph,en)
C = math_Voigt66to3333(C66)
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
S = math_Voigt6to33_stress(matmul(C66,math_33toVoigt6_strain(matmul(matmul(transpose(Fi),E),Fi))))!< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
do i =1,3; do j=1,3
dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko