diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 6a9ab9b03..08cfe7902 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -567,12 +567,13 @@ end function constitutive_damagedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el) +subroutine constitutive_microstructure(Tstar_v, Fe, Fp, Lp, subdt, ipc, ip, el) use prec, only: & pReal use material, only: & phase_plasticity, & phase_damage, & + phase_thermal, & phase_vacancy, & material_phase, & PLASTICITY_dislotwin_ID, & @@ -586,7 +587,8 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el) LOCAL_DAMAGE_anisoDuctile_ID, & LOCAL_DAMAGE_gurson_ID, & LOCAL_DAMAGE_phaseField_ID, & - LOCAL_VACANCY_generation_ID + LOCAL_VACANCY_generation_ID, & + LOCAL_THERMAL_adiabatic_ID use plastic_titanmod, only: & plastic_titanmod_microstructure use plastic_nonlocal, only: & @@ -611,6 +613,8 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el) damage_phaseField_microstructure use vacancy_generation, only: & vacancy_generation_microstructure + use thermal_adiabatic, only: & + thermal_adiabatic_microstructure implicit none integer(pInt), intent(in) :: & @@ -621,7 +625,8 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el) Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) real(pReal), intent(in), dimension(3,3) :: & Fe, & !< elastic deformation gradient - Fp !< plastic deformation gradient + Fp, & !< plastic deformation gradient + Lp real(pReal), intent(in) :: & subdt !< timestep @@ -659,6 +664,12 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el) end select + select case (phase_thermal(material_phase(ipc,ip,el))) + case (LOCAL_THERMAL_adiabatic_ID) + call thermal_adiabatic_microstructure(Tstar_v, Lp, subdt, ipc, ip, el) + + end select + select case (phase_vacancy(material_phase(ipc,ip,el))) case (LOCAL_VACANCY_generation_ID) call vacancy_generation_microstructure(Tstar_v, & diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 71f9ec9d3..848c90fca 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -444,6 +444,7 @@ subroutine crystallite_init crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_Lp(1:3,1:3,g,i,e), & crystallite_subdt(g,i,e), g,i,e) ! update dependent state variables to be consistent with basic states enddo enddo @@ -1720,6 +1721,7 @@ subroutine crystallite_integrateStateRK4() call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_Lp(1:3,1:3,g,i,e), & crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2054,6 +2056,7 @@ subroutine crystallite_integrateStateRKCK45() call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_Lp(1:3,1:3,g,i,e), & crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2333,6 +2336,7 @@ subroutine crystallite_integrateStateRKCK45() call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_Lp(1:3,1:3,g,i,e), & crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2583,6 +2587,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_Lp(1:3,1:3,g,i,e), & crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2941,6 +2946,7 @@ eIter = FEsolving_execElem(1:2) call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_Lp(1:3,1:3,g,i,e), & crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -3207,6 +3213,7 @@ subroutine crystallite_integrateStateFPI() call constitutive_microstructure(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & + crystallite_Lp(1:3,1:3,g,i,e), & crystallite_subdt(g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states p = mappingConstitutive(2,g,i,e) c = mappingConstitutive(1,g,i,e) @@ -3669,8 +3676,7 @@ logical function crystallite_integrateStress(& real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK - real(pReal), dimension(9,9) :: dLp_dT_constitutive99, & ! partial derivative of plastic velocity gradient calculated by constitutive law - dRLp_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) + real(pReal), dimension(9,9) :: dRLp_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) dRLp_dLp2, & ! working copy of dRdLp 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 @@ -3910,7 +3916,7 @@ logical function crystallite_integrateStress(& 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(math_Plain3333to99(dFe_dLp3333)) 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(math_Plain3333to99(dLp_dT3333)) 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 >> Lp_constitutive',math_transpose33(Lp_constitutive) diff --git a/code/thermal_adiabatic.f90 b/code/thermal_adiabatic.f90 index 581075f77..cacb90399 100644 --- a/code/thermal_adiabatic.f90 +++ b/code/thermal_adiabatic.f90 @@ -39,6 +39,7 @@ module thermal_adiabatic thermal_adiabatic_init, & thermal_adiabatic_stateInit, & thermal_adiabatic_aTolState, & + thermal_adiabatic_microstructure, & thermal_adiabatic_LTAndItsTangent, & thermal_adiabatic_getTemperature, & thermal_adiabatic_putTemperature, & @@ -247,6 +248,47 @@ subroutine thermal_adiabatic_aTolState(phase,instance) thermalState(phase)%aTolState = tempTol end subroutine thermal_adiabatic_aTolState +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine thermal_adiabatic_microstructure(Tstar_v, Lp, subdt, ipc, ip, el) + use lattice, only: & + lattice_massDensity, & + lattice_specificHeat, & + lattice_thermalExpansion33 + use material, only: & + mappingConstitutive, & + phase_thermalInstance, & + thermalState + use math, only: & + math_Mandel6to33 + + 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) :: & + subdt + integer(pInt) :: & + phase, & + constituent + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + + thermalState(phase)%state(1,constituent) = & + thermalState(phase)%subState0(1,constituent) + & + subdt* & + 0.95_pReal*sum(abs(math_Mandel6to33(Tstar_v))*Lp)/ & + (lattice_massDensity(phase)*lattice_specificHeat(phase)) + +end subroutine thermal_adiabatic_microstructure + !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !--------------------------------------------------------------------------------------------------