Merge branch 'development' into 38-introduce-rudimentary-PETSc-based-FEM-solver
This commit is contained in:
Martin Diehl
commit
cc262ae198
|
|
@ -5,15 +5,10 @@
|
|||
(output) orientation # quaternion
|
||||
(output) eulerangles # orientation as Bunge triple in degree
|
||||
(output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) in crystal reference coordinates
|
||||
(output) grainrotationx # deviation from initial orientation as angle in degrees around sample reference x axis
|
||||
(output) grainrotationy # deviation from initial orientation as angle in degrees around sample reference y axis
|
||||
(output) grainrotationz # deviation from initial orientation as angle in degrees around sample reference z axis
|
||||
(output) f # deformation gradient tensor; synonyms: "defgrad"
|
||||
(output) f # deformation gradient tensor
|
||||
(output) fe # elastic deformation gradient tensor
|
||||
(output) fp # plastic deformation gradient tensor
|
||||
(output) e # total strain as Green-Lagrange tensor
|
||||
(output) ee # elastic strain as Green-Lagrange tensor
|
||||
(output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola"
|
||||
(output) s # second Piola-Kichhoff stress tensor; synonyms: "tstar", "secondpiola", "2ndpiola"
|
||||
(output) p # first Piola-Kichhoff stress tensor
|
||||
(output) s # second Piola-Kichhoff stress tensor
|
||||
(output) lp # plastic velocity gradient tensor
|
||||
(output) elasmatrix # elastic stiffness matrix
|
||||
|
|
|
|||
|
|
@ -18,8 +18,6 @@ mech none
|
|||
(output) f # deformation gradient tensor; synonyms: "defgrad"
|
||||
(output) fe # elastic deformation gradient tensor
|
||||
(output) fp # plastic deformation gradient tensor
|
||||
(output) e # total strain as Green-Lagrange tensor
|
||||
(output) ee # elastic strain as Green-Lagrange tensor
|
||||
(output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola"
|
||||
(output) lp # plastic velocity gradient tensor
|
||||
|
||||
|
|
|
|||
|
|
@ -19,16 +19,16 @@ module constitutive
|
|||
constitutive_init, &
|
||||
constitutive_homogenizedC, &
|
||||
constitutive_microstructure, &
|
||||
constitutive_LpAndItsTangent, &
|
||||
constitutive_LiAndItsTangent, &
|
||||
constitutive_LpAndItsTangents, &
|
||||
constitutive_LiAndItsTangents, &
|
||||
constitutive_initialFi, &
|
||||
constitutive_TandItsTangent, &
|
||||
constitutive_SandItsTangents, &
|
||||
constitutive_collectDotState, &
|
||||
constitutive_collectDeltaState, &
|
||||
constitutive_postResults
|
||||
|
||||
private :: &
|
||||
constitutive_hooke_TandItsTangent
|
||||
constitutive_hooke_SandItsTangents
|
||||
|
||||
contains
|
||||
|
||||
|
|
@ -346,6 +346,7 @@ end subroutine constitutive_init
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns the homogenize elasticity matrix
|
||||
!> ToDo: homogenizedC66 would be more consistent
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function constitutive_homogenizedC(ipc,ip,el)
|
||||
use prec, only: &
|
||||
|
|
@ -430,7 +431,7 @@ end subroutine constitutive_microstructure
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the velocity gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v, Fi, ipc, ip, el)
|
||||
subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal
|
||||
use math, only: &
|
||||
|
|
@ -470,20 +471,21 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
|
|||
ip, & !< integration point
|
||||
el !< element
|
||||
real(pReal), intent(in), dimension(6) :: &
|
||||
Tstar_v !< 2nd Piola-Kirchhoff stress
|
||||
S6 !< 2nd Piola-Kirchhoff stress (vector notation)
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
Fi !< intermediate deformation gradient
|
||||
real(pReal), intent(out), dimension(3,3) :: &
|
||||
Lp !< plastic velocity gradient
|
||||
real(pReal), intent(out), dimension(3,3,3,3) :: &
|
||||
dLp_dTstar3333, & !< derivative of Lp with respect to Tstar (4th-order tensor)
|
||||
dLp_dFi3333 !< derivative of Lp with respect to Fi (4th-order tensor)
|
||||
real(pReal), dimension(6) :: &
|
||||
Mstar_v !< Mandel stress work conjugate with Lp
|
||||
dLp_dS, &
|
||||
dLp_dFi !< derivative of Lp with respect to Fi
|
||||
real(pReal), dimension(3,3,3,3) :: &
|
||||
dLp_dMp !< derivative of Lp with respect to Mandel stress
|
||||
real(pReal), dimension(9,9) :: &
|
||||
dLp_dMstar !< derivative of Lp with respect to Mstar (4th-order tensor)
|
||||
dLp_dMp99 !< derivative of Lp with respect to Mstar (matrix notation)
|
||||
real(pReal), dimension(3,3) :: &
|
||||
temp_33
|
||||
Mp, & !< Mandel stress work conjugate with Lp
|
||||
S !< 2nd Piola-Kirchhoff stress
|
||||
integer(pInt) :: &
|
||||
ho, & !< homogenization
|
||||
tme !< thermal member position
|
||||
|
|
@ -493,47 +495,57 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
|
|||
ho = material_homog(ip,el)
|
||||
tme = thermalMapping(ho)%p(ip,el)
|
||||
|
||||
Mstar_v = math_Mandel33to6(math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar_v)))
|
||||
S = math_Mandel6to33(S6)
|
||||
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
|
||||
|
||||
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
|
||||
|
||||
case (PLASTICITY_NONE_ID) plasticityType
|
||||
Lp = 0.0_pReal
|
||||
dLp_dMstar = 0.0_pReal
|
||||
dLp_dMp = 0.0_pReal
|
||||
|
||||
case (PLASTICITY_ISOTROPIC_ID) plasticityType
|
||||
call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
|
||||
call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
|
||||
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
|
||||
|
||||
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
|
||||
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
|
||||
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
|
||||
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
|
||||
|
||||
case (PLASTICITY_KINEHARDENING_ID) plasticityType
|
||||
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
|
||||
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
|
||||
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
|
||||
|
||||
case (PLASTICITY_NONLOCAL_ID) plasticityType
|
||||
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, &
|
||||
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
|
||||
temperature(ho)%p(tme),ip,el)
|
||||
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
|
||||
|
||||
case (PLASTICITY_DISLOTWIN_ID) plasticityType
|
||||
call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, &
|
||||
call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
|
||||
temperature(ho)%p(tme),ipc,ip,el)
|
||||
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
|
||||
|
||||
case (PLASTICITY_DISLOUCLA_ID) plasticityType
|
||||
call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, &
|
||||
call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
|
||||
temperature(ho)%p(tme), ipc,ip,el)
|
||||
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
|
||||
|
||||
end select plasticityType
|
||||
|
||||
dLp_dTstar3333 = math_Plain99to3333(dLp_dMstar)
|
||||
temp_33 = math_mul33x33(Fi,math_Mandel6to33(Tstar_v))
|
||||
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) &
|
||||
dLp_dFi3333(i,j,1:3,1:3) = math_mul33x33(temp_33,transpose(dLp_dTstar3333(i,j,1:3,1:3))) + &
|
||||
math_mul33x33(math_mul33x33(Fi,dLp_dTstar3333(i,j,1:3,1:3)),math_Mandel6to33(Tstar_v))
|
||||
do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
|
||||
dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
|
||||
math_mul33x33(math_mul33x33(Fi,dLp_dMp(i,j,1:3,1:3)),S)
|
||||
dLp_dS(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
|
||||
enddo
|
||||
|
||||
temp_33 = math_mul33x33(transpose(Fi),Fi)
|
||||
|
||||
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) &
|
||||
dLp_dTstar3333(i,j,1:3,1:3) = math_mul33x33(temp_33,dLp_dTstar3333(i,j,1:3,1:3))
|
||||
|
||||
end subroutine constitutive_LpAndItsTangent
|
||||
end subroutine constitutive_LpAndItsTangents
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the velocity gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v, Fi, ipc, ip, el)
|
||||
subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal
|
||||
use math, only: &
|
||||
|
|
@ -571,18 +583,18 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v
|
|||
ip, & !< integration point
|
||||
el !< element
|
||||
real(pReal), intent(in), dimension(6) :: &
|
||||
Tstar_v !< 2nd Piola-Kirchhoff stress
|
||||
S6 !< 2nd Piola-Kirchhoff stress (vector notation)
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
Fi !< intermediate deformation gradient
|
||||
real(pReal), intent(out), dimension(3,3) :: &
|
||||
Li !< intermediate velocity gradient
|
||||
real(pReal), intent(out), dimension(3,3,3,3) :: &
|
||||
dLi_dTstar3333, & !< derivative of Li with respect to Tstar (4th-order tensor)
|
||||
dLi_dFi3333
|
||||
dLi_dS, & !< derivative of Li with respect to S
|
||||
dLi_dFi
|
||||
real(pReal), dimension(3,3) :: &
|
||||
my_Li !< intermediate velocity gradient
|
||||
real(pReal), dimension(3,3,3,3) :: &
|
||||
my_dLi_dTstar
|
||||
my_dLi_dS
|
||||
real(pReal), dimension(3,3) :: &
|
||||
FiInv, &
|
||||
temp_33
|
||||
|
|
@ -594,51 +606,51 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v
|
|||
i, j
|
||||
|
||||
Li = 0.0_pReal
|
||||
dLi_dTstar3333 = 0.0_pReal
|
||||
dLi_dFi3333 = 0.0_pReal
|
||||
dLi_dS = 0.0_pReal
|
||||
dLi_dFi = 0.0_pReal
|
||||
|
||||
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
|
||||
case (PLASTICITY_isotropic_ID) plasticityType
|
||||
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el)
|
||||
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
|
||||
case default plasticityType
|
||||
my_Li = 0.0_pReal
|
||||
my_dLi_dTstar = 0.0_pReal
|
||||
my_dLi_dS = 0.0_pReal
|
||||
end select plasticityType
|
||||
|
||||
Li = Li + my_Li
|
||||
dLi_dTstar3333 = dLi_dTstar3333 + my_dLi_dTstar
|
||||
dLi_dS = dLi_dS + my_dLi_dS
|
||||
|
||||
KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el))
|
||||
kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el)))
|
||||
case (KINEMATICS_cleavage_opening_ID) kinematicsType
|
||||
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el)
|
||||
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
|
||||
case (KINEMATICS_slipplane_opening_ID) kinematicsType
|
||||
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el)
|
||||
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
|
||||
case (KINEMATICS_thermal_expansion_ID) kinematicsType
|
||||
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el)
|
||||
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
|
||||
case (KINEMATICS_vacancy_strain_ID) kinematicsType
|
||||
call kinematics_vacancy_strain_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el)
|
||||
call kinematics_vacancy_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
|
||||
case (KINEMATICS_hydrogen_strain_ID) kinematicsType
|
||||
call kinematics_hydrogen_strain_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el)
|
||||
call kinematics_hydrogen_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
|
||||
case default kinematicsType
|
||||
my_Li = 0.0_pReal
|
||||
my_dLi_dTstar = 0.0_pReal
|
||||
my_dLi_dS = 0.0_pReal
|
||||
end select kinematicsType
|
||||
Li = Li + my_Li
|
||||
dLi_dTstar3333 = dLi_dTstar3333 + my_dLi_dTstar
|
||||
dLi_dS = dLi_dS + my_dLi_dS
|
||||
enddo KinematicsLoop
|
||||
|
||||
FiInv = math_inv33(Fi)
|
||||
detFi = math_det33(Fi)
|
||||
Li = math_mul33x33(math_mul33x33(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration
|
||||
temp_33 = math_mul33x33(FiInv,Li)
|
||||
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
|
||||
dLi_dTstar3333(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dTstar3333(1:3,1:3,i,j)),FiInv)*detFi
|
||||
dLi_dFi3333 (1:3,1:3,i,j) = dLi_dFi3333(1:3,1:3,i,j) + Li*FiInv(j,i)
|
||||
dLi_dFi3333 (1:3,i,1:3,j) = dLi_dFi3333(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
|
||||
end forall
|
||||
do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
|
||||
dLi_dS(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
|
||||
dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i)
|
||||
dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
|
||||
end do
|
||||
|
||||
end subroutine constitutive_LiAndItsTangent
|
||||
end subroutine constitutive_LiAndItsTangents
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -696,10 +708,10 @@ end function constitutive_initialFi
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
|
||||
!> the elastic deformation gradient depending on the selected elastic law (so far no case switch
|
||||
!! because only Hooke is implemented
|
||||
!> the elastic/intermediate deformation gradients depending on the selected elastic law
|
||||
!! (so far no case switch because only Hooke is implemented)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine constitutive_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
|
||||
subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal
|
||||
|
||||
|
|
@ -712,30 +724,28 @@ subroutine constitutive_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
|
|||
Fe, & !< elastic deformation gradient
|
||||
Fi !< intermediate deformation gradient
|
||||
real(pReal), intent(out), dimension(3,3) :: &
|
||||
T !< 2nd Piola-Kirchhoff stress tensor
|
||||
S !< 2nd Piola-Kirchhoff stress tensor
|
||||
real(pReal), intent(out), dimension(3,3,3,3) :: &
|
||||
dT_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
|
||||
dT_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
|
||||
dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
|
||||
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
|
||||
|
||||
call constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
|
||||
call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
|
||||
|
||||
|
||||
end subroutine constitutive_TandItsTangent
|
||||
end subroutine constitutive_SandItsTangents
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
|
||||
!> the elastic deformation gradient using hookes law
|
||||
!> the elastic and intermeidate deformation gradients using Hookes law
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
|
||||
subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal
|
||||
use math, only : &
|
||||
math_mul3x3, &
|
||||
math_mul33x33, &
|
||||
math_mul3333xx33, &
|
||||
math_Mandel66to3333, &
|
||||
math_trace33, &
|
||||
math_I3
|
||||
use material, only: &
|
||||
material_phase, &
|
||||
|
|
@ -758,10 +768,10 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip,
|
|||
Fe, & !< elastic deformation gradient
|
||||
Fi !< intermediate deformation gradient
|
||||
real(pReal), intent(out), dimension(3,3) :: &
|
||||
T !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
|
||||
S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
|
||||
real(pReal), intent(out), dimension(3,3,3,3) :: &
|
||||
dT_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
|
||||
dT_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
|
||||
dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
|
||||
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
|
||||
real(pReal), dimension(3,3) :: E
|
||||
real(pReal), dimension(3,3,3,3) :: C
|
||||
integer(pInt) :: &
|
||||
|
|
@ -784,22 +794,22 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip,
|
|||
enddo DegradationLoop
|
||||
|
||||
E = 0.5_pReal*(math_mul33x33(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
|
||||
T = math_mul3333xx33(C,math_mul33x33(math_mul33x33(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
|
||||
S = math_mul3333xx33(C,math_mul33x33(math_mul33x33(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
|
||||
|
||||
dT_dFe = 0.0_pReal
|
||||
dS_dFe = 0.0_pReal
|
||||
forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt)
|
||||
dT_dFe(i,j,1:3,1:3) = &
|
||||
math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dT_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
|
||||
dT_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3)) !< dT_ij/dFi_kl = C_ijln * E_km * Fe_mn
|
||||
dS_dFe(i,j,1:3,1:3) = &
|
||||
math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
|
||||
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
|
||||
end forall
|
||||
|
||||
end subroutine constitutive_hooke_TandItsTangent
|
||||
end subroutine constitutive_hooke_SandItsTangents
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfracArray,ipc, ip, el)
|
||||
subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfracArray,ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal, &
|
||||
pLongInt
|
||||
|
|
@ -807,6 +817,10 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
|
|||
debug_level, &
|
||||
debug_constitutive, &
|
||||
debug_levelBasic
|
||||
use math, only: &
|
||||
math_Mandel6to33, &
|
||||
math_Mandel33to6, &
|
||||
math_mul33x33
|
||||
use mesh, only: &
|
||||
mesh_NcpElems, &
|
||||
mesh_maxNips
|
||||
|
|
@ -863,8 +877,12 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
|
|||
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
|
||||
FeArray, & !< elastic deformation gradient
|
||||
FpArray !< plastic deformation gradient
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
Fi !< intermediate deformation gradient
|
||||
real(pReal), intent(in), dimension(6) :: &
|
||||
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
|
||||
S6 !< 2nd Piola Kirchhoff stress (vector notation)
|
||||
real(pReal), dimension(3,3) :: &
|
||||
Mstar
|
||||
integer(pInt) :: &
|
||||
ho, & !< homogenization
|
||||
tme, & !< thermal member position
|
||||
|
|
@ -873,35 +891,50 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
|
|||
ho = material_homog( ip,el)
|
||||
tme = thermalMapping(ho)%p(ip,el)
|
||||
|
||||
Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
|
||||
|
||||
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
|
||||
|
||||
case (PLASTICITY_ISOTROPIC_ID) plasticityType
|
||||
call plastic_isotropic_dotState (Tstar_v,ipc,ip,el)
|
||||
call plastic_isotropic_dotState (math_Mandel33to6(Mstar),ipc,ip,el)
|
||||
|
||||
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
|
||||
call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
|
||||
call plastic_phenopowerlaw_dotState(math_Mandel33to6(Mstar),ipc,ip,el)
|
||||
|
||||
case (PLASTICITY_KINEHARDENING_ID) plasticityType
|
||||
call plastic_kinehardening_dotState(Tstar_v,ipc,ip,el)
|
||||
call plastic_kinehardening_dotState(math_Mandel33to6(Mstar),ipc,ip,el)
|
||||
|
||||
case (PLASTICITY_DISLOTWIN_ID) plasticityType
|
||||
call plastic_dislotwin_dotState (Tstar_v,temperature(ho)%p(tme), &
|
||||
call plastic_dislotwin_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), &
|
||||
ipc,ip,el)
|
||||
|
||||
case (PLASTICITY_DISLOUCLA_ID) plasticityType
|
||||
call plastic_disloucla_dotState (Tstar_v,temperature(ho)%p(tme), &
|
||||
call plastic_disloucla_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), &
|
||||
ipc,ip,el)
|
||||
|
||||
case (PLASTICITY_NONLOCAL_ID) plasticityType
|
||||
call plastic_nonlocal_dotState (Tstar_v,FeArray,FpArray,temperature(ho)%p(tme), &
|
||||
call plastic_nonlocal_dotState (math_Mandel33to6(Mstar),FeArray,FpArray,temperature(ho)%p(tme), &
|
||||
subdt,subfracArray,ip,el)
|
||||
end select plasticityType
|
||||
|
||||
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
|
||||
sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
|
||||
|
||||
sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
|
||||
|
||||
case (SOURCE_damage_anisoBrittle_ID) sourceType
|
||||
call source_damage_anisoBrittle_dotState (Tstar_v, ipc, ip, el)
|
||||
call source_damage_anisoBrittle_dotState (S6, ipc, ip, el) !< correct stress?
|
||||
|
||||
case (SOURCE_damage_isoDuctile_ID) sourceType
|
||||
call source_damage_isoDuctile_dotState ( ipc, ip, el)
|
||||
|
||||
case (SOURCE_damage_anisoDuctile_ID) sourceType
|
||||
call source_damage_anisoDuctile_dotState ( ipc, ip, el)
|
||||
|
||||
case (SOURCE_thermal_externalheat_ID) sourceType
|
||||
call source_thermal_externalheat_dotState( ipc, ip, el)
|
||||
|
||||
end select sourceType
|
||||
|
||||
enddo SourceLoop
|
||||
|
||||
end subroutine constitutive_collectDotState
|
||||
|
|
@ -910,7 +943,7 @@ end subroutine constitutive_collectDotState
|
|||
!> @brief for constitutive models having an instantaneous change of state
|
||||
!> will return false if delta state is not needed/supported by the constitutive model
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
|
||||
subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal, &
|
||||
pLongInt
|
||||
|
|
@ -918,6 +951,10 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
|
|||
debug_level, &
|
||||
debug_constitutive, &
|
||||
debug_levelBasic
|
||||
use math, only: &
|
||||
math_Mandel6to33, &
|
||||
math_Mandel33to6, &
|
||||
math_mul33x33
|
||||
use material, only: &
|
||||
phase_plasticity, &
|
||||
phase_source, &
|
||||
|
|
@ -945,29 +982,43 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
|
|||
ip, & !< integration point
|
||||
el !< element
|
||||
real(pReal), intent(in), dimension(6) :: &
|
||||
Tstar_v !< 2nd Piola-Kirchhoff stress
|
||||
S6 !< 2nd Piola Kirchhoff stress (vector notation)
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
Fe !< elastic deformation gradient
|
||||
Fe, & !< elastic deformation gradient
|
||||
Fi !< intermediate deformation gradient
|
||||
real(pReal), dimension(3,3) :: &
|
||||
Mstar
|
||||
integer(pInt) :: &
|
||||
s !< counter in source loop
|
||||
|
||||
Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
|
||||
|
||||
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
|
||||
|
||||
case (PLASTICITY_KINEHARDENING_ID) plasticityType
|
||||
call plastic_kinehardening_deltaState(Tstar_v,ipc,ip,el)
|
||||
call plastic_kinehardening_deltaState(math_Mandel33to6(Mstar),ipc,ip,el)
|
||||
|
||||
case (PLASTICITY_NONLOCAL_ID) plasticityType
|
||||
call plastic_nonlocal_deltaState(Tstar_v,ip,el)
|
||||
call plastic_nonlocal_deltaState(math_Mandel33to6(Mstar),ip,el)
|
||||
|
||||
end select plasticityType
|
||||
|
||||
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
|
||||
sourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
|
||||
|
||||
sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
|
||||
|
||||
case (SOURCE_damage_isoBrittle_ID) sourceType
|
||||
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, &
|
||||
ipc, ip, el)
|
||||
ipc, ip, el)
|
||||
|
||||
case (SOURCE_vacancy_irradiation_ID) sourceType
|
||||
call source_vacancy_irradiation_deltaState(ipc, ip, el)
|
||||
|
||||
case (SOURCE_vacancy_thermalfluc_ID) sourceType
|
||||
call source_vacancy_thermalfluc_deltaState(ipc, ip, el)
|
||||
|
||||
end select sourceType
|
||||
|
||||
enddo SourceLoop
|
||||
|
||||
end subroutine constitutive_collectDeltaState
|
||||
|
|
@ -976,7 +1027,7 @@ end subroutine constitutive_collectDeltaState
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns array of constitutive results
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
|
||||
function constitutive_postResults(S6, FeArray, ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal
|
||||
use mesh, only: &
|
||||
|
|
@ -1036,7 +1087,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
|
|||
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
|
||||
FeArray !< elastic deformation gradient
|
||||
real(pReal), intent(in), dimension(6) :: &
|
||||
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
|
||||
S6 !< 2nd Piola Kirchhoff stress (vector notation)
|
||||
integer(pInt) :: &
|
||||
startPos, endPos
|
||||
integer(pInt) :: &
|
||||
|
|
@ -1054,22 +1105,22 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
|
|||
|
||||
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
|
||||
case (PLASTICITY_ISOTROPIC_ID) plasticityType
|
||||
constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
|
||||
constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(S6,ipc,ip,el)
|
||||
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
|
||||
constitutive_postResults(startPos:endPos) = &
|
||||
plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
|
||||
plastic_phenopowerlaw_postResults(S6,ipc,ip,el)
|
||||
case (PLASTICITY_KINEHARDENING_ID) plasticityType
|
||||
constitutive_postResults(startPos:endPos) = &
|
||||
plastic_kinehardening_postResults(Tstar_v,ipc,ip,el)
|
||||
plastic_kinehardening_postResults(S6,ipc,ip,el)
|
||||
case (PLASTICITY_DISLOTWIN_ID) plasticityType
|
||||
constitutive_postResults(startPos:endPos) = &
|
||||
plastic_dislotwin_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el)
|
||||
plastic_dislotwin_postResults(S6,temperature(ho)%p(tme),ipc,ip,el)
|
||||
case (PLASTICITY_DISLOUCLA_ID) plasticityType
|
||||
constitutive_postResults(startPos:endPos) = &
|
||||
plastic_disloucla_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el)
|
||||
plastic_disloucla_postResults(S6,temperature(ho)%p(tme),ipc,ip,el)
|
||||
case (PLASTICITY_NONLOCAL_ID) plasticityType
|
||||
constitutive_postResults(startPos:endPos) = &
|
||||
plastic_nonlocal_postResults (Tstar_v,FeArray,ip,el)
|
||||
plastic_nonlocal_postResults (S6,FeArray,ip,el)
|
||||
end select plasticityType
|
||||
|
||||
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -448,8 +448,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
|
|||
subStepSizeHomog, &
|
||||
stepIncreaseHomog, &
|
||||
nMPstate
|
||||
use math, only: &
|
||||
math_transpose33
|
||||
use FEsolving, only: &
|
||||
FEsolving_execElem, &
|
||||
FEsolving_execIP, &
|
||||
|
|
@ -496,6 +494,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
|
|||
crystallite_converged, &
|
||||
crystallite_stressAndItsTangent, &
|
||||
crystallite_orientations
|
||||
#ifdef DEBUG
|
||||
use debug, only: &
|
||||
debug_level, &
|
||||
debug_homogenization, &
|
||||
|
|
@ -504,6 +503,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
|
|||
debug_levelSelective, &
|
||||
debug_e, &
|
||||
debug_i
|
||||
#endif
|
||||
|
||||
implicit none
|
||||
real(pReal), intent(in) :: dt !< time increment
|
||||
|
|
@ -517,18 +517,16 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
|
|||
mySource, &
|
||||
myNgrains
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! initialize to starting condition
|
||||
#ifdef DEBUG
|
||||
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then
|
||||
!$OMP CRITICAL (write2out)
|
||||
write(6,'(/a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i
|
||||
|
||||
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', &
|
||||
math_transpose33(materialpoint_F0(1:3,1:3,debug_i,debug_e))
|
||||
transpose(materialpoint_F0(1:3,1:3,debug_i,debug_e))
|
||||
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F', &
|
||||
math_transpose33(materialpoint_F(1:3,1:3,debug_i,debug_e))
|
||||
!$OMP END CRITICAL (write2out)
|
||||
transpose(materialpoint_F(1:3,1:3,debug_i,debug_e))
|
||||
endif
|
||||
#endif
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! initialize restoration points of ...
|
||||
|
|
@ -608,10 +606,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
|
|||
!---------------------------------------------------------------------------------------------------
|
||||
! calculate new subStep and new subFrac
|
||||
materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e)
|
||||
!$OMP FLUSH(materialpoint_subFrac)
|
||||
materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), &
|
||||
stepIncreaseHomog*materialpoint_subStep(i,e)) ! introduce flexibility for step increase/acceleration
|
||||
!$OMP FLUSH(materialpoint_subStep)
|
||||
|
||||
steppingNeeded: if (materialpoint_subStep(i,e) > subStepMinHomog) then
|
||||
|
||||
|
|
@ -671,7 +667,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
|
|||
hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
|
||||
hydrogenfluxState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e))! ...internal hydrogen transport state
|
||||
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad
|
||||
!$OMP FLUSH(materialpoint_subF0)
|
||||
endif steppingNeeded
|
||||
|
||||
else converged
|
||||
|
|
@ -689,7 +684,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
|
|||
!$OMP END CRITICAL (setTerminallyIll)
|
||||
else ! cutback makes sense
|
||||
materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback
|
||||
!$OMP FLUSH(materialpoint_subStep)
|
||||
|
||||
#ifdef DEBUG
|
||||
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt &
|
||||
|
|
@ -752,8 +746,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
|
|||
|
||||
if (materialpoint_subStep(i,e) > subStepMinHomog) then
|
||||
materialpoint_requested(i,e) = .true.
|
||||
materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) + &
|
||||
materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) - materialpoint_F0(1:3,1:3,i,e))
|
||||
materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) &
|
||||
+ materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) &
|
||||
- materialpoint_F0(1:3,1:3,i,e))
|
||||
materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt
|
||||
materialpoint_doneAndHappy(1:2,i,e) = [.false.,.true.]
|
||||
endif
|
||||
|
|
@ -810,13 +805,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
|
|||
materialpoint_doneAndHappy(1:2,i,e) = homogenization_updateState(i,e)
|
||||
materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(1:2,i,e)) ! converged if done and happy
|
||||
endif
|
||||
!$OMP FLUSH(materialpoint_converged)
|
||||
if (materialpoint_converged(i,e)) then
|
||||
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then
|
||||
!$OMP CRITICAL (distributionMPState)
|
||||
!$OMP END CRITICAL (distributionMPState)
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
enddo IpLooping3
|
||||
enddo elementLooping3
|
||||
|
|
@ -838,9 +826,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
|
|||
enddo elementLooping4
|
||||
!$OMP END PARALLEL DO
|
||||
else
|
||||
!$OMP CRITICAL (write2out)
|
||||
write(6,'(/,a,/)') '<< HOMOG >> Material Point terminally ill'
|
||||
!$OMP END CRITICAL (write2out)
|
||||
endif
|
||||
|
||||
end subroutine materialpoint_stressAndItsTangent
|
||||
|
|
|
|||
Loading…
Reference in New Issue