223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

changes related to intermediate configuration kinematics: 

- switched Fi and Li from state variables to crystallite variables

- Lp and Li are now work conjugate with the corresponding mandel stresses defined in their respective configuration

- T, Lp and Li need to return tangents wrt Fi arising from the convection of the material frame due to Fi

- Updated analytic jacobian to take into account tangents wrt Fi

- Updated Lp and Li residual jacobians to take into account tangents wrt Fi
This commit is contained in:
Pratheek Shanthraj 2015-03-06 13:09:00 +00:00
parent ecc8e2c5b4
commit bbb5ff6ae9
8 changed files with 449 additions and 981 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
code/CPFEM.f90
View File

@ -151,6 +151,8 @@ subroutine CPFEM_init
crystallite_F0, & crystallite_F0, &
crystallite_Fp0, & crystallite_Fp0, &
crystallite_Lp0, & crystallite_Lp0, &
crystallite_Fi0, &
crystallite_Li0, &
crystallite_dPdF0, & crystallite_dPdF0, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_localPlasticity crystallite_localPlasticity
@ -191,10 +193,18 @@ subroutine CPFEM_init
read (777,rec=1) crystallite_Fp0 read (777,rec=1) crystallite_Fp0
close (777) close (777)
call IO_read_realFile(777,'convergedFi',modelName,size(crystallite_Fi0))
read (777,rec=1) crystallite_Fi0
close (777)
call IO_read_realFile(777,'convergedLp',modelName,size(crystallite_Lp0)) call IO_read_realFile(777,'convergedLp',modelName,size(crystallite_Lp0))
read (777,rec=1) crystallite_Lp0 read (777,rec=1) crystallite_Lp0
close (777) close (777)
call IO_read_realFile(777,'convergedLi',modelName,size(crystallite_Li0))
read (777,rec=1) crystallite_Li0
close (777)
call IO_read_realFile(777,'convergeddPdF',modelName,size(crystallite_dPdF0)) call IO_read_realFile(777,'convergeddPdF',modelName,size(crystallite_dPdF0))
read (777,rec=1) crystallite_dPdF0 read (777,rec=1) crystallite_dPdF0
close (777) close (777)
@ -314,8 +324,12 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
crystallite_F0, & crystallite_F0, &
crystallite_Fp0, & crystallite_Fp0, &
crystallite_Fp, & crystallite_Fp, &
crystallite_Fi0, &
crystallite_Fi, &
crystallite_Lp0, & crystallite_Lp0, &
crystallite_Lp, & crystallite_Lp, &
crystallite_Li0, &
crystallite_Li, &
crystallite_dPdF0, & crystallite_dPdF0, &
crystallite_dPdF, & crystallite_dPdF, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
@ -400,6 +414,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...)
crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation
crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation
crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity
crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness
crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress
@ -439,10 +455,18 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
write (777,rec=1) crystallite_Fp0 write (777,rec=1) crystallite_Fp0
close (777) close (777)
call IO_write_jobRealFile(777,'convergedFi',size(crystallite_Fi0))
write (777,rec=1) crystallite_Fi0
close (777)
call IO_write_jobRealFile(777,'convergedLp',size(crystallite_Lp0)) call IO_write_jobRealFile(777,'convergedLp',size(crystallite_Lp0))
write (777,rec=1) crystallite_Lp0 write (777,rec=1) crystallite_Lp0
close (777) close (777)
call IO_write_jobRealFile(777,'convergedLi',size(crystallite_Li0))
write (777,rec=1) crystallite_Li0
close (777)
call IO_write_jobRealFile(777,'convergeddPdF',size(crystallite_dPdF0)) call IO_write_jobRealFile(777,'convergeddPdF',size(crystallite_dPdF0))
write (777,rec=1) crystallite_dPdF0 write (777,rec=1) crystallite_dPdF0
close (777) close (777)

348
code/constitutive.f90
View File

@ -28,10 +28,6 @@ module constitutive
constitutive_microstructure, & constitutive_microstructure, &
constitutive_LpAndItsTangent, & constitutive_LpAndItsTangent, &
constitutive_LiAndItsTangent, & constitutive_LiAndItsTangent, &
constitutive_getFi, &
constitutive_putFi, &
constitutive_getFi0, &
constitutive_getPartionedFi0, &
constitutive_TandItsTangent, & constitutive_TandItsTangent, &
constitutive_collectDotState, & constitutive_collectDotState, &
constitutive_collectDeltaState, & constitutive_collectDeltaState, &
@ -678,11 +674,15 @@ end subroutine constitutive_microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el) subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
math_identity2nd math_transpose33, &
math_mul33x33, &
math_Mandel6to33, &
math_Mandel33to6, &
math_Plain99to3333
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
phase_plasticityInstance, & phase_plasticityInstance, &
@ -719,42 +719,66 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el)
el !< element number el !< element number
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress Tstar_v !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
real(pReal), intent(out), dimension(9,9) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dLp_dTstar !< derivative of Lp with respect to Tstar (4th-order tensor) 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
real(pReal), dimension(9,9) :: &
dLp_dMstar !< derivative of Lp with respect to Mstar (4th-order tensor)
real(pReal), dimension(3,3) :: &
temp_33
integer(pInt) :: &
i, j
Mstar_v = math_Mandel33to6(math_mul33x33(math_mul33x33(math_transpose33(Fi),Fi), &
math_Mandel6to33(Tstar_v)))
select case (phase_plasticity(material_phase(ipc,ip,el))) select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_NONE_ID) case (PLASTICITY_NONE_ID)
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dTstar = 0.0_pReal dLp_dMstar = 0.0_pReal
case (PLASTICITY_J2_ID) case (PLASTICITY_J2_ID)
call plastic_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,ipc,ip,el) call plastic_j2_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) case (PLASTICITY_PHENOPOWERLAW_ID)
call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,ipc,ip,el) call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) case (PLASTICITY_NONLOCAL_ID)
call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v, &
constitutive_getTemperature(ipc,ip,el), & constitutive_getTemperature(ipc,ip,el), &
ipc,ip,el) ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID) case (PLASTICITY_DISLOTWIN_ID)
call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v, &
constitutive_getTemperature(ipc,ip,el), & constitutive_getTemperature(ipc,ip,el), &
ipc,ip,el) ipc,ip,el)
case (PLASTICITY_DISLOKMC_ID) case (PLASTICITY_DISLOKMC_ID)
call plastic_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & call plastic_dislokmc_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v, &
constitutive_getTemperature(ipc,ip,el), & constitutive_getTemperature(ipc,ip,el), &
ipc,ip,el) ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) case (PLASTICITY_DISLOUCLA_ID)
call plastic_disloucla_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & call plastic_disloucla_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v, &
constitutive_getTemperature(ipc,ip,el), & constitutive_getTemperature(ipc,ip,el), &
ipc,ip,el) ipc,ip,el)
case (PLASTICITY_TITANMOD_ID) case (PLASTICITY_TITANMOD_ID)
call plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & call plastic_titanmod_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v, &
constitutive_getTemperature(ipc,ip,el), & constitutive_getTemperature(ipc,ip,el), &
ipc,ip,el) ipc,ip,el)
end select end select
dLp_dTstar3333 = math_Plain99to3333(dLp_dMstar)
temp_33 = math_mul33x33(Fi,math_Mandel6to33(Tstar_v))
do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt
dLp_dFi3333(i,j,1:3,1:3) = math_mul33x33(temp_33,math_transpose33(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))
enddo; enddo
temp_33 = math_mul33x33(math_transpose33(Fi),Fi)
do i = 1_pInt, 3_pInt; do 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))
enddo; enddo
end subroutine constitutive_LpAndItsTangent end subroutine constitutive_LpAndItsTangent
@ -762,9 +786,15 @@ end subroutine constitutive_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, Tstar_v, Lp, ipc, ip, el) subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v, Fi, Lp, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: &
math_I3, &
math_inv33, &
math_det33, &
math_transpose33, &
math_mul33x33
use material, only: & use material, only: &
phase_damage, & phase_damage, &
phase_thermal, & phase_thermal, &
@ -787,254 +817,68 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, Tstar_v, Lp, ipc, ip, el
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress Tstar_v !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fi, & !< intermediate deformation gradient
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
Li !< intermediate velocity gradient Li !< intermediate velocity gradient
real(pReal), intent(out), dimension(9,9) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar !< derivative of Li with respect to Tstar (2nd-order tensor) dLi_dTstar3333, & !< derivative of Li with respect to Tstar (4th-order tensor)
dLi_dFi3333
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Li_temp !< intermediate velocity gradient Li_temp !< intermediate velocity gradient
real(pReal), dimension(9,9) :: & real(pReal), dimension(3,3,3,3) :: &
dLi_dTstar_temp !< derivative of Li with respect to Tstar (4th-order tensor) dLi_dTstar_temp
real(pReal), dimension(3,3) :: &
FiInv, &
temp_33
real(pReal) :: &
detFi
integer(pInt) :: &
i, j
Li = 0.0_pReal Li = 0.0_pReal
dLi_dTstar = 0.0_pReal dLi_dTstar3333 = 0.0_pReal
dLi_dFi3333 = 0.0_pReal
select case (phase_damage(material_phase(ipc,ip,el))) select case (phase_damage(material_phase(ipc,ip,el)))
case (LOCAL_DAMAGE_anisoBrittle_ID) case (LOCAL_DAMAGE_anisoBrittle_ID)
call damage_anisoBrittle_LdAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, ipc, ip, el) call damage_anisoBrittle_LdAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, ipc, ip, el)
Li = Li + Li_temp Li = Li + Li_temp
dLi_dTstar = dLi_dTstar + dLi_dTstar_temp dLi_dTstar3333 = dLi_dTstar3333 + dLi_dTstar_temp
case (LOCAL_DAMAGE_anisoDuctile_ID) case (LOCAL_DAMAGE_anisoDuctile_ID)
call damage_anisoDuctile_LdAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, ipc, ip, el) call damage_anisoDuctile_LdAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, ipc, ip, el)
Li = Li + Li_temp Li = Li + Li_temp
dLi_dTstar = dLi_dTstar + dLi_dTstar_temp dLi_dTstar3333 = dLi_dTstar3333 + dLi_dTstar_temp
end select end select
select case (phase_thermal(material_phase(ipc,ip,el))) select case (phase_thermal(material_phase(ipc,ip,el)))
case (LOCAL_THERMAL_adiabatic_ID) case (LOCAL_THERMAL_adiabatic_ID)
call thermal_adiabatic_LTAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, Lp, ipc, ip, el) call thermal_adiabatic_LTAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, Lp, ipc, ip, el)
Li = Li + Li_temp Li = Li + Li_temp
dLi_dTstar = dLi_dTstar + dLi_dTstar_temp dLi_dTstar3333 = dLi_dTstar3333 + dLi_dTstar_temp
end select end select
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)
do i = 1_pInt, 3_pInt; do 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)
enddo; enddo
end subroutine constitutive_LiAndItsTangent 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_DAMAGE_anisoDuctile_ID, &
LOCAL_THERMAL_adiabatic_ID
use damage_anisoBrittle, only: &
damage_anisoBrittle_getFd
use damage_anisoDuctile, only: &
damage_anisoDuctile_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))
case (LOCAL_DAMAGE_anisoDuctile_ID)
constitutive_getFi = math_mul33x33(constitutive_getFi,damage_anisoDuctile_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
!--------------------------------------------------------------------------------------------------
subroutine constitutive_putFi(Tstar_v, Lp, dt, ipc, ip, el)
use prec, only: &
pReal
use material, only: &
phase_damage, &
phase_thermal, &
material_phase, &
LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_DAMAGE_anisoDuctile_ID, &
LOCAL_THERMAL_adiabatic_ID
use damage_anisoBrittle, only: &
damage_anisoBrittle_putFd
use damage_anisoDuctile, only: &
damage_anisoDuctile_putFd
use thermal_adiabatic, only: &
thermal_adiabatic_putFT
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(in) :: &
dt
select case (phase_damage(material_phase(ipc,ip,el)))
case (LOCAL_DAMAGE_anisoBrittle_ID)
call damage_anisoBrittle_putFd (Tstar_v, dt, ipc, ip, el)
case (LOCAL_DAMAGE_anisoDuctile_ID)
call damage_anisoDuctile_putFd (Tstar_v, dt, ipc, ip, el)
end select
select case (phase_thermal(material_phase(ipc,ip,el)))
case (LOCAL_THERMAL_adiabatic_ID)
call thermal_adiabatic_putFT (Tstar_v, Lp, dt, ipc, ip, el)
end select
end subroutine constitutive_putFi
!--------------------------------------------------------------------------------------------------
!> @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_DAMAGE_anisoDuctile_ID, &
LOCAL_THERMAL_adiabatic_ID
use damage_anisoBrittle, only: &
damage_anisoBrittle_getFd0
use damage_anisoDuctile, only: &
damage_anisoDuctile_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))
case (LOCAL_DAMAGE_anisoDuctile_ID)
constitutive_getFi0 = math_mul33x33(constitutive_getFi0,damage_anisoDuctile_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_DAMAGE_anisoDuctile_ID, &
LOCAL_THERMAL_adiabatic_ID
use damage_anisoBrittle, only: &
damage_anisoBrittle_getPartionedFd0
use damage_anisoDuctile, only: &
damage_anisoDuctile_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))
case (LOCAL_DAMAGE_anisoDuctile_ID)
constitutive_getPartionedFi0 = math_mul33x33(constitutive_getPartionedFi0, &
damage_anisoDuctile_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 !> @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 !> the elastic deformation gradient depending on the selected elastic law (so far no case switch
!! because only hooke is implemented !! because only hooke is implemented
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el) subroutine constitutive_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
@ -1044,13 +888,15 @@ subroutine constitutive_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el)
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
T !< 2nd Piola-Kirchhoff stress tensor T !< 2nd Piola-Kirchhoff stress tensor
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dT_dFe !< derivative of 2nd P-K stress with respect to elastic deformation gradient 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
call constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el) call constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
end subroutine constitutive_TandItsTangent end subroutine constitutive_TandItsTangent
@ -1060,7 +906,7 @@ end subroutine constitutive_TandItsTangent
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic deformation gradient using hookes law !> the elastic deformation gradient using hookes law
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el) subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only : & use math, only : &
@ -1083,21 +929,28 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el)
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
T !< 2nd Piola-Kirchhoff stress tensor T !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dT_dFe !< dT/dFe 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
integer(pInt) :: i, j, k, l integer(pInt) :: i, j
real(pReal), dimension(3,3) :: E
real(pReal), dimension(3,3,3,3) :: C real(pReal), dimension(3,3,3,3) :: C
C = math_Mandel66to3333(constitutive_damagedC(ipc,ip,el)) C = math_Mandel66to3333(constitutive_damagedC(ipc,ip,el)) !< elastic stiffness in lattice configuration
T = math_mul3333xx33(C,0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)) E = 0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
T = math_mul3333xx33(C,math_mul33x33(math_mul33x33(math_transpose33(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
dT_dFe = 0.0_pReal dT_dFe = 0.0_pReal
forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt, k=1_pInt:3_pInt, l=1_pInt:3_pInt) & forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt)
dT_dFe(i,j,k,l) = sum(C(i,j,l,1:3)*Fe(k,1:3)) ! dT*_ij/dFe_kl dT_dFe(i,j,1:3,1:3) = &
math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),math_transpose33(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
end forall
end subroutine constitutive_hooke_TandItsTangent end subroutine constitutive_hooke_TandItsTangent
@ -1470,9 +1323,12 @@ function constitutive_getDamageDiffusion33(ipc, ip, el)
material_phase, & material_phase, &
phase_damage, & phase_damage, &
LOCAL_DAMAGE_isoBrittle_ID, & LOCAL_DAMAGE_isoBrittle_ID, &
LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_DAMAGE_phaseField_ID LOCAL_DAMAGE_phaseField_ID
use damage_isoBrittle, only: & use damage_isoBrittle, only: &
damage_isoBrittle_getDamageDiffusion33 damage_isoBrittle_getDamageDiffusion33
use damage_anisoBrittle, only: &
damage_anisoBrittle_getDamageDiffusion33
use damage_phaseField, only: & use damage_phaseField, only: &
damage_phaseField_getDamageDiffusion33 damage_phaseField_getDamageDiffusion33
@ -1487,6 +1343,8 @@ function constitutive_getDamageDiffusion33(ipc, ip, el)
select case(phase_damage(material_phase(ipc,ip,el))) select case(phase_damage(material_phase(ipc,ip,el)))
case (LOCAL_DAMAGE_isoBrittle_ID) case (LOCAL_DAMAGE_isoBrittle_ID)
constitutive_getDamageDiffusion33 = damage_isoBrittle_getDamageDiffusion33(ipc, ip, el) constitutive_getDamageDiffusion33 = damage_isoBrittle_getDamageDiffusion33(ipc, ip, el)
case (LOCAL_DAMAGE_anisoBrittle_ID)
constitutive_getDamageDiffusion33 = damage_anisoBrittle_getDamageDiffusion33(ipc, ip, el)
case (LOCAL_DAMAGE_phaseField_ID) case (LOCAL_DAMAGE_phaseField_ID)
constitutive_getDamageDiffusion33 = damage_phaseField_getDamageDiffusion33(ipc, ip, el) constitutive_getDamageDiffusion33 = damage_phaseField_getDamageDiffusion33(ipc, ip, el)
case default case default

355
code/crystallite.f90
View File

@ -43,21 +43,30 @@ module crystallite
crystallite_Fp, & !< current plastic def grad (end of converged time step) crystallite_Fp, & !< current plastic def grad (end of converged time step)
crystallite_Fp0, & !< plastic def grad at start of FE inc crystallite_Fp0, & !< plastic def grad at start of FE inc
crystallite_partionedFp0,& !< plastic def grad at start of homog inc crystallite_partionedFp0,& !< plastic def grad at start of homog inc
crystallite_Fi, & !< current intermediate def grad (end of converged time step)
crystallite_Fi0, & !< intermediate def grad at start of FE inc
crystallite_partionedFi0,& !< intermediate def grad at start of homog inc
crystallite_F0, & !< def grad at start of FE inc crystallite_F0, & !< def grad at start of FE inc
crystallite_partionedF, & !< def grad to be reached at end of homog inc crystallite_partionedF, & !< def grad to be reached at end of homog inc
crystallite_partionedF0, & !< def grad at start of homog inc crystallite_partionedF0, & !< def grad at start of homog inc
crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step) crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step)
crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc
crystallite_partionedLp0,& !< plastic velocity grad at start of homog inc crystallite_partionedLp0,& !< plastic velocity grad at start of homog inc
crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step)
crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc
crystallite_partionedLi0,& !< intermediate velocity grad at start of homog inc
crystallite_P !< 1st Piola-Kirchhoff stress per grain crystallite_P !< 1st Piola-Kirchhoff stress per grain
real(pReal), dimension(:,:,:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:,:,:), allocatable, private :: &
crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc
crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step)
crystallite_subFp0,& !< plastic def grad at start of crystallite inc crystallite_subFp0,& !< plastic def grad at start of crystallite inc
crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step)
crystallite_subFi0,& !< intermediate def grad at start of crystallite inc
crystallite_subF, & !< def grad to be reached at end of crystallite inc crystallite_subF, & !< def grad to be reached at end of crystallite inc
crystallite_subF0, & !< def grad at start of crystallite inc crystallite_subF0, & !< def grad at start of crystallite inc
crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc
crystallite_subLi0,& !< intermediate velocity grad at start of crystallite inc
crystallite_disorientation !< disorientation between two neighboring ips (only calculated for single grain IPs) crystallite_disorientation !< disorientation between two neighboring ips (only calculated for single grain IPs)
real(pReal), dimension(:,:,:,:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:,:,:,:), allocatable, public :: &
crystallite_dPdF, & !< current individual dPdF per grain (end of converged time step) crystallite_dPdF, & !< current individual dPdF per grain (end of converged time step)
@ -93,7 +102,9 @@ module crystallite
defgrad_ID, & defgrad_ID, &
fe_ID, & fe_ID, &
fp_ID, & fp_ID, &
fi_ID, &
lp_ID, & lp_ID, &
li_ID, &
e_ID, & e_ID, &
ee_ID, & ee_ID, &
p_ID, & p_ID, &
@ -225,12 +236,21 @@ subroutine crystallite_init
allocate(crystallite_subFp0(3,3,gMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subFp0(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Fp(3,3,gMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_Fp(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_invFp(3,3,gMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_invFp(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Fi0(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_partionedFi0(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subFi0(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Fi(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_invFi(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Fe(3,3,gMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_Fe(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subFe0(3,3,gMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subFe0(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Lp0(3,3,gMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_Lp0(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_partionedLp0(3,3,gMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_partionedLp0(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subLp0(3,3,gMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subLp0(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Lp(3,3,gMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_Lp(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Li0(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_partionedLi0(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subLi0(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Li(3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_dPdF0(3,3,3,3,gMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_dPdF0(3,3,3,3,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_partioneddPdF0(3,3,3,3,gMax,iMax,eMax),source=0.0_pReal) allocate(crystallite_partioneddPdF0(3,3,3,3,gMax,iMax,eMax),source=0.0_pReal)
@ -311,8 +331,12 @@ subroutine crystallite_init
crystallite_outputID(output,section) = fe_ID crystallite_outputID(output,section) = fe_ID
case ('fp') case ('fp')
crystallite_outputID(output,section) = fp_ID crystallite_outputID(output,section) = fp_ID
case ('fi')
crystallite_outputID(output,section) = fi_ID
case ('lp') case ('lp')
crystallite_outputID(output,section) = lp_ID crystallite_outputID(output,section) = lp_ID
case ('li')
crystallite_outputID(output,section) = li_ID
case ('e') case ('e')
crystallite_outputID(output,section) = e_ID crystallite_outputID(output,section) = e_ID
case ('ee') case ('ee')
@ -345,7 +369,7 @@ subroutine crystallite_init
mySize = 4_pInt mySize = 4_pInt
case(eulerangles_ID) case(eulerangles_ID)
mySize = 3_pInt mySize = 3_pInt
case(defgrad_ID,fe_ID,fp_ID,lp_ID,e_ID,ee_ID,p_ID,s_ID) case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,e_ID,ee_ID,p_ID,s_ID)
mySize = 9_pInt mySize = 9_pInt
case(elasmatrix_ID) case(elasmatrix_ID)
mySize = 36_pInt mySize = 36_pInt
@ -389,10 +413,12 @@ subroutine crystallite_init
myNgrains = homogenization_Ngrains(mesh_element(3,e)) ! look up homogenization-->grainCount myNgrains = homogenization_Ngrains(mesh_element(3,e)) ! look up homogenization-->grainCount
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1_pInt:myNgrains) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1_pInt:myNgrains)
crystallite_Fp0(1:3,1:3,g,i,e) = math_EulerToR(material_EulerAngles(1:3,g,i,e)) ! plastic def gradient reflects init orientation crystallite_Fp0(1:3,1:3,g,i,e) = math_EulerToR(material_EulerAngles(1:3,g,i,e)) ! plastic def gradient reflects init orientation
crystallite_Fi0(1:3,1:3,g,i,e) = math_I3
crystallite_F0(1:3,1:3,g,i,e) = math_I3 crystallite_F0(1:3,1:3,g,i,e) = math_I3
crystallite_localPlasticity(g,i,e) = phase_localPlasticity(material_phase(g,i,e)) crystallite_localPlasticity(g,i,e) = phase_localPlasticity(material_phase(g,i,e))
crystallite_Fe(1:3,1:3,g,i,e) = math_transpose33(crystallite_Fp0(1:3,1:3,g,i,e)) crystallite_Fe(1:3,1:3,g,i,e) = math_transpose33(crystallite_Fp0(1:3,1:3,g,i,e))
crystallite_Fp(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) crystallite_Fp(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e)
crystallite_Fi(1:3,1:3,g,i,e) = crystallite_Fi0(1:3,1:3,g,i,e)
crystallite_requested(g,i,e) = .true. crystallite_requested(g,i,e) = .true.
endforall endforall
enddo enddo
@ -401,6 +427,7 @@ subroutine crystallite_init
if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601_pInt) ! exit if nonlocal but no ping-pong if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601_pInt) ! exit if nonlocal but no ping-pong
crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFp0 = crystallite_Fp0
crystallite_partionedFi0 = crystallite_Fi0
crystallite_partionedF0 = crystallite_F0 crystallite_partionedF0 = crystallite_F0
crystallite_partionedF = crystallite_F0 crystallite_partionedF = crystallite_F0
@ -431,19 +458,27 @@ subroutine crystallite_init
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp: ', shape(crystallite_Fp) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp: ', shape(crystallite_Fp)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fi: ', shape(crystallite_Fi)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp: ', shape(crystallite_Lp) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp: ', shape(crystallite_Lp)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Li: ', shape(crystallite_Li)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_F0: ', shape(crystallite_F0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_F0: ', shape(crystallite_F0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp0: ', shape(crystallite_Fp0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp0: ', shape(crystallite_Fp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fi0: ', shape(crystallite_Fi0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp0: ', shape(crystallite_Lp0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp0: ', shape(crystallite_Lp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Li0: ', shape(crystallite_Li0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF: ', shape(crystallite_partionedF) write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF: ', shape(crystallite_partionedF)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFi0: ', shape(crystallite_partionedFi0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLi0: ', shape(crystallite_partionedLi0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF: ', shape(crystallite_subF) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF: ', shape(crystallite_subF)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF0: ', shape(crystallite_subF0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF0: ', shape(crystallite_subF0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFe0: ', shape(crystallite_subFe0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFe0: ', shape(crystallite_subFe0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFp0: ', shape(crystallite_subFp0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFp0: ', shape(crystallite_subFp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFi0: ', shape(crystallite_subFi0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subLp0: ', shape(crystallite_subLp0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subLp0: ', shape(crystallite_subLp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subLi0: ', shape(crystallite_subLi0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_P: ', shape(crystallite_P) write(6,'(a35,1x,7(i8,1x))') 'crystallite_P: ', shape(crystallite_P)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar0_v: ', shape(crystallite_Tstar0_v) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar0_v: ', shape(crystallite_Tstar0_v)
@ -514,6 +549,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
math_mul66x6, & math_mul66x6, &
math_Mandel6to33, & math_Mandel6to33, &
math_Mandel33to6, & math_Mandel33to6, &
math_Plain3333to99, &
math_Plain99to3333, &
math_I3, & math_I3, &
math_mul3333xx3333, & math_mul3333xx3333, &
math_mul33xx33, & math_mul33xx33, &
@ -541,10 +578,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
use constitutive, only: & use constitutive, only: &
constitutive_TandItsTangent, & constitutive_TandItsTangent, &
constitutive_LpAndItsTangent, & constitutive_LpAndItsTangent, &
constitutive_LiAndItsTangent, & constitutive_LiAndItsTangent
constitutive_getFi, &
constitutive_getFi0, &
constitutive_getPartionedFi0
implicit none implicit none
logical, intent(in) :: & logical, intent(in) :: &
@ -564,8 +598,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
F_backup, & F_backup, &
Fp_backup, & Fp_backup, &
InvFp_backup, & InvFp_backup, &
Fi_backup, &
InvFi_backup, &
Fe_backup, & Fe_backup, &
Lp_backup, & Lp_backup, &
Li_backup, &
P_backup P_backup
real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
Tstar_v_backup Tstar_v_backup
@ -586,19 +623,16 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
convergenceFlag_backup convergenceFlag_backup
! local variables used for calculating analytic Jacobian ! local variables used for calculating analytic Jacobian
real(pReal) :: detFi real(pReal), dimension(3,3) :: temp_33
real(pReal), dimension(3,3) :: temp_33, &
Fi, &
invFi, &
invFi0
real(pReal), dimension(3,3,3,3) :: dSdFe, & real(pReal), dimension(3,3,3,3) :: dSdFe, &
dSdF, & dSdF, &
dSdFiInv, & dSdFi, &
junk2, &
dLidS, & dLidS, &
dLidFi, &
dLpdS, & dLpdS, &
dLpdFi, &
dFidS, &
dFpinvdF, & dFpinvdF, &
dFiinvdF, &
rhs_3333, & rhs_3333, &
lhs_3333, & lhs_3333, &
temp_3333 temp_3333
@ -615,8 +649,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
math_transpose33(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e)) math_transpose33(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp0', & write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp0', &
math_transpose33(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e)) math_transpose33(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fi0', &
math_transpose33(crystallite_partionedFi0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp0', & write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp0', &
math_transpose33(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e)) math_transpose33(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Li0', &
math_transpose33(crystallite_partionedLi0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> F ', & write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> F ', &
math_transpose33(crystallite_partionedF(1:3,1:3,debug_g,debug_i,debug_e)) math_transpose33(crystallite_partionedF(1:3,1:3,debug_g,debug_i,debug_e))
endif endif
@ -640,12 +678,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
vacancyState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) vacancyState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e))
crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_partionedFp0(1:3,1:3,g,i,e) ! ...plastic def grad crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_partionedFp0(1:3,1:3,g,i,e) ! ...plastic def grad
crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_partionedLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_partionedLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad
crystallite_subFi0(1:3,1:3,g,i,e) = crystallite_partionedFi0(1:3,1:3,g,i,e) ! ...intermediate def grad
crystallite_subLi0(1:3,1:3,g,i,e) = crystallite_partionedLi0(1:3,1:3,g,i,e) ! ...intermediate velocity grad
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_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_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_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(math_mul33x33(crystallite_subF0(1:3,1:3,g,i,e), & 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(crystallite_subFp0(1:3,1:3,g,i,e))), &
math_inv33(constitutive_getFi0(g,i,e))) ! only needed later on for stiffness calculation math_inv33(crystallite_subFi0(1:3,1:3,g,i,e)))! only needed later on for stiffness calculation
crystallite_subFrac(g,i,e) = 0.0_pReal crystallite_subFrac(g,i,e) = 0.0_pReal
crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst
crystallite_todo(g,i,e) = .true. crystallite_todo(g,i,e) = .true.
@ -912,10 +952,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_subF0(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ...def grad crystallite_subF0(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ...def grad
!$OMP FLUSH(crystallite_subF0) !$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_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient
crystallite_subLi0(1:3,1:3,g,i,e) = crystallite_Li(1:3,1:3,g,i,e) ! ...intermediate velocity gradient
crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad
crystallite_subFi0(1:3,1:3,g,i,e) = crystallite_Fi(1:3,1:3,g,i,e) ! ...intermediate def grad
crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(crystallite_subF (1:3,1:3,g,i,e), & 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)), & crystallite_invFp(1:3,1:3,g,i,e)), &
math_inv33(constitutive_getFi(g,i,e))) ! only needed later on for stiffness calculation crystallite_invFi(1:3,1:3,g,i,e)) ! only needed later on for stiffness calculation
!if abbrevation, make c and p private in omp !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))%subState0(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) plasticState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e))
@ -967,7 +1009,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMP FLUSH(crystallite_Fp) !$OMP FLUSH(crystallite_Fp)
crystallite_invFp(1:3,1:3,g,i,e) = math_inv33(crystallite_Fp(1:3,1:3,g,i,e)) crystallite_invFp(1:3,1:3,g,i,e) = math_inv33(crystallite_Fp(1:3,1:3,g,i,e))
!$OMP FLUSH(crystallite_invFp) !$OMP FLUSH(crystallite_invFp)
crystallite_Fi(1:3,1:3,g,i,e) = crystallite_subFi0(1:3,1:3,g,i,e) ! ...intermediate def grad
!$OMP FLUSH(crystallite_Fi)
crystallite_invFi(1:3,1:3,g,i,e) = math_inv33(crystallite_Fi(1:3,1:3,g,i,e))
!$OMP FLUSH(crystallite_invFi)
crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad
crystallite_Li(1:3,1:3,g,i,e) = crystallite_subLi0(1:3,1:3,g,i,e) ! ...intermediate velocity grad
plasticState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = & plasticState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e)) plasticState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e))
damageState( mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = & damageState( mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
@ -1003,10 +1050,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
* (crystallite_partionedF(1:3,1:3,g,i,e) & * (crystallite_partionedF(1:3,1:3,g,i,e) &
- crystallite_partionedF0(1:3,1:3,g,i,e)) - crystallite_partionedF0(1:3,1:3,g,i,e))
!$OMP FLUSH(crystallite_subF) !$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_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)), & crystallite_invFp(1:3,1:3,g,i,e)), &
math_inv33(constitutive_getFi(g,i,e))) crystallite_invFi(1:3,1:3,g,i,e))
crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(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 crystallite_converged(g,i,e) = .false. ! start out non-converged
endif endif
@ -1083,11 +1129,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
e,'(',mesh_element(1,e),')',i,g e,'(',mesh_element(1,e),')',i,g
invFp = math_inv33(crystallite_partionedFp0(1:3,1:3,g,i,e)) invFp = math_inv33(crystallite_partionedFp0(1:3,1:3,g,i,e))
Fe_guess = math_mul33x33(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))) math_inv33(crystallite_partionedFi0(1:3,1:3,g,i,e)))
call constitutive_TandItsTangent(Tstar, junk2, Fe_guess,g,i,e) call constitutive_TandItsTangent(Tstar,dSdFe,dSdFi,Fe_guess,crystallite_partionedFi0(1:3,1:3,g,i,e),g,i,e)
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(Fe_guess,constitutive_getPartionedFi0(g,i,e)), & crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,g,i,e), invFp), &
math_mul33x33(Tstar,transpose(invFp)))* & math_mul33x33(Tstar,transpose(invFp)))
math_det33(constitutive_getPartionedFi0(g,i,e))
endif endif
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
@ -1097,8 +1142,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
math_transpose33(crystallite_P(1:3,1:3,g,i,e))*1.0e-6_pReal math_transpose33(crystallite_P(1:3,1:3,g,i,e))*1.0e-6_pReal
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp', & write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp', &
math_transpose33(crystallite_Fp(1:3,1:3,g,i,e)) math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fi', &
math_transpose33(crystallite_Fi(1:3,1:3,g,i,e))
write(6,'(a,/,3(12x,3(f14.9,1x)/),/)') '<< CRYST >> Lp', & write(6,'(a,/,3(12x,3(f14.9,1x)/),/)') '<< CRYST >> Lp', &
math_transpose33(crystallite_Lp(1:3,1:3,g,i,e)) math_transpose33(crystallite_Lp(1:3,1:3,g,i,e))
write(6,'(a,/,3(12x,3(f14.9,1x)/),/)') '<< CRYST >> Li', &
math_transpose33(crystallite_Li(1:3,1:3,g,i,e))
flush(6) flush(6)
endif endif
enddo enddo
@ -1113,37 +1162,44 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! --- ANALYTIC JACOBIAN --- ! --- ANALYTIC JACOBIAN ---
!$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFiInv,dLpdS,dFpinvdF,dFiinvdF,dLidS,rhs_3333,lhs_3333,& !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,&
!$OMP Fi,invFi,invFi0,detFi,temp_99,temp_33,temp_3333,myNgrains,error) !$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,myNgrains,error)
elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) 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 do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1_pInt,myNgrains do g = 1_pInt,myNgrains
Fi = constitutive_getFi(g,i,e) call constitutive_TandItsTangent(temp_33,dSdFe,dSdFi,crystallite_Fe(1:3,1:3,g,i,e), &
detFi = math_det33(Fi) crystallite_Fi(1:3,1:3,g,i,e),g,i,e) ! call constitutive law to calculate elastic stress tangent
invFi = math_inv33(Fi)
invFi0 = math_inv33(constitutive_getFi0(g,i,e))
call constitutive_TandItsTangent(temp_33,dSdFe,crystallite_Fe(1:3,1:3,g,i,e),g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
dSdFe(1:3,1:3,o,p) = math_mul33x33(invFi,math_mul33x33(dSdFe(1:3,1:3,o,p), &
math_transpose33(invFi)))*detFi
dSdFiInv = 0.0_pReal
temp_33 = math_mul33x33(temp_33,math_transpose33(invFi))*detFi
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
dSdFiInv(o,1:3,p,1:3) = dSdFiInv(o,1:3,p,1:3) + math_I3(o,p)*temp_33
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
dSdFiInv(1:3,o,p,1:3) = dSdFiInv(1:3,o,p,1:3) + math_I3(o,p)*math_transpose33(temp_33)
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
dSdFiInv(1:3,1:3,o,p) = dSdFiInv(1:3,1:3,o,p) - math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e))*invFi(p,o)
call constitutive_LpAndItsTangent(temp_33,temp_99,crystallite_Tstar_v(1:6,g,i,e), & call constitutive_LiAndItsTangent(temp_33,dLidS,dLidFi,crystallite_Tstar_v(1:6,g,i,e), &
g,i,e) crystallite_Fi(1:3,1:3,g,i,e),crystallite_Lp(1:3,1:3,g,i,e), &
dLpdS = reshape(temp_99,shape=[3,3,3,3]) g,i,e) ! call constitutive law to calculate Li tangent in lattice configuration
call constitutive_LiAndItsTangent(temp_33,temp_99,crystallite_Tstar_v(1:6,g,i,e), & temp_33 = math_inv33(crystallite_subFi0(1:3,1:3,g,i,e))
crystallite_Lp(1:3,1:3,g,i,e), g,i,e) lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
dLidS = reshape(temp_99,shape=[3,3,3,3]) do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) + &
temp_33 = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e),invFi)) crystallite_subdt(g,i,e)*math_mul33x33(temp_33,dLidFi(1:3,1:3,o,p))
lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) + &
crystallite_invFi(1:3,1:3,g,i,e)*crystallite_invFi(p,o,g,i,e)
rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) - &
crystallite_subdt(g,i,e)*math_mul33x33(temp_33,dLidS(1:3,1:3,o,p))
enddo; enddo
call math_invert(9_pInt,math_Plain3333to99(lhs_3333),temp_99,error)
if (error) then
call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=g, &
ext_msg='inversion error in analytic tangent calculation')
dFidS = 0.0_pReal
else
dFidS = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333)
endif
dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
call constitutive_LpAndItsTangent(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fi(1:3,1:3,g,i,e),g,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
temp_33 = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e), &
crystallite_invFi(1:3,1:3,g,i,e)))
rhs_3333 = 0.0_pReal rhs_3333 = 0.0_pReal
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & 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),temp_33) rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),temp_33)
@ -1152,29 +1208,25 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
math_inv33(crystallite_subFp0(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) & forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(temp_33,dLpdS(1:3,1:3,p,o)),invFi) temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(temp_33,dLpdS(1:3,1:3,p,o)), &
crystallite_invFi(1:3,1:3,g,i,e))
temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
crystallite_invFp(1:3,1:3,g,i,e)), invFi0) crystallite_invFp(1:3,1:3,g,i,e)), &
math_inv33(crystallite_subFi0(1:3,1:3,g,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & 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(temp_33,dLidS(1:3,1:3,p,o)) temp_3333(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + math_mul33x33(temp_33,dLidS(1:3,1:3,p,o))
lhs_3333 = crystallite_subdt(g,i,e)*math_mul3333xx3333(dSdFe,temp_3333) lhs_3333 = crystallite_subdt(g,i,e)*math_mul3333xx3333(dSdFe,temp_3333) + &
math_mul3333xx3333(dSdFi,dFidS)
temp_3333 = 0.0_pReal call math_invert(9_pInt,math_identity2nd(9_pInt)+math_Plain3333to99(lhs_3333),temp_99,error)
temp_33 = math_mul33x33(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) = math_mul33x33(temp_33,dLidS(1:3,1:3,p,o))
lhs_3333 = lhs_3333 + crystallite_subdt(g,i,e)*math_mul3333xx3333(dSdFiInv,temp_3333)
call math_invert(9_pInt,math_identity2nd(9_pInt)+reshape(lhs_3333,shape=[9,9]),temp_99,error)
if (error) then if (error) then
call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=g, & call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=g, &
ext_msg='inversion error in analytic tangent calculation') ext_msg='inversion error in analytic tangent calculation')
dSdF = rhs_3333 dSdF = rhs_3333
else else
dSdF = math_mul3333xx3333(reshape(temp_99,shape=[3,3,3,3]),rhs_3333) dSdF = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333)
endif endif
dFpinvdF = 0.0_pReal dFpinvdF = 0.0_pReal
@ -1182,30 +1234,24 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(g,i,e)* & 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(math_inv33(crystallite_subFp0(1:3,1:3,g,i,e)), &
math_mul33x33(temp_3333(1:3,1:3,p,o),invFi)) math_mul33x33(temp_3333(1:3,1:3,p,o), &
crystallite_invFi(1:3,1:3,g,i,e)))
dFiinvdF = 0.0_pReal
temp_3333 = math_mul3333xx3333(dLidS,dSdF)
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
dFiinvdF(1:3,1:3,p,o) = -crystallite_subdt(g,i,e)* &
math_mul33x33(math_inv33(crystallite_Fp(1:3,1:3,g,i,e)), &
math_mul33x33(invFi0,temp_3333(1:3,1:3,p,o)))
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.0_pReal crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.0_pReal
temp_33 = math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e), & temp_33 = math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e), &
math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,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))))*detFi math_transpose33(crystallite_invFp(1:3,1:3,g,i,e))))
forall(p=1_pInt:3_pInt) & forall(p=1_pInt:3_pInt) &
crystallite_dPdF(p,1:3,p,1:3,g,i,e) = math_transpose33(temp_33) crystallite_dPdF(p,1:3,p,1:3,g,i,e) = math_transpose33(temp_33)
temp_33 = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), & temp_33 = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), &
math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))*detFi math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & 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) + & 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)),temp_33) math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33)
temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
crystallite_invFp(1:3,1:3,g,i,e))*detFi crystallite_invFp(1:3,1:3,g,i,e))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & 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) + & 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(temp_33,dSdF(1:3,1:3,p,o)), & math_mul33x33(math_mul33x33(temp_33,dSdF(1:3,1:3,p,o)), &
@ -1213,15 +1259,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
crystallite_invFp(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)))*detFi math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & 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) + & crystallite_dPdF(1:3,1:3,p,o,g,i,e) = crystallite_dPdF(1:3,1:3,p,o,g,i,e) + &
math_mul33x33(temp_33,math_transpose33(dFpinvdF(1:3,1:3,p,o))) math_mul33x33(temp_33,math_transpose33(dFpinvdF(1:3,1:3,p,o)))
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) - &
crystallite_P(1:3,1:3,g,i,e)*sum(math_transpose33(Fi)*dFiinvdF(1:3,1:3,p,o))
enddo; enddo enddo; enddo
enddo elementLooping6 enddo elementLooping6
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -1259,8 +1301,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
F_backup(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ... and kinematics F_backup(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ... and kinematics
Fp_backup(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) Fp_backup(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e)
InvFp_backup(1:3,1:3,g,i,e) = crystallite_invFp(1:3,1:3,g,i,e) InvFp_backup(1:3,1:3,g,i,e) = crystallite_invFp(1:3,1:3,g,i,e)
Fi_backup(1:3,1:3,g,i,e) = crystallite_Fi(1:3,1:3,g,i,e)
InvFi_backup(1:3,1:3,g,i,e) = crystallite_invFi(1:3,1:3,g,i,e)
Fe_backup(1:3,1:3,g,i,e) = crystallite_Fe(1:3,1:3,g,i,e) Fe_backup(1:3,1:3,g,i,e) = crystallite_Fe(1:3,1:3,g,i,e)
Lp_backup(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) Lp_backup(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e)
Li_backup(1:3,1:3,g,i,e) = crystallite_Li(1:3,1:3,g,i,e)
Tstar_v_backup(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e) Tstar_v_backup(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e)
P_backup(1:3,1:3,g,i,e) = crystallite_P(1:3,1:3,g,i,e) P_backup(1:3,1:3,g,i,e) = crystallite_P(1:3,1:3,g,i,e)
convergenceFlag_backup(g,i,e) = crystallite_converged(g,i,e) convergenceFlag_backup(g,i,e) = crystallite_converged(g,i,e)
@ -1308,8 +1353,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_Fp(1:3,1:3,g,i,e) = Fp_backup(1:3,1:3,g,i,e) crystallite_Fp(1:3,1:3,g,i,e) = Fp_backup(1:3,1:3,g,i,e)
crystallite_invFp(1:3,1:3,g,i,e) = InvFp_backup(1:3,1:3,g,i,e) crystallite_invFp(1:3,1:3,g,i,e) = InvFp_backup(1:3,1:3,g,i,e)
crystallite_Fi(1:3,1:3,g,i,e) = Fi_backup(1:3,1:3,g,i,e)
crystallite_invFi(1:3,1:3,g,i,e) = InvFi_backup(1:3,1:3,g,i,e)
crystallite_Fe(1:3,1:3,g,i,e) = Fe_backup(1:3,1:3,g,i,e) crystallite_Fe(1:3,1:3,g,i,e) = Fe_backup(1:3,1:3,g,i,e)
crystallite_Lp(1:3,1:3,g,i,e) = Lp_backup(1:3,1:3,g,i,e) crystallite_Lp(1:3,1:3,g,i,e) = Lp_backup(1:3,1:3,g,i,e)
crystallite_Li(1:3,1:3,g,i,e) = Li_backup(1:3,1:3,g,i,e)
crystallite_Tstar_v(1:6,g,i,e) = Tstar_v_backup(1:6,g,i,e) crystallite_Tstar_v(1:6,g,i,e) = Tstar_v_backup(1:6,g,i,e)
endforall endforall
enddo enddo
@ -1339,8 +1387,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
vacancyState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) vacancyState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e)
crystallite_Fi(1:3,1:3,g,i,e) = crystallite_subFi0(1:3,1:3,g,i,e)
crystallite_Fe(1:3,1:3,g,i,e) = crystallite_subFe0(1:3,1:3,g,i,e) crystallite_Fe(1:3,1:3,g,i,e) = crystallite_subFe0(1:3,1:3,g,i,e)
crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e)
crystallite_Li(1:3,1:3,g,i,e) = crystallite_subLi0(1:3,1:3,g,i,e)
crystallite_Tstar_v(1:6,g,i,e) = crystallite_subTstar0_v(1:6,g,i,e) crystallite_Tstar_v(1:6,g,i,e) = crystallite_subTstar0_v(1:6,g,i,e)
endforall endforall
enddo enddo
@ -1445,8 +1495,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_subF(1:3,1:3,g,i,e) = F_backup(1:3,1:3,g,i,e) crystallite_subF(1:3,1:3,g,i,e) = F_backup(1:3,1:3,g,i,e)
crystallite_Fp(1:3,1:3,g,i,e) = Fp_backup(1:3,1:3,g,i,e) crystallite_Fp(1:3,1:3,g,i,e) = Fp_backup(1:3,1:3,g,i,e)
crystallite_invFp(1:3,1:3,g,i,e) = InvFp_backup(1:3,1:3,g,i,e) crystallite_invFp(1:3,1:3,g,i,e) = InvFp_backup(1:3,1:3,g,i,e)
crystallite_Fi(1:3,1:3,g,i,e) = Fi_backup(1:3,1:3,g,i,e)
crystallite_invFi(1:3,1:3,g,i,e) = InvFi_backup(1:3,1:3,g,i,e)
crystallite_Fe(1:3,1:3,g,i,e) = Fe_backup(1:3,1:3,g,i,e) crystallite_Fe(1:3,1:3,g,i,e) = Fe_backup(1:3,1:3,g,i,e)
crystallite_Lp(1:3,1:3,g,i,e) = Lp_backup(1:3,1:3,g,i,e) crystallite_Lp(1:3,1:3,g,i,e) = Lp_backup(1:3,1:3,g,i,e)
crystallite_Li(1:3,1:3,g,i,e) = Li_backup(1:3,1:3,g,i,e)
crystallite_Tstar_v(1:6,g,i,e) = Tstar_v_backup(1:6,g,i,e) crystallite_Tstar_v(1:6,g,i,e) = Tstar_v_backup(1:6,g,i,e)
crystallite_P(1:3,1:3,g,i,e) = P_backup(1:3,1:3,g,i,e) crystallite_P(1:3,1:3,g,i,e) = P_backup(1:3,1:3,g,i,e)
crystallite_converged(g,i,e) = convergenceFlag_backup(g,i,e) crystallite_converged(g,i,e) = convergenceFlag_backup(g,i,e)
@ -3552,11 +3605,10 @@ logical function crystallite_integrateStress(&
debug_g, & debug_g, &
debug_cumLpCalls, & debug_cumLpCalls, &
debug_cumLpTicks, & debug_cumLpTicks, &
debug_StressLoopDistribution debug_StressLoopLpDistribution, &
debug_StressLoopLiDistribution
use constitutive, only: constitutive_LpAndItsTangent, & use constitutive, only: constitutive_LpAndItsTangent, &
constitutive_LiAndItsTangent, & constitutive_LiAndItsTangent, &
constitutive_getFi0, &
constitutive_putFi, &
constitutive_TandItsTangent constitutive_TandItsTangent
use math, only: math_mul33x33, & use math, only: math_mul33x33, &
math_mul33xx33, & math_mul33xx33, &
@ -3576,7 +3628,8 @@ logical function crystallite_integrateStress(&
math_Mandel33to6, & math_Mandel33to6, &
math_Plain3333to99, & math_Plain3333to99, &
math_Plain33to9, & math_Plain33to9, &
math_Plain9to33 math_Plain9to33, &
math_Plain99to3333
use mesh, only: mesh_element use mesh, only: mesh_element
implicit none implicit none
@ -3592,6 +3645,8 @@ logical function crystallite_integrateStress(&
Fp_new, & ! plastic deformation gradient at end of timestep Fp_new, & ! plastic deformation gradient at end of timestep
Fe_new, & ! elastic deformation gradient at end of timestep Fe_new, & ! elastic deformation gradient at end of timestep
invFp_new, & ! inverse of Fp_new invFp_new, & ! inverse of Fp_new
Fi_new, & ! gradient of intermediate deformation stages
invFi_new, &
invFp_current, & ! inverse of Fp_current invFp_current, & ! inverse of Fp_current
invFi_current, & ! inverse of Fp_current invFi_current, & ! inverse of Fp_current
Lpguess, & ! current guess for plastic velocity gradient Lpguess, & ! current guess for plastic velocity gradient
@ -3607,32 +3662,26 @@ logical function crystallite_integrateStress(&
residuumLi_old, & ! last residuum of intermediate velocity gradient residuumLi_old, & ! last residuum of intermediate velocity gradient
deltaLi, & ! direction of next guess deltaLi, & ! direction of next guess
Tstar, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration Tstar, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration
Tstar_unloaded, & ! 2nd Piola-Kirchhoff Stress in unloaded configuration
A, & A, &
B, & B, &
Fe, & ! elastic deformation gradient Fe, & ! elastic deformation gradient
Fi, & ! gradient of intermediate deformation stages
invFi, &
temp_33 temp_33
real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation
real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK
integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK
real(pReal), dimension(9,9) :: dLp_dT_constitutive99, & ! partial derivative of plastic velocity gradient calculated by constitutive law real(pReal), dimension(9,9) :: dLp_dT_constitutive99, & ! partial derivative of plastic velocity gradient calculated by constitutive law
dLi_dT_constitutive99, & ! partial derivative of intermediate velocity gradient calculated by constitutive law
dT_dFe99, & ! partial derivative of 2nd Piola-Kirchhoff stress calculated by constitutive law
dFe_dLp99, & ! partial derivative of elastic deformation gradient
dFe_dLi99, &
dFiInv_dLi99, &
dT_dFiInv99, &
dRLp_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) dRLp_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme)
dRLp_dLp2, & ! working copy of dRdLp dRLp_dLp2, & ! working copy of dRdLp
dRLi_dLi ! partial derivative of residuumI (Jacobian for NEwton-Raphson scheme) dRLi_dLi ! partial derivative of residuumI (Jacobian for NEwton-Raphson scheme)
real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress
dT_dFe3333_unloaded, & dT_dFi3333, &
dFe_dLp3333, & ! partial derivative of elastic deformation gradient dFe_dLp3333, & ! partial derivative of elastic deformation gradient
dFe_dLi3333, & dFe_dLi3333, &
dFiInv_dLi3333, & dFi_dLi3333, &
dT_dFiInv3333 dLp_dFi3333, &
dLi_dFi3333, &
dLp_dT3333, &
dLi_dT3333
real(pReal) det, & ! determinant real(pReal) det, & ! determinant
detInvFi, & detInvFi, &
steplengthLp0, & steplengthLp0, &
@ -3687,8 +3736,11 @@ logical function crystallite_integrateStress(&
!* feed local variables !* feed local variables
Fp_current = crystallite_subFp0(1:3,1:3,g,i,e) ! "Fp_current" is only used as temp var here... Fp_current = crystallite_subFp0(1:3,1:3,g,i,e) ! "Fp_current" is only used as temp var here...
Lpguess = crystallite_Lp (1:3,1:3,g,i,e) ! ... and take it as first guess Lpguess = crystallite_Lp (1:3,1:3,g,i,e) ! ... and take it as first guess
Fi_current = crystallite_subFi0(1:3,1:3,g,i,e) ! intermediate configuration, assume decomposition as F = Fe Fi Fp
Liguess = crystallite_Li (1:3,1:3,g,i,e) ! ... and take it as first guess
Liguess_old = Liguess
!* inversion of Fp_current... !* inversion of Fp_current...
@ -3705,19 +3757,12 @@ logical function crystallite_integrateStress(&
#endif #endif
return return
endif endif
A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
!* feed local variables
Fi_current = constitutive_getFi0(g,i,e) ! intermediate configuration, assume decomposition as F = Fe Fi Fp
call constitutive_LiAndItsTangent(Liguess, dLi_dT_constitutive99, &
crystallite_Tstar_v(1:6,g,i,e), Lpguess, g, i, e)
Liguess_old = Liguess
!* inversion of Fi_current... !* inversion of Fi_current...
invFi_current = math_inv33(Fi_current) invFi_current = math_inv33(Fi_current)
if (all(invFi_current == 0.0_pReal)) then ! ... failed? if (all(invFi_current == 0.0_pReal)) then ! ... failed?
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then 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 ',& 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 ',&
@ -3748,9 +3793,9 @@ logical function crystallite_integrateStress(&
return return
endif IloopsExeced endif IloopsExeced
invFi = math_mul33x33(invFi_current,math_I3 - dt*Liguess) invFi_new = math_mul33x33(invFi_current,math_I3 - dt*Liguess)
detInvFi = math_det33(invFi) Fi_new = math_inv33(invFi_new)
Fi = math_inv33(invFi) detInvFi = math_det33(invFi_new)
NiterationStressLp = 0_pInt NiterationStressLp = 0_pInt
jacoCounterLp = 0_pInt jacoCounterLp = 0_pInt
@ -3773,32 +3818,18 @@ logical function crystallite_integrateStress(&
!* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
B = math_I3 - dt*Lpguess B = math_I3 - dt*Lpguess
Fe = math_mul33x33(math_mul33x33(A,B),invFi) ! current elastic deformation tensor Fe = math_mul33x33(math_mul33x33(A,B), invFi_new) ! current elastic deformation tensor
call constitutive_TandItsTangent(Tstar_unloaded, dT_dFe3333_unloaded, Fe, g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration call constitutive_TandItsTangent(Tstar, dT_dFe3333, dT_dFi3333, Fe, Fi_new, 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
dT_dFe3333 = 0.0_pReal
dT_dFiInv3333 = 0.0_pReal
temp_33 = math_mul33x33(Tstar_unloaded,math_transpose33(invFi))/detInvFi ! push Tstar forward from unloaded to plastic (lattice) configuration
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
dT_dFe3333 (1:3,1:3,o,p) = &
math_mul33x33(invFi,math_mul33x33(dT_dFe3333_unloaded(1:3,1:3,o,p), &
math_transpose33(invFi)))/detInvFi
dT_dFiInv3333(o,1:3,p,1:3) = dT_dFiInv3333(o,1:3,p,1:3) + math_I3(o,p)*temp_33
dT_dFiInv3333(1:3,o,p,1:3) = dT_dFiInv3333(1:3,o,p,1:3) + math_I3(o,p)*math_transpose33(temp_33)
dT_dFiInv3333(1:3,1:3,o,p) = dT_dFiInv3333(1:3,1:3,o,p) - Tstar*invFi(p,o)
enddo; enddo
Tstar_v = math_Mandel33to6(Tstar) Tstar_v = math_Mandel33to6(Tstar)
!* calculate plastic velocity gradient and its tangent from constitutive law !* calculate plastic velocity gradient and its tangent from constitutive law
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
endif endif
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive99, Tstar_v, & call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT3333, dLp_dFi3333, &
g, i, e) Tstar_v, Fi_new, g, i, e)
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
@ -3855,13 +3886,11 @@ logical function crystallite_integrateStress(&
if (mod(jacoCounterLp, iJacoLpresiduum) == 0_pInt) then if (mod(jacoCounterLp, iJacoLpresiduum) == 0_pInt) then
dFe_dLp3333 = 0.0_pReal dFe_dLp3333 = 0.0_pReal
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
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) dFe_dLp3333(o,1:3,p,1:3) = A(o,p)*math_transpose33(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
enddo; enddo enddo; enddo
dFe_dLp3333 = - dt * dFe_dLp3333 dFe_dLp3333 = - dt * dFe_dLp3333
dFe_dLp99 = math_Plain3333to99(dFe_dLp3333)
dT_dFe99 = math_Plain3333to99(dT_dFe3333)
dRLp_dLp = math_identity2nd(9_pInt) & dRLp_dLp = math_identity2nd(9_pInt) &
- math_mul99x99(dLp_dT_constitutive99, math_mul99x99(dT_dFe99 , dFe_dLp99)) - math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333))
dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine
work = math_plain33to9(residuumLp) work = math_plain33to9(residuumLp)
#if(FLOAT==8) #if(FLOAT==8)
@ -3879,8 +3908,8 @@ logical function crystallite_integrateStress(&
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,*) write(6,*)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLp',transpose(dRLp_dLp) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLp',transpose(dRLp_dLp)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLp',transpose(dFe_dLp99) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLp',transpose(math_Plain3333to99(dFe_dLp3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFe_constitutive',transpose(dT_dFe99) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFe_constitutive',transpose(math_Plain3333to99(dT_dFe3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(dLp_dT_constitutive99) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(dLp_dT_constitutive99)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> A',math_transpose33(A) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> A',math_transpose33(A)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> B',math_transpose33(B) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> B',math_transpose33(B)
@ -3899,10 +3928,17 @@ logical function crystallite_integrateStress(&
enddo LpLoop enddo LpLoop
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionStress)
debug_StressLoopLpDistribution(NiterationStressLp,numerics_integrationMode) = &
debug_StressLoopLpDistribution(NiterationStressLp,numerics_integrationMode) + 1_pInt
!$OMP END CRITICAL (distributionStress)
endif
!* calculate intermediate velocity gradient and its tangent from constitutive law !* calculate intermediate velocity gradient and its tangent from constitutive law
call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT_constitutive99, Tstar_v, Lpguess, & call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT3333, dLi_dFi3333, &
g, i, e) Tstar_v, Fi_new, Lpguess, g, i, e)
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
@ -3928,18 +3964,20 @@ logical function crystallite_integrateStress(&
if (mod(jacoCounterLi, iJacoLpresiduum) == 0_pInt) then if (mod(jacoCounterLi, iJacoLpresiduum) == 0_pInt) then
temp_33 = math_mul33x33(math_mul33x33(A,B),invFi_current) temp_33 = math_mul33x33(math_mul33x33(A,B),invFi_current)
dFe_dLi3333 = 0.0_pReal dFe_dLi3333 = 0.0_pReal
dFiInv_dLi3333 = 0.0_pReal dFi_dLi3333 = 0.0_pReal
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
dFe_dLi3333 (1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) dFe_dLi3333(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
dFiInv_dLi3333(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current dFi_dLi3333(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current
enddo; enddo enddo; enddo
dT_dFiInv99 = math_Plain3333to99(dT_dFiInv3333) do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
dFe_dLi99 = math_Plain3333to99(dFe_dLi3333) dFi_dLi3333(1:3,1:3,o,p) = math_mul33x33(math_mul33x33(Fi_new,dFi_dLi3333(1:3,1:3,o,p)),Fi_new)
dFiInv_dLi99 = math_Plain3333to99(dFiInv_dLi3333) enddo; enddo
dRLi_dLi = math_identity2nd(9_pInt) &
- math_mul99x99(dLi_dT_constitutive99, math_mul99x99(dT_dFe99, dFe_dLi99) + & dRLi_dLi = math_identity2nd(9_pInt) &
math_mul99x99(dT_dFiInv99, dFiInv_dLi99)) - math_Plain3333to99(math_mul3333xx3333(dLi_dT3333, math_mul3333xx3333(dT_dFe3333, dFe_dLi3333) + &
math_mul3333xx3333(dT_dFi3333, dFi_dLi3333))) &
- math_Plain3333to99(math_mul3333xx3333(dLi_dFi3333, dFi_dLi3333))
work = math_plain33to9(residuumLi) work = math_plain33to9(residuumLi)
#if(FLOAT==8) #if(FLOAT==8)
call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
@ -3956,6 +3994,13 @@ logical function crystallite_integrateStress(&
Liguess = Liguess + steplengthLi * deltaLi Liguess = Liguess + steplengthLi * deltaLi
enddo LiLoop enddo LiLoop
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionStress)
debug_StressLoopLiDistribution(NiterationStressLi,numerics_integrationMode) = &
debug_StressLoopLiDistribution(NiterationStressLi,numerics_integrationMode) + 1_pInt
!$OMP END CRITICAL (distributionStress)
endif
!* calculate new plastic and elastic deformation gradient !* calculate new plastic and elastic deformation gradient
invFp_new = math_mul33x33(invFp_current,B) invFp_new = math_mul33x33(invFp_current,B)
@ -3974,22 +4019,24 @@ logical function crystallite_integrateStress(&
#endif #endif
return return
endif endif
Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi) ! calc resulting Fe Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi_new) ! calc resulting Fe
!* calculate 1st Piola-Kirchhoff stress !* calculate 1st Piola-Kirchhoff stress
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(Fg_new,invFp_new), & crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(Fg_new,invFp_new), &
math_mul33x33(math_Mandel6to33(Tstar_v), & math_mul33x33(math_Mandel6to33(Tstar_v), &
math_transpose33(invFp_new)))/detInvFi math_transpose33(invFp_new)))
!* store local values in global variables !* store local values in global variables
crystallite_Lp(1:3,1:3,g,i,e) = Lpguess crystallite_Lp(1:3,1:3,g,i,e) = Lpguess
crystallite_Li(1:3,1:3,g,i,e) = Liguess
crystallite_Tstar_v(1:6,g,i,e) = Tstar_v crystallite_Tstar_v(1:6,g,i,e) = Tstar_v
crystallite_Fp(1:3,1:3,g,i,e) = Fp_new crystallite_Fp(1:3,1:3,g,i,e) = Fp_new
crystallite_Fi(1:3,1:3,g,i,e) = Fi_new
crystallite_Fe(1:3,1:3,g,i,e) = Fe_new crystallite_Fe(1:3,1:3,g,i,e) = Fe_new
crystallite_invFp(1:3,1:3,g,i,e) = invFp_new crystallite_invFp(1:3,1:3,g,i,e) = invFp_new
call constitutive_putFi(Tstar_v, Lpguess, dt, g, i, e) crystallite_invFi(1:3,1:3,g,i,e) = invFi_new
!* set return flag to true !* set return flag to true
@ -4004,16 +4051,10 @@ logical function crystallite_integrateStress(&
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fe Lp Fe^-1', & write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fe Lp Fe^-1', &
math_transpose33(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv33(Fe_new)))) ! transpose to get correct print out order math_transpose33(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv33(Fe_new)))) ! transpose to get correct print out order
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp',math_transpose33(crystallite_Fp(1:3,1:3,g,i,e)) write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp',math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fi',math_transpose33(crystallite_Fi(1:3,1:3,g,i,e))
endif endif
#endif #endif
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionStress)
debug_StressLoopDistribution(NiterationStressLp,numerics_integrationMode) = &
debug_StressLoopDistribution(NiterationStressLp,numerics_integrationMode) + 1_pInt
!$OMP END CRITICAL (distributionStress)
endif
end function crystallite_integrateStress end function crystallite_integrateStress
@ -4275,10 +4316,18 @@ function crystallite_postResults(ipc, ip, el)
mySize = 9_pInt mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = & crystallite_postResults(c+1:c+mySize) = &
reshape(math_transpose33(crystallite_Fp(1:3,1:3,ipc,ip,el)),[mySize]) reshape(math_transpose33(crystallite_Fp(1:3,1:3,ipc,ip,el)),[mySize])
case (fi_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &
reshape(math_transpose33(crystallite_Fi(1:3,1:3,ipc,ip,el)),[mySize])
case (lp_ID) case (lp_ID)
mySize = 9_pInt mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = & crystallite_postResults(c+1:c+mySize) = &
reshape(math_transpose33(crystallite_Lp(1:3,1:3,ipc,ip,el)),[mySize]) reshape(math_transpose33(crystallite_Lp(1:3,1:3,ipc,ip,el)),[mySize])
case (li_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &
reshape(math_transpose33(crystallite_Li(1:3,1:3,ipc,ip,el)),[mySize])
case (p_ID) case (p_ID)
mySize = 9_pInt mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = & crystallite_postResults(c+1:c+mySize) = &

268
code/damage_anisoBrittle.f90
View File

@ -56,13 +56,10 @@ module damage_anisoBrittle
damage_anisoBrittle_aTolState, & damage_anisoBrittle_aTolState, &
damage_anisoBrittle_microstructure, & damage_anisoBrittle_microstructure, &
damage_anisoBrittle_LdAndItsTangent, & damage_anisoBrittle_LdAndItsTangent, &
damage_anisoBrittle_getFd, &
damage_anisoBrittle_putFd, &
damage_anisoBrittle_getFd0, &
damage_anisoBrittle_getPartionedFd0, &
damage_anisoBrittle_getDamage, & damage_anisoBrittle_getDamage, &
damage_anisoBrittle_putLocalDamage, & damage_anisoBrittle_putLocalDamage, &
damage_anisoBrittle_getLocalDamage, & damage_anisoBrittle_getLocalDamage, &
damage_anisoBrittle_getDamageDiffusion33, &
damage_anisoBrittle_postResults damage_anisoBrittle_postResults
contains contains
@ -133,20 +130,20 @@ subroutine damage_anisoBrittle_init(fileUnit)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(damage_anisoBrittle_sizePostResults(maxNinstance), source=0_pInt) allocate(damage_anisoBrittle_sizePostResults(maxNinstance), source=0_pInt)
allocate(damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) allocate(damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),maxNinstance), source=0_pInt)
allocate(damage_anisoBrittle_output(maxval(phase_Noutput),maxNinstance)) allocate(damage_anisoBrittle_output(maxval(phase_Noutput),maxNinstance))
damage_anisoBrittle_output = '' damage_anisoBrittle_output = ''
allocate(damage_anisoBrittle_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) allocate(damage_anisoBrittle_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
allocate(damage_anisoBrittle_Noutput(maxNinstance), source=0_pInt) allocate(damage_anisoBrittle_Noutput(maxNinstance), source=0_pInt)
allocate(damage_anisoBrittle_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) 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_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
allocate(damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,maxNinstance),source=0_pInt) allocate(damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt)
allocate(damage_anisoBrittle_totalNcleavage(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_damage(maxNinstance), source=0.0_pReal)
allocate(damage_anisoBrittle_aTol_disp(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_sdot_0(maxNinstance), source=0.0_pReal)
allocate(damage_anisoBrittle_N(maxNinstance), source=0.0_pReal) allocate(damage_anisoBrittle_N(maxNinstance), source=0.0_pReal)
rewind(fileUnit) rewind(fileUnit)
phase = 0_pInt phase = 0_pInt
@ -226,9 +223,9 @@ subroutine damage_anisoBrittle_init(fileUnit)
damage_anisoBrittle_aTol_disp(instance) = 1.0e-3_pReal ! default absolute tolerance 1e-3 damage_anisoBrittle_aTol_disp(instance) = 1.0e-3_pReal ! default absolute tolerance 1e-3
if (damage_anisoBrittle_sdot_0(instance) <= 0.0_pReal) & if (damage_anisoBrittle_sdot_0(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//LOCAL_DAMAGE_anisoBrittle_LABEL//')') call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//LOCAL_DAMAGE_anisoBrittle_LABEL//')')
if (any(damage_anisoBrittle_critDisp(:,instance) < 0.0_pReal)) & if (any(damage_anisoBrittle_critDisp(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='critical_displacement ('//LOCAL_DAMAGE_anisoBrittle_LABEL//')') call IO_error(211_pInt,el=instance,ext_msg='critical_displacement ('//LOCAL_DAMAGE_anisoBrittle_LABEL//')')
if (any(damage_anisoBrittle_critLoad(:,instance) < 0.0_pReal)) & if (any(damage_anisoBrittle_critLoad(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='critical_load ('//LOCAL_DAMAGE_anisoBrittle_LABEL//')') call IO_error(211_pInt,el=instance,ext_msg='critical_load ('//LOCAL_DAMAGE_anisoBrittle_LABEL//')')
if (damage_anisoBrittle_N(instance) <= 0.0_pReal) & if (damage_anisoBrittle_N(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity_damage ('//LOCAL_DAMAGE_anisoBrittle_LABEL//')') call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity_damage ('//LOCAL_DAMAGE_anisoBrittle_LABEL//')')
@ -256,9 +253,8 @@ subroutine damage_anisoBrittle_init(fileUnit)
! Determine size of state array ! Determine size of state array
sizeDotState = 0_pInt sizeDotState = 0_pInt
sizeState = sizeDotState + & sizeState = sizeDotState + &
1_pInt + & ! non-local damage 1_pInt + & ! non-local damage
1_pInt + & ! opening on each damage system 1_pInt ! opening on each damage system
9_pInt ! Fd
damageState(phase)%sizeState = sizeState damageState(phase)%sizeState = sizeState
damageState(phase)%sizeDotState = sizeDotState damageState(phase)%sizeDotState = sizeDotState
@ -304,8 +300,7 @@ subroutine damage_anisoBrittle_stateInit(phase,instance)
real(pReal), dimension(damageState(phase)%sizeState) :: tempState real(pReal), dimension(damageState(phase)%sizeState) :: tempState
tempState(1) = 1.0_pReal tempState(1) = 1.0_pReal
tempState(2) = 0.0_pReal tempState(2) = 1.0_pReal
tempState(3:11) = reshape(math_I3, shape=[9])
damageState(phase)%state0 = spread(tempState,2,size(damageState(phase)%state(1,:))) damageState(phase)%state0 = spread(tempState,2,size(damageState(phase)%state(1,:)))
@ -326,7 +321,6 @@ subroutine damage_anisoBrittle_aTolState(phase,instance)
tempTol(1) = damage_anisoBrittle_aTol_damage(instance) tempTol(1) = damage_anisoBrittle_aTol_damage(instance)
tempTol(2) = damage_anisoBrittle_aTol_disp (instance) tempTol(2) = damage_anisoBrittle_aTol_disp (instance)
tempTol(3:11) = damage_anisoBrittle_aTol_damage(instance)
damageState(phase)%aTolState = tempTol damageState(phase)%aTolState = tempTol
@ -385,27 +379,25 @@ subroutine damage_anisoBrittle_microstructure(Tstar_v, subdt, ipc, ip, el)
damageState(phase)%state(2,constituent) = & damageState(phase)%state(2,constituent) = &
damageState(phase)%state(2,constituent) + & damageState(phase)%state(2,constituent) + &
subdt*damage_anisoBrittle_sdot_0(instance)* & subdt*damage_anisoBrittle_sdot_0(instance)* &
( (abs(traction_d)/traction_crit)**damage_anisoBrittle_N(instance) + & ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**damage_anisoBrittle_N(instance) + &
(abs(traction_t)/traction_crit)**damage_anisoBrittle_N(instance) + & (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**damage_anisoBrittle_N(instance) + &
(max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance))/ & (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**damage_anisoBrittle_N(instance))/ &
damage_anisoBrittle_critDisp(f,instance) damage_anisoBrittle_critDisp(f,instance)
enddo enddo
enddo enddo
localDamage = min(1.0_pReal,max(residualStiffness,1.0/damageState(phase)%state(2,constituent))) localDamage = max(residualStiffness,1.0_pReal/damageState(phase)%state(2,constituent))
damageState(phase)%state(1,constituent) = & damageState(phase)%state(1,constituent) = &
localDamage + & localDamage
(damageState(phase)%subState0(1,constituent) - localDamage)* &
exp(-subdt/lattice_DamageMobility(phase))
end subroutine damage_anisoBrittle_microstructure end subroutine damage_anisoBrittle_microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el) subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar3333, Tstar_v, ipc, ip, el)
use material, only: & use material, only: &
mappingConstitutive, & mappingConstitutive, &
phase_damageInstance, & phase_damageInstance, &
@ -427,9 +419,7 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
Tstar_v !< 2nd Piola-Kirchhoff stress Tstar_v !< 2nd Piola-Kirchhoff stress
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient Ld !< damage velocity gradient
real(pReal), intent(out), dimension(9,9) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
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) dLd_dTstar3333 !< derivative of Ld with respect to Tstar (4th-order tensor)
integer(pInt) :: & integer(pInt) :: &
phase, & phase, &
@ -457,10 +447,11 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
udotd = & udotd = &
sign(1.0_pReal,traction_d)* & sign(1.0_pReal,traction_d)* &
damage_anisoBrittle_sdot_0(instance)* & damage_anisoBrittle_sdot_0(instance)* &
(abs(traction_d)/traction_crit)**damage_anisoBrittle_N(instance) (max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**damage_anisoBrittle_N(instance)
if (udotd /= 0.0_pReal) then if (udotd /= 0.0_pReal) then
Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase) Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)
dudotd_dt = udotd*damage_anisoBrittle_N(instance)/traction_d dudotd_dt = sign(1.0_pReal,traction_d)*udotd*damage_anisoBrittle_N(instance)/ &
max(0.0_pReal, abs(traction_d) - traction_crit)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & 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) + & dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* & dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* &
@ -470,10 +461,11 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
udott = & udott = &
sign(1.0_pReal,traction_t)* & sign(1.0_pReal,traction_t)* &
damage_anisoBrittle_sdot_0(instance)* & damage_anisoBrittle_sdot_0(instance)* &
(abs(traction_t)/traction_crit)**damage_anisoBrittle_N(instance) (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**damage_anisoBrittle_N(instance)
if (udott /= 0.0_pReal) then if (udott /= 0.0_pReal) then
Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase) Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)
dudott_dt = udott*damage_anisoBrittle_N(instance)/traction_t dudott_dt = sign(1.0_pReal,traction_t)*udott*damage_anisoBrittle_N(instance)/ &
max(0.0_pReal, abs(traction_t) - traction_crit)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & 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) + & dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* & dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* &
@ -481,11 +473,13 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
endif endif
udotn = & udotn = &
sign(1.0_pReal,traction_n)* &
damage_anisoBrittle_sdot_0(instance)* & damage_anisoBrittle_sdot_0(instance)* &
(max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance) (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**damage_anisoBrittle_N(instance)
if (udotn /= 0.0_pReal) then if (udotn /= 0.0_pReal) then
Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase) Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
dudotn_dt = udotn*damage_anisoBrittle_N(instance)/traction_n dudotn_dt = sign(1.0_pReal,traction_n)*udotn*damage_anisoBrittle_N(instance)/ &
max(0.0_pReal, abs(traction_n) - traction_crit)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & 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) + & dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* & dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* &
@ -494,174 +488,9 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
enddo enddo
enddo enddo
dLd_dTstar = math_Plain3333to99(dLd_dTstar3333)
end subroutine damage_anisoBrittle_LdAndItsTangent 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:11,constituent),shape=[3,3])
end function damage_anisoBrittle_getFd
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el)
use material, only: &
mappingConstitutive, &
phase_damageInstance, &
damageState
use math, only: &
math_mul33x33, &
math_inv33, &
math_I3
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(in) :: &
dt
real(pReal), dimension(3,3) :: &
Ld !< damage velocity gradient
integer(pInt) :: &
phase, &
constituent, &
instance, &
f, i, index_myFamily
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, udott, udotn
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase)
Ld = 0.0_pReal
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
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 = damage_anisoBrittle_critLoad(f,instance)* &
damage_anisoBrittle_getDamage(ipc, ip, el)
udotd = &
sign(1.0_pReal,traction_d)* &
damage_anisoBrittle_sdot_0(instance)* &
(abs(traction_d)/traction_crit)**damage_anisoBrittle_N(instance)
Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)
udott = &
sign(1.0_pReal,traction_t)* &
damage_anisoBrittle_sdot_0(instance)* &
(abs(traction_t)/traction_crit)**damage_anisoBrittle_N(instance)
Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)
udotn = &
damage_anisoBrittle_sdot_0(instance)* &
(max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance)
Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
enddo
enddo
damageState(phase)%state(3:11,constituent) = &
reshape(math_mul33x33(math_inv33(math_I3 - dt*Ld), &
damage_anisoBrittle_getFd0(ipc, ip, el)), shape=[9])
end subroutine damage_anisoBrittle_putFd
!--------------------------------------------------------------------------------------------------
!> @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: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:11,constituent),shape=[3,3])
end function damage_anisoBrittle_getPartionedFd0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns damage !> @brief returns damage
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -738,6 +567,33 @@ function damage_anisoBrittle_getLocalDamage(ipc, ip, el)
end function damage_anisoBrittle_getLocalDamage end function damage_anisoBrittle_getLocalDamage
!--------------------------------------------------------------------------------------------------
!> @brief returns brittle damage diffusion tensor
!--------------------------------------------------------------------------------------------------
function damage_anisoBrittle_getDamageDiffusion33(ipc, ip, el)
use lattice, only: &
lattice_DamageDiffusion33
use material, only: &
mappingConstitutive, &
damageState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
damage_anisoBrittle_getDamageDiffusion33
integer(pInt) :: &
phase, constituent
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
damage_anisoBrittle_getDamageDiffusion33 = &
lattice_DamageDiffusion33(1:3,1:3,phase)
end function damage_anisoBrittle_getDamageDiffusion33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results !> @brief return array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

226
code/damage_anisoDuctile.f90
View File

@ -59,10 +59,6 @@ module damage_anisoDuctile
damage_anisoDuctile_aTolState, & damage_anisoDuctile_aTolState, &
damage_anisoDuctile_microstructure, & damage_anisoDuctile_microstructure, &
damage_anisoDuctile_LdAndItsTangent, & damage_anisoDuctile_LdAndItsTangent, &
damage_anisoDuctile_getFd, &
damage_anisoDuctile_putFd, &
damage_anisoDuctile_getFd0, &
damage_anisoDuctile_getPartionedFd0, &
damage_anisoDuctile_getDamage, & damage_anisoDuctile_getDamage, &
damage_anisoDuctile_putLocalDamage, & damage_anisoDuctile_putLocalDamage, &
damage_anisoDuctile_getLocalDamage, & damage_anisoDuctile_getLocalDamage, &
@ -253,8 +249,7 @@ subroutine damage_anisoDuctile_init(fileUnit)
sizeDotState = 0_pInt sizeDotState = 0_pInt
sizeState = sizeDotState + & sizeState = sizeDotState + &
1_pInt + & ! time regularised damage 1_pInt + & ! time regularised damage
1_pInt + & ! damaged plasticity 1_pInt ! damaged plasticity
9 ! Fd
damageState(phase)%sizeState = sizeState damageState(phase)%sizeState = sizeState
damageState(phase)%sizeDotState = sizeDotState damageState(phase)%sizeDotState = sizeDotState
damageState(phase)%sizePostResults = damage_anisoDuctile_sizePostResults(instance) damageState(phase)%sizePostResults = damage_anisoDuctile_sizePostResults(instance)
@ -298,7 +293,6 @@ subroutine damage_anisoDuctile_stateInit(phase, instance)
tempState(1) = 1.0_pReal tempState(1) = 1.0_pReal
tempState(2) = 0.0_pReal tempState(2) = 0.0_pReal
tempState(3:11) = reshape(math_I3, shape=[9])
damageState(phase)%state0 = spread(tempState,2,size(damageState(phase)%state(1,:))) damageState(phase)%state0 = spread(tempState,2,size(damageState(phase)%state(1,:)))
@ -384,7 +378,7 @@ end subroutine damage_anisoDuctile_microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el) subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar3333, Tstar_v, ipc, ip, el)
use lattice, only: & use lattice, only: &
lattice_maxNslipFamily, & lattice_maxNslipFamily, &
lattice_NslipSystem, & lattice_NslipSystem, &
@ -401,7 +395,10 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
math_identity4th, & math_identity4th, &
math_symmetric33, & math_symmetric33, &
math_Mandel33to6, & math_Mandel33to6, &
math_tensorproduct math_tensorproduct, &
math_det33, &
math_mul33x33, &
math_transpose33
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -412,9 +409,7 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
Tstar_v !< 2nd Piola-Kirchhoff stress Tstar_v !< 2nd Piola-Kirchhoff stress
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient Ld !< damage velocity gradient
real(pReal), intent(out), dimension(9,9) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
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) dLd_dTstar3333 !< derivative of Ld with respect to Tstar (4th-order tensor)
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
projection_d, projection_t, projection_n !< projection modes 3x3 tensor projection_d, projection_t, projection_n !< projection modes 3x3 tensor
@ -466,7 +461,7 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
dudotd_dt = udotd*damage_anisoDuctile_N(instance)/traction_d dudotd_dt = udotd*damage_anisoDuctile_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) & 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) + & dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudotd_dt*projection_d(k,l)* projection_d(m,n) dudotd_dt*projection_d(k,l)*projection_d(m,n)
endif endif
udott = & udott = &
@ -486,216 +481,17 @@ subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
(max(0.0_pReal,traction_n)/traction_crit - & (max(0.0_pReal,traction_n)/traction_crit - &
max(0.0_pReal,traction_n)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance) max(0.0_pReal,traction_n)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance)
if (udotn /= 0.0_pReal) then if (udotn /= 0.0_pReal) then
Ld = Ld + udotn*projection_n(1:3,1:3) Ld = Ld + udotn*projection_n
dudotn_dt = udotn*damage_anisoDuctile_N(instance)/traction_n dudotn_dt = udotn*damage_anisoDuctile_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) & 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) + & dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudotn_dt*projection_n(k,l)* projection_n(m,n) dudotn_dt*projection_n(k,l)*projection_n(m,n)
endif endif
enddo enddo
enddo enddo
dLd_dTstar = math_Plain3333to99(dLd_dTstar3333)
end subroutine damage_anisoDuctile_LdAndItsTangent end subroutine damage_anisoDuctile_LdAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief returns local damage deformation gradient
!--------------------------------------------------------------------------------------------------
pure function damage_anisoDuctile_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_anisoDuctile_getFd
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase)
damage_anisoDuctile_getFd = &
reshape(damageState(phase)%state(3:11,constituent),shape=[3,3])
end function damage_anisoDuctile_getFd
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el)
use material, only: &
mappingConstitutive, &
phase_damageInstance, &
damageState
use math, only: &
math_mul33x33, &
math_inv33, &
math_I3, &
math_symmetric33, &
math_Mandel33to6, &
math_tensorproduct
use lattice, only: &
lattice_maxNslipFamily, &
lattice_NslipSystem, &
lattice_sd, &
lattice_st, &
lattice_sn
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) :: &
dt
real(pReal), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), dimension(3,3) :: &
projection_d, projection_t, projection_n !< projection modes 3x3 tensor
real(pReal), dimension(6) :: &
projection_d_v, projection_t_v, projection_n_v !< projection modes 3x3 vector
integer(pInt) :: &
phase, &
constituent, &
instance, &
f, i, index_myFamily
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, udott, udotn
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase)
Ld = 0.0_pReal
do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,damage_anisoDuctile_Nslip(f,instance) ! process each (active) slip system in family
projection_d = math_tensorproduct(lattice_sd(1:3,index_myFamily+i,phase),&
lattice_sn(1:3,index_myFamily+i,phase))
projection_t = math_tensorproduct(lattice_st(1:3,index_myFamily+i,phase),&
lattice_sn(1:3,index_myFamily+i,phase))
projection_n = math_tensorproduct(lattice_sn(1:3,index_myFamily+i,phase),&
lattice_sn(1:3,index_myFamily+i,phase))
projection_d_v(1:6) = math_Mandel33to6(math_symmetric33(projection_d(1:3,1:3)))
projection_t_v(1:6) = math_Mandel33to6(math_symmetric33(projection_t(1:3,1:3)))
projection_n_v(1:6) = math_Mandel33to6(math_symmetric33(projection_n(1:3,1:3)))
traction_d = dot_product(Tstar_v,projection_d_v(1:6))
traction_t = dot_product(Tstar_v,projection_t_v(1:6))
traction_n = dot_product(Tstar_v,projection_n_v(1:6))
traction_crit = damage_anisoDuctile_critLoad(f,instance)* &
damage_anisoDuctile_getDamage(ipc, ip, el) ! degrading critical load carrying capacity by damage
udotd = &
sign(1.0_pReal,traction_d)* &
damage_anisoDuctile_sdot_0(instance)* &
(abs(traction_d)/traction_crit - &
abs(traction_d)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance)
if (udotd /= 0.0_pReal) then
Ld = Ld + udotd*projection_d
endif
udott = &
sign(1.0_pReal,traction_t)* &
damage_anisoDuctile_sdot_0(instance)* &
(abs(traction_t)/traction_crit) - &
abs(traction_t)/damage_anisoDuctile_critLoad(f,instance)**damage_anisoDuctile_N(instance)
if (udott /= 0.0_pReal) then
Ld = Ld + udott*projection_t
endif
udotn = &
damage_anisoDuctile_sdot_0(instance)* &
(max(0.0_pReal,traction_n)/traction_crit - &
max(0.0_pReal,traction_n)/damage_anisoDuctile_critLoad(f,instance))**damage_anisoDuctile_N(instance)
if (udotn /= 0.0_pReal) then
Ld = Ld + udotn*projection_n(1:3,1:3)
endif
enddo
enddo
damageState(phase)%state(3:11,constituent) = reshape(math_mul33x33(math_inv33(math_I3 - dt*Ld), &
damage_anisoDuctile_getFd0(ipc, ip, el)), &
shape=[9])
end subroutine damage_anisoDuctile_putFd
!--------------------------------------------------------------------------------------------------
!> @brief returns local damage deformation gradient
!--------------------------------------------------------------------------------------------------
pure function damage_anisoDuctile_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_anisoDuctile_getFd0
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase)
damage_anisoDuctile_getFd0 = &
reshape(damageState(phase)%subState0(3:11,constituent),shape=[3,3])
end function damage_anisoDuctile_getFd0
!--------------------------------------------------------------------------------------------------
!> @brief returns local damage deformation gradient
!--------------------------------------------------------------------------------------------------
pure function damage_anisoDuctile_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_anisoDuctile_getPartionedFd0
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_damageInstance(phase)
damage_anisoDuctile_getPartionedFd0 = &
reshape(damageState(phase)%partionedState0(3:11,constituent),shape=[3,3])
end function damage_anisoDuctile_getPartionedFd0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns damage !> @brief returns damage
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

46
code/debug.f90
View File

@ -73,7 +73,8 @@ module debug
debug_MaterialpointLoopDistribution debug_MaterialpointLoopDistribution
integer(pInt), dimension(:,:), allocatable, public :: & integer(pInt), dimension(:,:), allocatable, public :: &
debug_StressLoopDistribution, & !< distribution of stress iterations until convergence debug_StressLoopLiDistribution, & !< distribution of stress iterations until convergence
debug_StressLoopLpDistribution, & !< distribution of stress iterations until convergence
debug_StateLoopDistribution !< distribution of state iterations until convergence debug_StateLoopDistribution !< distribution of state iterations until convergence
real(pReal), public :: & real(pReal), public :: &
@ -137,10 +138,14 @@ subroutine debug_init
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess endif mainProcess
if (allocated(debug_StressLoopDistribution)) & if (allocated(debug_StressLoopLpDistribution)) &
deallocate(debug_StressLoopDistribution) deallocate(debug_StressLoopLpDistribution)
allocate(debug_StressLoopDistribution(nStress+1,2)) allocate(debug_StressLoopLpDistribution(nStress+1,2))
debug_StressLoopDistribution = 0_pInt debug_StressLoopLpDistribution = 0_pInt
if (allocated(debug_StressLoopLiDistribution)) &
deallocate(debug_StressLoopLiDistribution)
allocate(debug_StressLoopLiDistribution(nStress+1,2))
debug_StressLoopLiDistribution = 0_pInt
if (allocated(debug_StateLoopDistribution)) & if (allocated(debug_StateLoopDistribution)) &
deallocate(debug_StateLoopDistribution) deallocate(debug_StateLoopDistribution)
allocate(debug_StateLoopDistribution(nState+1,2)) allocate(debug_StateLoopDistribution(nState+1,2))
@ -311,7 +316,8 @@ subroutine debug_reset
implicit none implicit none
debug_StressLoopDistribution = 0_pInt debug_StressLoopLpDistribution = 0_pInt
debug_StressLoopLiDistribution = 0_pInt
debug_StateLoopDistribution = 0_pInt debug_StateLoopDistribution = 0_pInt
debug_CrystalliteLoopDistribution = 0_pInt debug_CrystalliteLoopDistribution = 0_pInt
debug_MaterialpointStateLoopDistribution = 0_pInt debug_MaterialpointStateLoopDistribution = 0_pInt
@ -378,18 +384,32 @@ subroutine debug_info
endif endif
integral = 0_pInt integral = 0_pInt
write(6,'(3/,a)') 'distribution_StressLoop : stress stiffness' write(6,'(3/,a)') 'distribution_StressLoopLp : stress stiffness'
do j=1_pInt,nStress+1_pInt do j=1_pInt,nStress+1_pInt
if (any(debug_StressLoopDistribution(j,:) /= 0_pInt )) then if (any(debug_StressLoopLpDistribution(j,:) /= 0_pInt )) then
integral = integral + j*(debug_StressLoopDistribution(j,1) + debug_StressLoopDistribution(j,2)) integral = integral + j*(debug_StressLoopLpDistribution(j,1) + debug_StressLoopLpDistribution(j,2))
exceed = ' ' exceed = ' '
if (j > nStress) exceed = '+' ! last entry gets "+" if (j > nStress) exceed = '+' ! last entry gets "+"
write(6,'(i25,a1,i10,1x,i10)') min(nStress,j),exceed,debug_StressLoopDistribution(j,1),& write(6,'(i25,a1,i10,1x,i10)') min(nStress,j),exceed,debug_StressLoopLpDistribution(j,1),&
debug_StressLoopDistribution(j,2) debug_StressLoopLpDistribution(j,2)
endif endif
enddo enddo
write(6,'(a15,i10,2(1x,i10))') ' total',integral,sum(debug_StressLoopDistribution(:,1)), & write(6,'(a15,i10,2(1x,i10))') ' total',integral,sum(debug_StressLoopLpDistribution(:,1)), &
sum(debug_StressLoopDistribution(:,2)) sum(debug_StressLoopLpDistribution(:,2))
integral = 0_pInt
write(6,'(3/,a)') 'distribution_StressLoopLi : stress stiffness'
do j=1_pInt,nStress+1_pInt
if (any(debug_StressLoopLiDistribution(j,:) /= 0_pInt )) then
integral = integral + j*(debug_StressLoopLiDistribution(j,1) + debug_StressLoopLiDistribution(j,2))
exceed = ' '
if (j > nStress) exceed = '+' ! last entry gets "+"
write(6,'(i25,a1,i10,1x,i10)') min(nStress,j),exceed,debug_StressLoopLiDistribution(j,1),&
debug_StressLoopLiDistribution(j,2)
endif
enddo
write(6,'(a15,i10,2(1x,i10))') ' total',integral,sum(debug_StressLoopLiDistribution(:,1)), &
sum(debug_StressLoopLiDistribution(:,2))
integral = 0_pInt integral = 0_pInt
write(6,'(2/,a)') 'distribution_CrystalliteStateLoop :' write(6,'(2/,a)') 'distribution_CrystalliteStateLoop :'

12
code/homogenization.f90
View File

@ -393,8 +393,12 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_F0, & crystallite_F0, &
crystallite_Fp0, & crystallite_Fp0, &
crystallite_Fp, & crystallite_Fp, &
crystallite_Fi0, &
crystallite_Fi, &
crystallite_Lp0, & crystallite_Lp0, &
crystallite_Lp, & crystallite_Lp, &
crystallite_Li0, &
crystallite_Li, &
crystallite_dPdF, & crystallite_dPdF, &
crystallite_dPdF0, & crystallite_dPdF0, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
@ -403,6 +407,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_partionedF, & crystallite_partionedF, &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_partionedLp0, & crystallite_partionedLp0, &
crystallite_partionedFi0, &
crystallite_partionedLi0, &
crystallite_partioneddPdF0, & crystallite_partioneddPdF0, &
crystallite_partionedTstar0_v, & crystallite_partionedTstar0_v, &
crystallite_dt, & crystallite_dt, &
@ -463,6 +469,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) ! ...plastic def grads crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) ! ...plastic def grads
crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads
crystallite_partionedFi0(1:3,1:3,g,i,e) = crystallite_Fi0(1:3,1:3,g,i,e) ! ...intermediate def grads
crystallite_partionedLi0(1:3,1:3,g,i,e) = crystallite_Li0(1:3,1:3,g,i,e) ! ...intermediate velocity grads
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness
crystallite_partionedF0(1:3,1:3,g,i,e) = crystallite_F0(1:3,1:3,g,i,e) ! ...def grads crystallite_partionedF0(1:3,1:3,g,i,e) = crystallite_F0(1:3,1:3,g,i,e) ! ...def grads
crystallite_partionedTstar0_v(1:6,g,i,e) = crystallite_Tstar0_v(1:6,g,i,e) ! ...2nd PK stress crystallite_partionedTstar0_v(1:6,g,i,e) = crystallite_Tstar0_v(1:6,g,i,e) ! ...2nd PK stress
@ -515,6 +523,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_partionedF0(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) ! ...def grads crystallite_partionedF0(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) ! ...def grads
crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fp(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fp(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fi(1:3,1:3,1:myNgrains,i,e) ! ...intermediate def grads
crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) = crystallite_Li(1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e)! ...stiffness crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e)! ...stiffness
crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) = crystallite_Tstar_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) = crystallite_Tstar_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
forall (g = 1:myNgrains) forall (g = 1:myNgrains)
@ -572,6 +582,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
! restore... ! restore...
crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads
crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_Fi(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate def grads
crystallite_Li(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads
crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness
crystallite_Tstar_v(1:6,1:myNgrains,i,e) = crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress crystallite_Tstar_v(1:6,1:myNgrains,i,e) = crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
forall (g = 1:myNgrains) forall (g = 1:myNgrains)

151
code/thermal_adiabatic.f90
View File

@ -40,10 +40,6 @@ module thermal_adiabatic
thermal_adiabatic_stateInit, & thermal_adiabatic_stateInit, &
thermal_adiabatic_aTolState, & thermal_adiabatic_aTolState, &
thermal_adiabatic_LTAndItsTangent, & thermal_adiabatic_LTAndItsTangent, &
thermal_adiabatic_getFT, &
thermal_adiabatic_putFT, &
thermal_adiabatic_getFT0, &
thermal_adiabatic_getPartionedFT0, &
thermal_adiabatic_getTemperature, & thermal_adiabatic_getTemperature, &
thermal_adiabatic_putTemperature, & thermal_adiabatic_putTemperature, &
thermal_adiabatic_getHeatGeneration, & thermal_adiabatic_getHeatGeneration, &
@ -254,7 +250,7 @@ end subroutine thermal_adiabatic_aTolState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_LTAndItsTangent(LT, dLT_dTstar, Tstar_v, Lp, ipc, ip, el) subroutine thermal_adiabatic_LTAndItsTangent(LT, dLT_dTstar3333, Tstar_v, Lp, ipc, ip, el)
use lattice, only: & use lattice, only: &
lattice_massDensity, & lattice_massDensity, &
lattice_specificHeat, & lattice_specificHeat, &
@ -278,9 +274,7 @@ subroutine thermal_adiabatic_LTAndItsTangent(LT, dLT_dTstar, Tstar_v, Lp, ipc, i
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
LT !< thermal velocity gradient LT !< thermal velocity gradient
real(pReal), intent(out), dimension(9,9) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
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) dLT_dTstar3333 !< derivative of LT with respect to Tstar (4th-order tensor)
integer(pInt) :: & integer(pInt) :: &
phase, & phase, &
@ -302,150 +296,9 @@ subroutine thermal_adiabatic_LTAndItsTangent(LT, dLT_dTstar, Tstar_v, Lp, ipc, i
dLT_dTstar3333(i,j,k,l) = Lp(k,l)*lattice_thermalExpansion33(i,j,phase) dLT_dTstar3333(i,j,k,l) = Lp(k,l)*lattice_thermalExpansion33(i,j,phase)
dLT_dTstar3333 = 0.95_pReal*dLT_dTstar3333/(lattice_massDensity(phase)*lattice_specificHeat(phase)) dLT_dTstar3333 = 0.95_pReal*dLT_dTstar3333/(lattice_massDensity(phase)*lattice_specificHeat(phase))
dLT_dTstar = math_Plain3333to99(dLT_dTstar3333)
end subroutine thermal_adiabatic_LTAndItsTangent 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
!--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_putFT(Tstar_v, Lp, dt, ipc, ip, el)
use material, only: &
mappingConstitutive, &
thermalState
use math, only: &
math_Mandel6to33
use lattice, only: &
lattice_massDensity, &
lattice_specificHeat, &
lattice_thermalExpansion33
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(in) :: &
dt
integer(pInt) :: &
phase, &
constituent
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))
thermalState(phase)%state(1,constituent) = thermalState(phase)%subState0(1,constituent) + Tdot*dt
end subroutine thermal_adiabatic_putFT
!--------------------------------------------------------------------------------------------------
!> @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 !> @brief returns temperature based on local damage model state layout
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------