overhaul of kinematics from FeFp to FeFiFp decomposition
This commit is contained in:
Pratheek Shanthraj
parent
26a609dbee
commit
f9cec3079e
|
|
@ -27,6 +27,10 @@ module constitutive
|
|||
constitutive_damagedC, &
|
||||
constitutive_microstructure, &
|
||||
constitutive_LpAndItsTangent, &
|
||||
constitutive_LiAndItsTangent, &
|
||||
constitutive_getFi, &
|
||||
constitutive_getFi0, &
|
||||
constitutive_getPartionedFi0, &
|
||||
constitutive_TandItsTangent, &
|
||||
constitutive_collectDotState, &
|
||||
constitutive_collectDeltaState, &
|
||||
|
|
@ -34,12 +38,10 @@ module constitutive
|
|||
constitutive_putLocalDamage, &
|
||||
constitutive_getDamage, &
|
||||
constitutive_getSlipDamage, &
|
||||
constitutive_getDamageStrain, &
|
||||
constitutive_getDamageDiffusion33, &
|
||||
constitutive_getAdiabaticTemperature, &
|
||||
constitutive_putAdiabaticTemperature, &
|
||||
constitutive_getTemperature, &
|
||||
constitutive_getThermalStrain, &
|
||||
constitutive_getLocalVacancyConcentration, &
|
||||
constitutive_putLocalVacancyConcentration, &
|
||||
constitutive_getVacancyConcentration, &
|
||||
|
|
@ -701,6 +703,199 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el)
|
|||
end subroutine constitutive_LpAndItsTangent
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the velocity gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, Tstar_v, Lp, ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal
|
||||
use material, only: &
|
||||
phase_damage, &
|
||||
phase_thermal, &
|
||||
material_phase, &
|
||||
LOCAL_DAMAGE_anisoBrittle_ID, &
|
||||
LOCAL_THERMAL_adiabatic_ID
|
||||
use damage_anisoBrittle, only: &
|
||||
damage_anisoBrittle_LdAndItsTangent
|
||||
use thermal_adiabatic, only: &
|
||||
thermal_adiabatic_LTAndItsTangent
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), intent(in), dimension(6) :: &
|
||||
Tstar_v !< 2nd Piola-Kirchhoff stress
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
Lp !< plastic velocity gradient
|
||||
real(pReal), intent(out), dimension(3,3) :: &
|
||||
Li !< intermediate velocity gradient
|
||||
real(pReal), intent(out), dimension(9,9) :: &
|
||||
dLi_dTstar !< derivative of Li with respect to Tstar (2nd-order tensor)
|
||||
real(pReal), dimension(3,3) :: &
|
||||
Li_temp !< intermediate velocity gradient
|
||||
real(pReal), dimension(9,9) :: &
|
||||
dLi_dTstar_temp !< derivative of Li with respect to Tstar (4th-order tensor)
|
||||
|
||||
Li = 0.0_pReal
|
||||
dLi_dTstar = 0.0_pReal
|
||||
|
||||
select case (phase_damage(material_phase(ipc,ip,el)))
|
||||
case (LOCAL_DAMAGE_anisoBrittle_ID)
|
||||
call damage_anisoBrittle_LdAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, ipc, ip, el)
|
||||
Li = Li + Li_temp
|
||||
dLi_dTstar = dLi_dTstar + dLi_dTstar_temp
|
||||
|
||||
end select
|
||||
|
||||
select case (phase_thermal(material_phase(ipc,ip,el)))
|
||||
case (LOCAL_THERMAL_adiabatic_ID)
|
||||
call thermal_adiabatic_LTAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, Lp, ipc, ip, el)
|
||||
Li = Li + Li_temp
|
||||
dLi_dTstar = dLi_dTstar + dLi_dTstar_temp
|
||||
|
||||
end select
|
||||
|
||||
end subroutine constitutive_LiAndItsTangent
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the intermediate deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function constitutive_getFi(ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal
|
||||
use math, only: &
|
||||
math_I3, &
|
||||
math_mul33x33
|
||||
use material, only: &
|
||||
phase_damage, &
|
||||
phase_thermal, &
|
||||
material_phase, &
|
||||
LOCAL_DAMAGE_anisoBrittle_ID, &
|
||||
LOCAL_THERMAL_adiabatic_ID
|
||||
use damage_anisoBrittle, only: &
|
||||
damage_anisoBrittle_getFd
|
||||
use thermal_adiabatic, only: &
|
||||
thermal_adiabatic_getFT
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), dimension(3,3) :: &
|
||||
constitutive_getFi !< intermediate deformation gradient
|
||||
|
||||
constitutive_getFi = math_I3
|
||||
|
||||
select case (phase_damage(material_phase(ipc,ip,el)))
|
||||
case (LOCAL_DAMAGE_anisoBrittle_ID)
|
||||
constitutive_getFi = math_mul33x33(constitutive_getFi,damage_anisoBrittle_getFd (ipc, ip, el))
|
||||
|
||||
end select
|
||||
|
||||
select case (phase_thermal(material_phase(ipc,ip,el)))
|
||||
case (LOCAL_THERMAL_adiabatic_ID)
|
||||
constitutive_getFi = math_mul33x33(constitutive_getFi,thermal_adiabatic_getFT (ipc, ip, el))
|
||||
|
||||
end select
|
||||
|
||||
end function constitutive_getFi
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the intermediate deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function constitutive_getFi0(ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal
|
||||
use math, only: &
|
||||
math_I3, &
|
||||
math_mul33x33
|
||||
use material, only: &
|
||||
phase_damage, &
|
||||
phase_thermal, &
|
||||
material_phase, &
|
||||
LOCAL_DAMAGE_anisoBrittle_ID, &
|
||||
LOCAL_THERMAL_adiabatic_ID
|
||||
use damage_anisoBrittle, only: &
|
||||
damage_anisoBrittle_getFd0
|
||||
use thermal_adiabatic, only: &
|
||||
thermal_adiabatic_getFT0
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), dimension(3,3) :: &
|
||||
constitutive_getFi0 !< intermediate deformation gradient
|
||||
|
||||
constitutive_getFi0 = math_I3
|
||||
|
||||
select case (phase_damage(material_phase(ipc,ip,el)))
|
||||
case (LOCAL_DAMAGE_anisoBrittle_ID)
|
||||
constitutive_getFi0 = math_mul33x33(constitutive_getFi0,damage_anisoBrittle_getFd0 (ipc, ip, el))
|
||||
|
||||
end select
|
||||
|
||||
select case (phase_thermal(material_phase(ipc,ip,el)))
|
||||
case (LOCAL_THERMAL_adiabatic_ID)
|
||||
constitutive_getFi0 = math_mul33x33(constitutive_getFi0,thermal_adiabatic_getFT0 (ipc, ip, el))
|
||||
|
||||
end select
|
||||
|
||||
end function constitutive_getFi0
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the intermediate deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function constitutive_getPartionedFi0(ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal
|
||||
use math, only: &
|
||||
math_I3, &
|
||||
math_mul33x33
|
||||
use material, only: &
|
||||
phase_damage, &
|
||||
phase_thermal, &
|
||||
material_phase, &
|
||||
LOCAL_DAMAGE_anisoBrittle_ID, &
|
||||
LOCAL_THERMAL_adiabatic_ID
|
||||
use damage_anisoBrittle, only: &
|
||||
damage_anisoBrittle_getPartionedFd0
|
||||
use thermal_adiabatic, only: &
|
||||
thermal_adiabatic_getPartionedFT0
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), dimension(3,3) :: &
|
||||
constitutive_getPartionedFi0 !< intermediate deformation gradient
|
||||
|
||||
constitutive_getPartionedFi0 = math_I3
|
||||
|
||||
select case (phase_damage(material_phase(ipc,ip,el)))
|
||||
case (LOCAL_DAMAGE_anisoBrittle_ID)
|
||||
constitutive_getPartionedFi0 = math_mul33x33(constitutive_getPartionedFi0, &
|
||||
damage_anisoBrittle_getPartionedFd0(ipc, ip, el))
|
||||
|
||||
end select
|
||||
|
||||
select case (phase_thermal(material_phase(ipc,ip,el)))
|
||||
case (LOCAL_THERMAL_adiabatic_ID)
|
||||
constitutive_getPartionedFi0 = math_mul33x33(constitutive_getPartionedFi0, &
|
||||
thermal_adiabatic_getPartionedFT0(ipc, ip, el))
|
||||
|
||||
end select
|
||||
|
||||
end function constitutive_getPartionedFi0
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
|
||||
|
|
@ -1160,41 +1355,6 @@ function constitutive_getSlipDamage(nSlip, Tstar_v, ipc, ip, el)
|
|||
|
||||
end function constitutive_getSlipDamage
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns damage deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function constitutive_getDamageStrain(ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal
|
||||
use math, only: &
|
||||
math_I3
|
||||
use material, only: &
|
||||
material_phase, &
|
||||
LOCAL_DAMAGE_anisoBrittle_ID, &
|
||||
phase_damage
|
||||
use damage_anisoBrittle, only: &
|
||||
damage_anisoBrittle_getDamageStrain
|
||||
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), dimension(3,3) :: &
|
||||
constitutive_getDamageStrain
|
||||
|
||||
|
||||
select case (phase_damage(material_phase(ipc,ip,el)))
|
||||
case (LOCAL_DAMAGE_anisoBrittle_ID)
|
||||
constitutive_getDamageStrain = damage_anisoBrittle_getDamageStrain(ipc, ip, el)
|
||||
case default
|
||||
constitutive_getDamageStrain = math_I3
|
||||
|
||||
end select
|
||||
|
||||
end function constitutive_getDamageStrain
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns damage diffusion tensor
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -1328,35 +1488,6 @@ function constitutive_getTemperature(ipc, ip, el)
|
|||
|
||||
end function constitutive_getTemperature
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns thermal deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function constitutive_getThermalStrain(ipc, ip, el)
|
||||
use prec, only: &
|
||||
pReal
|
||||
use math, only: &
|
||||
math_I3
|
||||
use lattice, only: &
|
||||
lattice_referenceTemperature, &
|
||||
lattice_thermalExpansion33
|
||||
use material, only: &
|
||||
mappingConstitutive
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), dimension(3,3) :: &
|
||||
constitutive_getThermalStrain
|
||||
|
||||
constitutive_getThermalStrain = math_I3 + &
|
||||
(constitutive_getTemperature(ipc, ip, el) - &
|
||||
lattice_referenceTemperature(mappingConstitutive(2,ipc,ip,el)))* &
|
||||
lattice_thermalExpansion33(1:3,1:3,mappingConstitutive(2,ipc,ip,el))
|
||||
|
||||
end function constitutive_getThermalStrain
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns local vacancy concentration
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -516,7 +516,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
|
|||
math_I3, &
|
||||
math_mul3333xx3333, &
|
||||
math_mul33xx33, &
|
||||
math_invert
|
||||
math_invert, &
|
||||
math_det33
|
||||
use FEsolving, only: &
|
||||
FEsolving_execElem, &
|
||||
FEsolving_execIP
|
||||
|
|
@ -537,7 +538,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
|
|||
homogenization_maxNgrains
|
||||
use constitutive, only: &
|
||||
constitutive_TandItsTangent, &
|
||||
constitutive_LpAndItsTangent
|
||||
constitutive_LpAndItsTangent, &
|
||||
constitutive_LiAndItsTangent, &
|
||||
constitutive_getFi, &
|
||||
constitutive_getFi0, &
|
||||
constitutive_getPartionedFi0
|
||||
|
||||
implicit none
|
||||
logical, intent(in) :: &
|
||||
|
|
@ -579,12 +584,21 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
|
|||
logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
|
||||
convergenceFlag_backup
|
||||
! local variables used for calculating analytic Jacobian
|
||||
real(pReal), dimension(3,3):: junk
|
||||
real(pReal) :: detInvFi
|
||||
real(pReal), dimension(3,3) :: junk, &
|
||||
Fi, &
|
||||
invFi, &
|
||||
Fi0, &
|
||||
invFi0
|
||||
real(pReal), dimension(3,3,3,3) :: dSdFe, &
|
||||
dFedF, &
|
||||
dSdF, &
|
||||
junk2, &
|
||||
dLpdS,dFpinvdF,rhs_3333,lhs_3333,temp_3333
|
||||
dLidS, &
|
||||
dLpdS, &
|
||||
dFpinvdF, &
|
||||
rhs_3333, &
|
||||
lhs_3333, &
|
||||
temp_3333
|
||||
real(pReal), dimension(9,9):: temp_99
|
||||
logical :: error
|
||||
|
||||
|
|
@ -624,8 +638,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
|
|||
crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness
|
||||
crystallite_subF0(1:3,1:3,g,i,e) = crystallite_partionedF0(1:3,1:3,g,i,e) ! ...def grad
|
||||
crystallite_subTstar0_v(1:6,g,i,e) = crystallite_partionedTstar0_v(1:6,g,i,e) !...2nd PK stress
|
||||
crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF0(1:3,1:3,g,i,e), &
|
||||
math_inv33(crystallite_subFp0(1:3,1:3,g,i,e))) ! only needed later on for stiffness calculation
|
||||
crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(crystallite_subF0(1:3,1:3,g,i,e), &
|
||||
math_inv33(crystallite_subFp0(1:3,1:3,g,i,e))), &
|
||||
math_inv33(constitutive_getFi0(g,i,e))) ! only needed later on for stiffness calculation
|
||||
crystallite_subFrac(g,i,e) = 0.0_pReal
|
||||
crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst
|
||||
crystallite_todo(g,i,e) = .true.
|
||||
|
|
@ -893,8 +908,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
|
|||
!$OMP FLUSH(crystallite_subF0)
|
||||
crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient
|
||||
crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad
|
||||
crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
|
||||
crystallite_invFp(1:3,1:3,g,i,e)) ! only needed later on for stiffness calculation
|
||||
crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(crystallite_subF (1:3,1:3,g,i,e), &
|
||||
crystallite_invFp(1:3,1:3,g,i,e)), &
|
||||
math_inv33(constitutive_getFi(g,i,e))) ! only needed later on for stiffness calculation
|
||||
!if abbrevation, make c and p private in omp
|
||||
plasticState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e)) = &
|
||||
plasticState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e))
|
||||
|
|
@ -979,6 +995,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
|
|||
- crystallite_partionedF0(1:3,1:3,g,i,e))
|
||||
!$OMP FLUSH(crystallite_subF)
|
||||
crystallite_Fe(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e))
|
||||
crystallite_Fe(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(crystallite_subF (1:3,1:3,g,i,e), &
|
||||
crystallite_invFp(1:3,1:3,g,i,e)), &
|
||||
math_inv33(constitutive_getFi(g,i,e)))
|
||||
crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e)
|
||||
crystallite_converged(g,i,e) = .false. ! start out non-converged
|
||||
endif
|
||||
|
|
@ -1054,9 +1073,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
|
|||
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> no convergence: respond fully elastic at el (elFE) ip g ', &
|
||||
e,'(',mesh_element(1,e),')',i,g
|
||||
invFp = math_inv33(crystallite_partionedFp0(1:3,1:3,g,i,e))
|
||||
Fe_guess = math_mul33x33(crystallite_partionedF(1:3,1:3,g,i,e), invFp)
|
||||
Fe_guess = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,g,i,e), invFp), &
|
||||
math_inv33(constitutive_getPartionedFi0(g,i,e)))
|
||||
call constitutive_TandItsTangent(Tstar, junk2, Fe_guess,g,i,e)
|
||||
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp)))
|
||||
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(Fe_guess,constitutive_getPartionedFi0(g,i,e)), &
|
||||
math_mul33x33(Tstar,transpose(invFp)))* &
|
||||
math_det33(constitutive_getPartionedFi0(g,i,e))
|
||||
endif
|
||||
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
|
||||
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
|
||||
|
|
@ -1082,7 +1104,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
|
|||
|
||||
! --- ANALYTIC JACOBIAN ---
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(dFedF,dSdF,dSdFe,dLpdS,dFpinvdF,rhs_3333,lhs_3333,temp_3333,temp_99,junk,myNgrains)
|
||||
!$OMP PARALLEL DO PRIVATE(dFedF,dSdF,dSdFe,dLpdS,dFpinvdF,rhs_3333,lhs_3333,temp_99,junk,dLidS,Fi,invFi,Fi0,invFi0,detInvFi,temp_3333,myNgrains)
|
||||
elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2)
|
||||
myNgrains = homogenization_Ngrains(mesh_element(3,e))
|
||||
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
|
||||
|
|
@ -1091,46 +1113,71 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
|
|||
call constitutive_LpAndItsTangent(junk,temp_99,crystallite_Tstar_v(1:6,g,i,e), &
|
||||
g,i,e)
|
||||
dLpdS = reshape(temp_99,shape=[3,3,3,3])
|
||||
call constitutive_LiAndItsTangent(junk,temp_99,crystallite_Tstar_v(1:6,g,i,e), &
|
||||
crystallite_Lp(1:3,1:3,g,i,e), g,i,e)
|
||||
dLidS = reshape(temp_99,shape=[3,3,3,3])
|
||||
Fi = constitutive_getFi(g,i,e)
|
||||
invFi = math_inv33(Fi)
|
||||
Fi0 = constitutive_getFi0(g,i,e)
|
||||
invFi0 = math_inv33(Fi0)
|
||||
detInvFi = math_det33(invFi)
|
||||
|
||||
junk = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e),invFi))
|
||||
rhs_3333 = 0.0_pReal
|
||||
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
||||
rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3), &
|
||||
math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))
|
||||
junk = math_mul33x33(math_mul33x33(crystallite_Fe(1:3,1:3,g,i,e), &
|
||||
crystallite_Fp(1:3,1:3,g,i,e)), &
|
||||
math_inv33(crystallite_Fp0(1:3,1:3,g,i,e)))
|
||||
rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),junk)
|
||||
|
||||
temp_3333 = 0.0_pReal
|
||||
junk = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
|
||||
math_inv33(crystallite_subFp0(1:3,1:3,g,i,e)))
|
||||
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
||||
temp_3333(1:3,1:3,p,o) = math_mul33x33(junk,dLpdS(1:3,1:3,p,o))
|
||||
lhs_3333 = crystallite_dt(g,i,e)*math_mul3333xx3333(dSdFe,temp_3333)
|
||||
temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(junk,dLpdS(1:3,1:3,p,o)),invFi)
|
||||
|
||||
junk = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
|
||||
crystallite_invFp(1:3,1:3,g,i,e)), invFi0)
|
||||
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
||||
temp_3333(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + math_mul33x33(junk,dLidS(1:3,1:3,p,o))
|
||||
|
||||
lhs_3333 = crystallite_subdt(g,i,e)*math_mul3333xx3333(dSdFe,temp_3333)
|
||||
|
||||
call math_invert(9_pInt,math_identity2nd(9_pInt)+reshape(lhs_3333,shape=[9,9]),temp_99,error)
|
||||
if (error) call IO_error(error_ID=400_pInt,ext_msg='analytic tangent inversion')
|
||||
dSdF = math_mul3333xx3333(reshape(temp_99,shape=[3,3,3,3]),rhs_3333)
|
||||
temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
|
||||
|
||||
dFpinvdF = 0.0_pReal
|
||||
temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
|
||||
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
||||
dFpinvdF(1:3,1:3,p,o) = -crystallite_dt(g,i,e)* &
|
||||
math_mul33x33(math_inv33(crystallite_Fp0(1:3,1:3,g,i,e)), &
|
||||
temp_3333(1:3,1:3,p,o))
|
||||
temp_3333 = 0.0_pReal
|
||||
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
||||
temp_3333(o,p,o,1:3) = crystallite_invFp(1:3,p,g,i,e) ! dFe^T_ij/dF_kl = delta_jk * (Fp current^-1)_li
|
||||
dFedF = 0.0_pReal
|
||||
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
||||
dFedF(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + &
|
||||
math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
|
||||
dFpinvdF(1:3,1:3,p,o))
|
||||
dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(g,i,e)* &
|
||||
math_mul33x33(math_inv33(crystallite_subFp0(1:3,1:3,g,i,e)), &
|
||||
math_mul33x33(temp_3333(1:3,1:3,p,o),invFi))
|
||||
|
||||
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.0_pReal
|
||||
junk = math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e), &
|
||||
math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), &
|
||||
math_transpose33(crystallite_invFp(1:3,1:3,g,i,e))))/detInvFi
|
||||
forall(p=1_pInt:3_pInt) &
|
||||
crystallite_dPdF(p,1:3,p,1:3,g,i,e) = math_transpose33(junk)
|
||||
|
||||
junk = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), &
|
||||
math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))/detInvFi
|
||||
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
||||
crystallite_dPdF(1:3,1:3,o,p,g,i,e) = &
|
||||
math_mul33x33(math_mul33x33(dFedF(1:3,1:3,o,p),&
|
||||
math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e))), &
|
||||
math_transpose33(crystallite_invFp(1:3,1:3,g,i,e))) + & ! dP/dF = dFe/dF * S * Fp^-T...
|
||||
math_mul33x33(crystallite_Fe(1:3,1:3,g,i,e), &
|
||||
math_mul33x33(dSdF(1:3,1:3,o,p), &
|
||||
math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))) + & ! + Fe * dS/dF * Fp^-T
|
||||
math_mul33x33(crystallite_Fe(1:3,1:3,g,i,e), &
|
||||
math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), &
|
||||
math_transpose33(dFpinvdF(1:3,1:3,p,o)))) ! + Fe * S * dFp^-T/dF
|
||||
crystallite_dPdF(1:3,1:3,p,o,g,i,e) = crystallite_dPdF(1:3,1:3,p,o,g,i,e) + &
|
||||
math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e),dFpinvdF(1:3,1:3,p,o)),junk)
|
||||
|
||||
junk = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
|
||||
crystallite_invFp(1:3,1:3,g,i,e))/detInvFi
|
||||
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
||||
crystallite_dPdF(1:3,1:3,p,o,g,i,e) = crystallite_dPdF(1:3,1:3,p,o,g,i,e) + &
|
||||
math_mul33x33(math_mul33x33(junk,dSdF(1:3,1:3,p,o)), &
|
||||
math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))
|
||||
|
||||
junk = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
|
||||
crystallite_invFp(1:3,1:3,g,i,e)), &
|
||||
math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)))/detInvFi
|
||||
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
||||
crystallite_dPdF(1:3,1:3,p,o,g,i,e) = crystallite_dPdF(1:3,1:3,p,o,g,i,e) + &
|
||||
math_mul33x33(junk,math_transpose33(dFpinvdF(1:3,1:3,p,o)))
|
||||
|
||||
enddo; enddo
|
||||
enddo elementLooping6
|
||||
!$OMP END PARALLEL DO
|
||||
|
|
@ -3296,11 +3343,12 @@ logical function crystallite_integrateStress(&
|
|||
debug_cumLpTicks, &
|
||||
debug_StressLoopDistribution
|
||||
use constitutive, only: constitutive_LpAndItsTangent, &
|
||||
constitutive_TandItsTangent, &
|
||||
constitutive_getThermalStrain, &
|
||||
constitutive_getDamageStrain
|
||||
constitutive_LiAndItsTangent, &
|
||||
constitutive_getFi0, &
|
||||
constitutive_TandItsTangent
|
||||
use math, only: math_mul33x33, &
|
||||
math_mul33xx33, &
|
||||
math_mul3333xx3333, &
|
||||
math_mul66x6, &
|
||||
math_mul99x99, &
|
||||
math_transpose33, &
|
||||
|
|
@ -3328,45 +3376,65 @@ logical function crystallite_integrateStress(&
|
|||
!*** local variables ***!
|
||||
real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep
|
||||
Fp_current, & ! plastic deformation gradient at start of timestep
|
||||
Fi_current, & ! intermediate deformation gradient at start of timestep
|
||||
Fp_new, & ! plastic deformation gradient at end of timestep
|
||||
Fe_new, & ! elastic deformation gradient at end of timestep
|
||||
invFp_new, & ! inverse of Fp_new
|
||||
invFp_current, & ! inverse of Fp_current
|
||||
invFi_current, & ! inverse of Fp_current
|
||||
Lpguess, & ! current guess for plastic velocity gradient
|
||||
Lpguess_old, & ! known last good guess for plastic velocity gradient
|
||||
Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law
|
||||
residuum, & ! current residuum of plastic velocity gradient
|
||||
residuum_old, & ! last residuum of plastic velocity gradient
|
||||
deltaLp, & ! direction of next guess
|
||||
Tstar,& ! 2nd Piola-Kirchhoff Stress
|
||||
A,&
|
||||
Liguess, & ! current guess for plastic velocity gradient
|
||||
Liguess_old, & ! known last good guess for plastic velocity gradient
|
||||
Li_constitutive, & ! plastic velocity gradient resulting from constitutive law
|
||||
residuumI, & ! current residuum of plastic velocity gradient
|
||||
residuumI_old, & ! last residuum of plastic velocity gradient
|
||||
deltaLi, & ! direction of next guess
|
||||
Tstar, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration
|
||||
Tstar_unloaded, & ! 2nd Piola-Kirchhoff Stress in unloaded configuration
|
||||
A, &
|
||||
B, &
|
||||
Fe, & ! elastic deformation gradient
|
||||
Fi, & ! gradient of intermediate deformation stages
|
||||
Ci, & ! stretch of intermediate deformation stages
|
||||
invFi
|
||||
invFi, &
|
||||
temp_33
|
||||
real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation
|
||||
real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK
|
||||
integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK
|
||||
real(pReal), dimension(9,9) :: dLp_dT_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law
|
||||
dLi_dT_constitutive, & ! partial derivative of intermediate velocity gradient calculated by constitutive law
|
||||
dT_dFe_constitutive, & ! partial derivative of 2nd Piola-Kirchhoff stress calculated by constitutive law
|
||||
dFe_dLp, & ! partial derivative of elastic deformation gradient
|
||||
dR_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme)
|
||||
dR_dLp2 ! working copy of dRdLp
|
||||
dR_dLp2, & ! working copy of dRdLp
|
||||
dRI_dLp, & ! partial derivative of residuumI (Jacobian for NEwton-Raphson scheme)
|
||||
dRI_dLp2 ! working copy of dRIdLp
|
||||
real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress
|
||||
dFe_dLp3333 ! partial derivative of elastic deformation gradient
|
||||
dFe_dLp3333, & ! partial derivative of elastic deformation gradient
|
||||
dInvFi_dLp3333, &
|
||||
dFe_dInvFi3333, &
|
||||
dT_dInvFi3333, &
|
||||
temp_3333
|
||||
real(pReal) det, & ! determinant
|
||||
detInvFi, &
|
||||
steplength0, &
|
||||
steplength, &
|
||||
steplengthI0, &
|
||||
steplengthI, &
|
||||
dt, & ! time increment
|
||||
aTol, &
|
||||
detFi
|
||||
aTol
|
||||
logical error ! flag indicating an error
|
||||
integer(pInt) NiterationStress, & ! number of stress integrations
|
||||
NiterationStressI, & ! number of inner stress integrations
|
||||
ierr, & ! error indicator for LAPACK
|
||||
o, &
|
||||
p, &
|
||||
jacoCounter ! counter to check for Jacobian update
|
||||
jacoCounter, &
|
||||
jacoICounter ! counters to check for Jacobian update
|
||||
integer(pLongInt) tick, &
|
||||
tock, &
|
||||
tickrate, &
|
||||
|
|
@ -3424,20 +3492,28 @@ logical function crystallite_integrateStress(&
|
|||
return
|
||||
endif
|
||||
A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
|
||||
Fi = constitutive_getDamageStrain(g,i,e) ! intermediate configuration, assume decomposition as F = Fe Fi Fp
|
||||
Fi = math_mul33x33(Fi,constitutive_getThermalStrain(g,i,e)) ! Fi = Ft Fd
|
||||
Ci = math_mul33x33(math_transpose33(Fi),Fi) ! non-plastic stretch tensor (neglecting elastic contributions)
|
||||
invFi = math_inv33(Fi)
|
||||
detFi = math_det33(Fi)
|
||||
|
||||
Fi_current = constitutive_getFi0(g,i,e) ! intermediate configuration, assume decomposition as F = Fe Fi Fp
|
||||
invFi_current = math_inv33(Fi_current)
|
||||
if (all(invFi_current == 0.0_pReal)) then ! ... failed?
|
||||
#ifndef _OPENMP
|
||||
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
|
||||
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip g ',&
|
||||
e,'(',mesh_element(1,e),')',i,g
|
||||
if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0_pInt) &
|
||||
write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp_current',math_transpose33(Fi_current(1:3,1:3))
|
||||
endif
|
||||
#endif
|
||||
return
|
||||
endif
|
||||
|
||||
!* start LpLoop with normal step length
|
||||
|
||||
NiterationStress = 0_pInt
|
||||
jacoCounter = 0_pInt
|
||||
steplength0 = 1.0_pReal
|
||||
steplength = steplength0
|
||||
residuum_old = 0.0_pReal
|
||||
jacoCounter = 0_pInt
|
||||
steplength0 = 1.0_pReal
|
||||
steplength = steplength0
|
||||
residuum_old = 0.0_pReal
|
||||
|
||||
LpLoop: do
|
||||
NiterationStress = NiterationStress + 1_pInt
|
||||
|
|
@ -3449,17 +3525,106 @@ logical function crystallite_integrateStress(&
|
|||
#endif
|
||||
return
|
||||
endif loopsExeced
|
||||
B = math_I3 - dt*Lpguess
|
||||
NiterationStressI = 0_pInt
|
||||
Liguess = 0.0_pReal
|
||||
Liguess_old = Liguess
|
||||
jacoICounter = 0_pInt
|
||||
steplengthI0 = 1.0_pReal
|
||||
steplengthI = steplengthI0
|
||||
residuumI_old = 0.0_pReal
|
||||
LiLoop: do ! inner stress integration loop for consistency with Fi
|
||||
NiterationStressI = NiterationStressI + 1_pInt
|
||||
IloopsExeced: if (NiterationStressI > nStress) then
|
||||
#ifndef _OPENMP
|
||||
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
|
||||
write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached inelastic loop limit',nStress, &
|
||||
' at el (elFE) ip g ', e,mesh_element(1,e),i,g
|
||||
#endif
|
||||
return
|
||||
endif IloopsExeced
|
||||
|
||||
!* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
|
||||
|
||||
invFi = math_mul33x33(invFi_current,math_I3 - dt*Liguess)
|
||||
detInvFi = math_det33(invFi)
|
||||
Fi = math_inv33(invFi)
|
||||
Fe = math_mul33x33(math_mul33x33(A,B),invFi) ! current elastic deformation tensor
|
||||
call constitutive_TandItsTangent(Tstar_unloaded, dT_dFe3333, Fe, g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration
|
||||
Tstar = math_mul33x33(invFi, &
|
||||
math_mul33x33(Tstar_unloaded,math_transpose33(invFi)))/detInvFi ! push Tstar forward from unloaded to plastic (lattice) configuration
|
||||
Tstar_v = math_Mandel33to6(Tstar)
|
||||
|
||||
!* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
|
||||
!* calculate intermediate velocity gradient and its tangent from constitutive law
|
||||
|
||||
B = math_I3 - dt*Lpguess
|
||||
Fe = math_mul33x33(math_mul33x33(A,B),invFi) ! current elastic deformation tensor
|
||||
call constitutive_TandItsTangent(Tstar, dT_dFe3333, Fe, g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration
|
||||
Tstar = math_mul33x33(Ci, math_mul33x33(invFi, &
|
||||
math_mul33x33(Tstar,math_transpose33(invFi))))/detFi ! push Tstar forward to plastic (lattice) configuration
|
||||
Tstar_v = math_Mandel33to6(Tstar)
|
||||
call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT_constitutive, Tstar_v, Lpguess, &
|
||||
g, i, e)
|
||||
|
||||
!* update current residuum and check for convergence of loop
|
||||
|
||||
aTol = max(rTol_crystalliteStress * max(math_norm33(Liguess),math_norm33(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
|
||||
aTol_crystalliteStress) ! minimum lower cutoff
|
||||
residuumI = Liguess - Li_constitutive
|
||||
if (any(residuumI /= residuumI)) then ! NaN in residuum...
|
||||
return ! ...me = .false. to inform integrator about problem
|
||||
elseif (math_norm33(residuumI) < aTol) then ! converged if below absolute tolerance
|
||||
exit LiLoop ! ...leave iteration loop
|
||||
elseif (math_norm33(residuumI) < math_norm33(residuumI_old) .or. NiterationStress == 1_pInt ) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
|
||||
residuumI_old = residuumI ! ...remember old values and...
|
||||
Liguess_old = Liguess
|
||||
steplengthI = steplengthI0 ! ...proceed with normal step length (calculate new search direction)
|
||||
else ! not converged and residuum not improved...
|
||||
steplengthI = 0.5_pReal * steplengthI ! ...try with smaller step length in same direction
|
||||
Liguess = Liguess_old + steplengthI * deltaLi
|
||||
cycle LiLoop
|
||||
endif
|
||||
|
||||
!* calculate Jacobian for correction term
|
||||
|
||||
if (mod(jacoICounter, iJacoLpresiduum) == 0_pInt) then
|
||||
dInvFi_dLp3333 = 0.0_pReal
|
||||
do o=1_pInt,3_pInt
|
||||
dInvFi_dLp3333(1:3,o,1:3,o) = invFi_current
|
||||
enddo
|
||||
dInvFi_dLp3333 = -dt * dInvFi_dLp3333
|
||||
dFe_dInvFi3333 = 0.0_pReal
|
||||
temp_33 = math_mul33x33(A,B)
|
||||
do o=1_pInt,3_pInt
|
||||
dFe_dInvFi3333(1:3,o,1:3,o) = temp_33
|
||||
enddo
|
||||
dT_dInvFi3333 = 0.0_pReal
|
||||
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
|
||||
dT_dInvFi3333(1:3,1:3,p,o) = -Tstar * Fi(o,p)
|
||||
enddo; enddo
|
||||
temp_33 = math_mul33x33(invFi,Tstar_unloaded)/detInvFi
|
||||
do o=1_pInt,3_pInt
|
||||
dT_dInvFi3333(1:3,o,o,1:3) = dT_dInvFi3333(1:3,o,o,1:3) + temp_33
|
||||
dT_dInvFi3333(o,1:3,o,1:3) = dT_dInvFi3333(o,1:3,o,1:3) + temp_33
|
||||
enddo
|
||||
temp_3333 = math_mul3333xx3333(dT_dFe3333,dFe_dInvFi3333)
|
||||
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
|
||||
dT_dInvFi3333(1:3,1:3,o,p) = dT_dInvFi3333(1:3,1:3,o,p) + &
|
||||
math_mul33x33(math_mul33x33(invFi,temp_3333(1:3,1:3,o,p)), &
|
||||
math_transpose33(invFi))/detInvFi
|
||||
enddo; enddo
|
||||
dRI_dLp = math_identity2nd(9_pInt) - &
|
||||
math_mul99x99(dLi_dT_constitutive,math_Plain3333to99(math_mul3333xx3333(dT_dInvFi3333,dInvFi_dLp3333)))
|
||||
dRI_dLp2 = dRI_dLp ! will be overwritten in first call to LAPACK routine
|
||||
work = math_plain33to9(residuumI)
|
||||
#if(FLOAT==8)
|
||||
call dgesv(9,1,dRI_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp
|
||||
#elif(FLOAT==4)
|
||||
call sgesv(9,1,dRI_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp
|
||||
#endif
|
||||
if (ierr /= 0_pInt) then
|
||||
return
|
||||
endif
|
||||
deltaLi = - math_plain9to33(work)
|
||||
endif
|
||||
jacoICounter = jacoICounter + 1_pInt ! increase counter for jaco update
|
||||
|
||||
Liguess = Liguess + steplengthI * deltaLi
|
||||
enddo LiLoop
|
||||
|
||||
!* calculate plastic velocity gradient and its tangent from constitutive law
|
||||
|
||||
|
|
@ -3524,7 +3689,7 @@ logical function crystallite_integrateStress(&
|
|||
if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then
|
||||
dFe_dLp3333 = 0.0_pReal
|
||||
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
|
||||
dFe_dLp3333(o,p,1:3,p) = A(o,1:3) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) delta(l,j)
|
||||
dFe_dLp3333(o,1:3,p,1:3) = A(o,p)*math_transpose33(invFi) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
|
||||
enddo; enddo
|
||||
dFe_dLp3333 = -dt * dFe_dLp3333
|
||||
dFe_dLp = math_Plain3333to99(dFe_dLp3333)
|
||||
|
|
@ -3594,7 +3759,7 @@ logical function crystallite_integrateStress(&
|
|||
|
||||
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(Fe_new,Fi), &
|
||||
math_mul33x33(math_Mandel6to33(Tstar_v), &
|
||||
math_transpose33(invFp_new)))*detFi
|
||||
math_transpose33(invFp_new)))/detInvFi
|
||||
|
||||
|
||||
!* store local values in global variables
|
||||
|
|
|
|||
|
|
@ -32,7 +32,8 @@ module damage_anisoBrittle
|
|||
damage_anisoBrittle_Ncleavage !< number of cleavage systems per family
|
||||
|
||||
real(pReal), dimension(:), allocatable, private :: &
|
||||
damage_anisoBrittle_aTol, &
|
||||
damage_anisoBrittle_aTol_damage, &
|
||||
damage_anisoBrittle_aTol_disp, &
|
||||
damage_anisoBrittle_sdot_0, &
|
||||
damage_anisoBrittle_N
|
||||
|
||||
|
|
@ -54,10 +55,13 @@ module damage_anisoBrittle
|
|||
damage_anisoBrittle_stateInit, &
|
||||
damage_anisoBrittle_aTolState, &
|
||||
damage_anisoBrittle_dotState, &
|
||||
damage_anisoBrittle_LdAndItsTangent, &
|
||||
damage_anisoBrittle_getFd, &
|
||||
damage_anisoBrittle_getFd0, &
|
||||
damage_anisoBrittle_getPartionedFd0, &
|
||||
damage_anisoBrittle_getDamage, &
|
||||
damage_anisoBrittle_putLocalDamage, &
|
||||
damage_anisoBrittle_getLocalDamage, &
|
||||
damage_anisoBrittle_getDamageStrain, &
|
||||
damage_anisoBrittle_postResults
|
||||
|
||||
contains
|
||||
|
|
@ -139,9 +143,10 @@ subroutine damage_anisoBrittle_init(fileUnit)
|
|||
allocate(damage_anisoBrittle_Noutput(maxNinstance), source=0_pInt)
|
||||
allocate(damage_anisoBrittle_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
|
||||
allocate(damage_anisoBrittle_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
|
||||
allocate(damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt)
|
||||
allocate(damage_anisoBrittle_totalNcleavage(maxNinstance), source=0_pInt)
|
||||
allocate(damage_anisoBrittle_aTol(maxNinstance), source=0.0_pReal)
|
||||
allocate(damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,maxNinstance),source=0_pInt)
|
||||
allocate(damage_anisoBrittle_totalNcleavage(maxNinstance), source=0_pInt)
|
||||
allocate(damage_anisoBrittle_aTol_damage(maxNinstance), source=0.0_pReal)
|
||||
allocate(damage_anisoBrittle_aTol_disp(maxNinstance), source=0.0_pReal)
|
||||
allocate(damage_anisoBrittle_sdot_0(maxNinstance), source=0.0_pReal)
|
||||
allocate(damage_anisoBrittle_N(maxNinstance), source=0.0_pReal)
|
||||
|
||||
|
|
@ -177,7 +182,10 @@ subroutine damage_anisoBrittle_init(fileUnit)
|
|||
end select
|
||||
|
||||
case ('atol_damage')
|
||||
damage_anisoBrittle_aTol(instance) = IO_floatValue(line,positions,2_pInt)
|
||||
damage_anisoBrittle_aTol_damage(instance) = IO_floatValue(line,positions,2_pInt)
|
||||
|
||||
case ('atol_disp')
|
||||
damage_anisoBrittle_aTol_disp(instance) = IO_floatValue(line,positions,2_pInt)
|
||||
|
||||
case ('sdot_0')
|
||||
damage_anisoBrittle_sdot_0(instance) = IO_floatValue(line,positions,2_pInt)
|
||||
|
|
@ -226,8 +234,8 @@ subroutine damage_anisoBrittle_init(fileUnit)
|
|||
enddo outputsLoop
|
||||
! Determine size of state array
|
||||
sizeDotState = 2_pInt + & ! viscous and non-viscous damage values
|
||||
9_pInt + & ! damage deformation gradient
|
||||
3_pInt*damage_anisoBrittle_totalNcleavage(instance) ! opening on each damage system
|
||||
3_pInt*damage_anisoBrittle_totalNcleavage(instance) + & ! opening on each damage system
|
||||
9_pInt ! Fd
|
||||
sizeState = sizeDotState
|
||||
|
||||
damageState(phase)%sizeState = sizeState
|
||||
|
|
@ -252,7 +260,7 @@ subroutine damage_anisoBrittle_init(fileUnit)
|
|||
if (any(numerics_integrator == 5_pInt)) &
|
||||
allocate(damageState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
|
||||
|
||||
call damage_anisoBrittle_stateInit(phase)
|
||||
call damage_anisoBrittle_stateInit(phase,instance)
|
||||
call damage_anisoBrittle_aTolState(phase,instance)
|
||||
endif
|
||||
|
||||
|
|
@ -262,27 +270,22 @@ end subroutine damage_anisoBrittle_init
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief sets the relevant state values for a given instance of this damage
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine damage_anisoBrittle_stateInit(phase)
|
||||
subroutine damage_anisoBrittle_stateInit(phase,instance)
|
||||
use material, only: &
|
||||
damageState
|
||||
use math, only: &
|
||||
math_I3
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: phase !< number specifying the phase of the damage
|
||||
integer(pInt), intent(in) :: phase, instance !< number specifying the phase of the damage
|
||||
|
||||
real(pReal), dimension(damageState(phase)%sizeState) :: tempState
|
||||
|
||||
tempState = 0.0_pReal
|
||||
tempState(1) = 1.0_pReal
|
||||
tempState(2) = 1.0_pReal
|
||||
tempState(3) = 1.0_pReal
|
||||
tempState(4) = 0.0_pReal
|
||||
tempState(5) = 0.0_pReal
|
||||
tempState(6) = 0.0_pReal
|
||||
tempState(7) = 1.0_pReal
|
||||
tempState(8) = 0.0_pReal
|
||||
tempState(9) = 0.0_pReal
|
||||
tempState(10) = 0.0_pReal
|
||||
tempState(11) = 1.0_pReal
|
||||
tempState(3:2+3*damage_anisoBrittle_totalNcleavage(instance)) = 0.0_pReal
|
||||
tempState(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
|
||||
3*damage_anisoBrittle_totalNcleavage(instance)+11) = reshape(math_I3, shape=[9])
|
||||
damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:)))
|
||||
damageState(phase)%state0 = damageState(phase)%state
|
||||
damageState(phase)%partionedState0 = damageState(phase)%state
|
||||
|
|
@ -301,7 +304,11 @@ subroutine damage_anisoBrittle_aTolState(phase,instance)
|
|||
instance ! number specifying the current instance of the damage
|
||||
real(pReal), dimension(damageState(phase)%sizeState) :: tempTol
|
||||
|
||||
tempTol = damage_anisoBrittle_aTol(instance)
|
||||
tempTol(1) = damage_anisoBrittle_aTol_damage(instance)
|
||||
tempTol(2) = damage_anisoBrittle_aTol_damage(instance)
|
||||
tempTol(3:2+3*damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_disp(instance)
|
||||
tempTol(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
|
||||
3*damage_anisoBrittle_totalNcleavage(instance)+11) = damage_anisoBrittle_aTol_disp(instance)
|
||||
damageState(phase)%aTolState = tempTol
|
||||
end subroutine damage_anisoBrittle_aTolState
|
||||
|
||||
|
|
@ -330,13 +337,13 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el)
|
|||
el !< element
|
||||
real(pReal), intent(in), dimension(6) :: &
|
||||
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
|
||||
real(pReal), dimension(3,3) :: &
|
||||
Ld, Fd
|
||||
integer(pInt) :: &
|
||||
phase, &
|
||||
constituent, &
|
||||
instance, &
|
||||
f, i, index, index_myFamily
|
||||
real(pReal), dimension(3,3) :: &
|
||||
Ld, Fd
|
||||
real(pReal) :: &
|
||||
traction_d, traction_t, traction_n, traction_crit, &
|
||||
opening, &
|
||||
|
|
@ -350,20 +357,20 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el)
|
|||
(1.0_pReal/lattice_DamageMobility(phase))* &
|
||||
(damageState(phase)%state(2,constituent) - &
|
||||
damageState(phase)%state(1,constituent))
|
||||
|
||||
damageState(phase)%dotState(2,constituent) = 0.0_pReal
|
||||
nonLocalFactor = 1.0_pReal + &
|
||||
(damage_anisoBrittle_getDamage(ipc, ip, el) - &
|
||||
damageState(phase)%state(1,constituent))
|
||||
Ld = 0.0_pReal
|
||||
index = 11_pInt
|
||||
(damage_anisoBrittle_getDamage (ipc, ip, el) - &
|
||||
damage_anisoBrittle_getLocalDamage(ipc, ip, el))
|
||||
Ld = 0.0_pReal
|
||||
index = 2_pInt
|
||||
do f = 1_pInt,lattice_maxNcleavageFamily
|
||||
index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family
|
||||
do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
|
||||
opening = math_norm3(damageState(phase)%state(index+1:index+3,constituent))
|
||||
opening = math_norm3(damageState(phase)%state0(index+1:index+3,constituent))
|
||||
traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
|
||||
traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
|
||||
traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
|
||||
traction_crit = (1.0_pReal - opening/damage_anisoBrittle_critDisp(f,instance))* &
|
||||
traction_crit = max(0.001_pReal,(1.0_pReal - opening/damage_anisoBrittle_critDisp(f,instance)))* &
|
||||
damage_anisoBrittle_critLoad(f,instance)*nonLocalFactor
|
||||
damageState(phase)%dotState(index+1,constituent) = &
|
||||
sign(1.0_pReal,traction_d)* &
|
||||
|
|
@ -393,11 +400,218 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el)
|
|||
index = index + 3_pInt
|
||||
enddo
|
||||
enddo
|
||||
Fd = reshape(damageState(phase)%state(3:11,constituent),shape=[3,3])
|
||||
damageState(phase)%dotState(3:11,constituent) = reshape(math_mul33x33(Ld,Fd),shape=[9])
|
||||
Fd = reshape(damageState(phase)%state(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
|
||||
3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), &
|
||||
shape=[3,3])
|
||||
damageState(phase)%dotState(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
|
||||
3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent) = &
|
||||
reshape(math_mul33x33(Ld,Fd), shape=[9])
|
||||
|
||||
end subroutine damage_anisoBrittle_dotState
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the velocity gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el)
|
||||
use material, only: &
|
||||
mappingConstitutive, &
|
||||
phase_damageInstance, &
|
||||
damageState
|
||||
use math, only: &
|
||||
math_norm3, &
|
||||
math_Plain3333to99
|
||||
use lattice, only: &
|
||||
lattice_Scleavage, &
|
||||
lattice_Scleavage_v, &
|
||||
lattice_maxNcleavageFamily, &
|
||||
lattice_NcleavageSystem
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), intent(in), dimension(6) :: &
|
||||
Tstar_v !< 2nd Piola-Kirchhoff stress
|
||||
real(pReal), intent(out), dimension(3,3) :: &
|
||||
Ld !< damage velocity gradient
|
||||
real(pReal), intent(out), dimension(9,9) :: &
|
||||
dLd_dTstar !< derivative of Ld with respect to Tstar (2nd-order tensor)
|
||||
real(pReal), dimension(3,3,3,3) :: &
|
||||
dLd_dTstar3333 !< derivative of Ld with respect to Tstar (4th-order tensor)
|
||||
integer(pInt) :: &
|
||||
phase, &
|
||||
constituent, &
|
||||
instance, &
|
||||
f, i, index, index_myFamily, k, l, m, n
|
||||
real(pReal) :: &
|
||||
traction_d, traction_t, traction_n, traction_crit, &
|
||||
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt, &
|
||||
opening, &
|
||||
nonLocalFactor
|
||||
|
||||
phase = mappingConstitutive(2,ipc,ip,el)
|
||||
constituent = mappingConstitutive(1,ipc,ip,el)
|
||||
instance = phase_damageInstance(phase)
|
||||
|
||||
Ld = 0.0_pReal
|
||||
dLd_dTstar3333 = 0.0_pReal
|
||||
nonLocalFactor = 1.0_pReal + &
|
||||
(damage_anisoBrittle_getDamage (ipc, ip, el) - &
|
||||
damage_anisoBrittle_getLocalDamage(ipc, ip, el))
|
||||
index = 2_pInt
|
||||
do f = 1_pInt,lattice_maxNcleavageFamily
|
||||
index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family
|
||||
do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
|
||||
opening = math_norm3(damageState(phase)%state0(index+1:index+3,constituent))
|
||||
traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
|
||||
traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
|
||||
traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
|
||||
traction_crit = max(0.001_pReal,(1.0_pReal - opening/damage_anisoBrittle_critDisp(f,instance)))* &
|
||||
damage_anisoBrittle_critLoad(f,instance)*nonLocalFactor
|
||||
udotd = &
|
||||
sign(1.0_pReal,traction_d)* &
|
||||
damage_anisoBrittle_sdot_0(instance)* &
|
||||
(abs(traction_d)/traction_crit)**damage_anisoBrittle_N(instance)
|
||||
if (udotd /= 0.0_pReal) then
|
||||
Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)
|
||||
dudotd_dt = udotd*damage_anisoBrittle_N(instance)/traction_d
|
||||
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
|
||||
dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
|
||||
dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* &
|
||||
lattice_Scleavage(m,n,1,index_myFamily+i,phase)
|
||||
endif
|
||||
|
||||
udott = &
|
||||
sign(1.0_pReal,traction_t)* &
|
||||
damage_anisoBrittle_sdot_0(instance)* &
|
||||
(abs(traction_t)/traction_crit)**damage_anisoBrittle_N(instance)
|
||||
if (udott /= 0.0_pReal) then
|
||||
Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)
|
||||
dudott_dt = udott*damage_anisoBrittle_N(instance)/traction_t
|
||||
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
|
||||
dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
|
||||
dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* &
|
||||
lattice_Scleavage(m,n,2,index_myFamily+i,phase)
|
||||
endif
|
||||
|
||||
udotn = &
|
||||
sign(1.0_pReal,traction_n)* &
|
||||
damage_anisoBrittle_sdot_0(instance)* &
|
||||
(abs(traction_n)/traction_crit)**damage_anisoBrittle_N(instance)
|
||||
if (udotn /= 0.0_pReal) then
|
||||
Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
|
||||
dudotn_dt = udotn*damage_anisoBrittle_N(instance)/traction_n
|
||||
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
|
||||
dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
|
||||
dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* &
|
||||
lattice_Scleavage(m,n,3,index_myFamily+i,phase)
|
||||
endif
|
||||
|
||||
index = index + 3_pInt
|
||||
enddo
|
||||
enddo
|
||||
dLd_dTstar = math_Plain3333to99(dLd_dTstar3333)
|
||||
|
||||
end subroutine damage_anisoBrittle_LdAndItsTangent
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns local damage deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function damage_anisoBrittle_getFd(ipc, ip, el)
|
||||
use material, only: &
|
||||
mappingConstitutive, &
|
||||
phase_damageInstance, &
|
||||
damageState
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
integer(pInt) :: &
|
||||
phase, &
|
||||
constituent, &
|
||||
instance
|
||||
real(pReal), dimension(3,3) :: &
|
||||
damage_anisoBrittle_getFd
|
||||
|
||||
phase = mappingConstitutive(2,ipc,ip,el)
|
||||
constituent = mappingConstitutive(1,ipc,ip,el)
|
||||
instance = phase_damageInstance(phase)
|
||||
|
||||
damage_anisoBrittle_getFd = &
|
||||
reshape(damageState(phase)%state(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
|
||||
3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), &
|
||||
shape=[3,3])
|
||||
|
||||
end function damage_anisoBrittle_getFd
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns local damage deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function damage_anisoBrittle_getFd0(ipc, ip, el)
|
||||
use material, only: &
|
||||
mappingConstitutive, &
|
||||
phase_damageInstance, &
|
||||
damageState
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
integer(pInt) :: &
|
||||
phase, &
|
||||
constituent, &
|
||||
instance
|
||||
real(pReal), dimension(3,3) :: &
|
||||
damage_anisoBrittle_getFd0
|
||||
|
||||
phase = mappingConstitutive(2,ipc,ip,el)
|
||||
constituent = mappingConstitutive(1,ipc,ip,el)
|
||||
instance = phase_damageInstance(phase)
|
||||
|
||||
damage_anisoBrittle_getFd0 = &
|
||||
reshape(damageState(phase)%subState0(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
|
||||
3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), &
|
||||
shape=[3,3])
|
||||
|
||||
end function damage_anisoBrittle_getFd0
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns local damage deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function damage_anisoBrittle_getPartionedFd0(ipc, ip, el)
|
||||
use material, only: &
|
||||
mappingConstitutive, &
|
||||
phase_damageInstance, &
|
||||
damageState
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
integer(pInt) :: &
|
||||
phase, &
|
||||
constituent, &
|
||||
instance
|
||||
real(pReal), dimension(3,3) :: &
|
||||
damage_anisoBrittle_getPartionedFd0
|
||||
|
||||
phase = mappingConstitutive(2,ipc,ip,el)
|
||||
constituent = mappingConstitutive(1,ipc,ip,el)
|
||||
instance = phase_damageInstance(phase)
|
||||
|
||||
damage_anisoBrittle_getPartionedFd0 = &
|
||||
reshape(damageState(phase)%partionedState0(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
|
||||
3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), &
|
||||
shape=[3,3])
|
||||
|
||||
end function damage_anisoBrittle_getPartionedFd0
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns damage
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -473,29 +687,6 @@ function damage_anisoBrittle_getLocalDamage(ipc, ip, el)
|
|||
|
||||
end function damage_anisoBrittle_getLocalDamage
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns local damage deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function damage_anisoBrittle_getDamageStrain(ipc, ip, el)
|
||||
use material, only: &
|
||||
mappingConstitutive, &
|
||||
phase_damageInstance, &
|
||||
damageState
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), dimension(3,3) :: &
|
||||
damage_anisoBrittle_getDamageStrain
|
||||
|
||||
damage_anisoBrittle_getDamageStrain = &
|
||||
reshape(damageState(mappingConstitutive(2,ipc,ip,el))%state(3:11,mappingConstitutive(1,ipc,ip,el)), &
|
||||
shape=[3,3])
|
||||
|
||||
end function damage_anisoBrittle_getDamageStrain
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief return array of constitutive results
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -40,6 +40,10 @@ module thermal_adiabatic
|
|||
thermal_adiabatic_stateInit, &
|
||||
thermal_adiabatic_aTolState, &
|
||||
thermal_adiabatic_dotState, &
|
||||
thermal_adiabatic_LTAndItsTangent, &
|
||||
thermal_adiabatic_getFT, &
|
||||
thermal_adiabatic_getFT0, &
|
||||
thermal_adiabatic_getPartionedFT0, &
|
||||
thermal_adiabatic_getTemperature, &
|
||||
thermal_adiabatic_putTemperature, &
|
||||
thermal_adiabatic_postResults
|
||||
|
|
@ -283,6 +287,160 @@ subroutine thermal_adiabatic_dotState(Tstar_v, Lp, ipc, ip, el)
|
|||
|
||||
end subroutine thermal_adiabatic_dotState
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the velocity gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine thermal_adiabatic_LTAndItsTangent(LT, dLT_dTstar, Tstar_v, Lp, ipc, ip, el)
|
||||
use lattice, only: &
|
||||
lattice_massDensity, &
|
||||
lattice_specificHeat
|
||||
use material, only: &
|
||||
mappingConstitutive, &
|
||||
phase_thermalInstance, &
|
||||
thermalState
|
||||
use math, only: &
|
||||
math_Plain3333to99, &
|
||||
math_Mandel6to33, &
|
||||
math_I3
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), intent(in), dimension(6) :: &
|
||||
Tstar_v !< 2nd Piola-Kirchhoff stress
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
Lp !< plastic velocity gradient
|
||||
real(pReal), intent(out), dimension(3,3) :: &
|
||||
LT !< thermal velocity gradient
|
||||
real(pReal), intent(out), dimension(9,9) :: &
|
||||
dLT_dTstar !< derivative of LT with respect to Tstar (2nd-order tensor)
|
||||
real(pReal), dimension(3,3,3,3) :: &
|
||||
dLT_dTstar3333 !< derivative of LT with respect to Tstar (4th-order tensor)
|
||||
integer(pInt) :: &
|
||||
phase, &
|
||||
constituent, &
|
||||
i, j, k, l
|
||||
real(pReal) :: &
|
||||
Tdot
|
||||
|
||||
phase = mappingConstitutive(2,ipc,ip,el)
|
||||
constituent = mappingConstitutive(1,ipc,ip,el)
|
||||
|
||||
Tdot = 0.95_pReal &
|
||||
* sum(abs(math_Mandel6to33(Tstar_v))*Lp) &
|
||||
/ (lattice_massDensity(phase)*lattice_specificHeat(phase))
|
||||
|
||||
LT = Tdot*math_I3
|
||||
dLT_dTstar3333 = 0.0_pReal
|
||||
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
|
||||
dLT_dTstar3333(i,j,k,l) = dLT_dTstar3333(i,j,k,l) + Lp(k,l)*math_I3(i,j)
|
||||
|
||||
dLT_dTstar3333 = 0.95_pReal*dLT_dTstar3333/(lattice_massDensity(phase)*lattice_specificHeat(phase))
|
||||
dLT_dTstar = math_Plain3333to99(dLT_dTstar3333)
|
||||
|
||||
end subroutine thermal_adiabatic_LTAndItsTangent
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns local thermal deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function thermal_adiabatic_getFT(ipc, ip, el)
|
||||
use material, only: &
|
||||
mappingConstitutive, &
|
||||
thermalState
|
||||
use math, only: &
|
||||
math_I3
|
||||
use lattice, only: &
|
||||
lattice_referenceTemperature, &
|
||||
lattice_thermalExpansion33
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
integer(pInt) :: &
|
||||
phase, &
|
||||
constituent
|
||||
real(pReal), dimension(3,3) :: &
|
||||
thermal_adiabatic_getFT
|
||||
|
||||
phase = mappingConstitutive(2,ipc,ip,el)
|
||||
constituent = mappingConstitutive(1,ipc,ip,el)
|
||||
|
||||
thermal_adiabatic_getFT = math_I3 + &
|
||||
lattice_thermalExpansion33(1:3,1:3,phase)* &
|
||||
(thermalState(phase)%state(1,constituent) - lattice_referenceTemperature(phase))
|
||||
|
||||
end function thermal_adiabatic_getFT
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns local thermal deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function thermal_adiabatic_getFT0(ipc, ip, el)
|
||||
use material, only: &
|
||||
mappingConstitutive, &
|
||||
thermalState
|
||||
use math, only: &
|
||||
math_I3
|
||||
use lattice, only: &
|
||||
lattice_referenceTemperature, &
|
||||
lattice_thermalExpansion33
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
integer(pInt) :: &
|
||||
phase, &
|
||||
constituent
|
||||
real(pReal), dimension(3,3) :: &
|
||||
thermal_adiabatic_getFT0
|
||||
|
||||
phase = mappingConstitutive(2,ipc,ip,el)
|
||||
constituent = mappingConstitutive(1,ipc,ip,el)
|
||||
|
||||
thermal_adiabatic_getFT0 = math_I3 + &
|
||||
lattice_thermalExpansion33(1:3,1:3,phase)* &
|
||||
(thermalState(phase)%subState0(1,constituent) - lattice_referenceTemperature(phase))
|
||||
|
||||
end function thermal_adiabatic_getFT0
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns local thermal deformation gradient
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function thermal_adiabatic_getPartionedFT0(ipc, ip, el)
|
||||
use material, only: &
|
||||
mappingConstitutive, &
|
||||
thermalState
|
||||
use math, only: &
|
||||
math_I3
|
||||
use lattice, only: &
|
||||
lattice_referenceTemperature, &
|
||||
lattice_thermalExpansion33
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ipc, & !< grain number
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
integer(pInt) :: &
|
||||
phase, &
|
||||
constituent
|
||||
real(pReal), dimension(3,3) :: &
|
||||
thermal_adiabatic_getPartionedFT0
|
||||
|
||||
phase = mappingConstitutive(2,ipc,ip,el)
|
||||
constituent = mappingConstitutive(1,ipc,ip,el)
|
||||
|
||||
thermal_adiabatic_getPartionedFT0 = math_I3 + &
|
||||
lattice_thermalExpansion33(1:3,1:3,phase)* &
|
||||
(thermalState(phase)%partionedState0(1,constituent) - lattice_referenceTemperature(phase))
|
||||
|
||||
end function thermal_adiabatic_getPartionedFT0
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns temperature based on local damage model state layout
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
Loading…
Reference in New Issue