diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 621a091d8..6480f7382 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -32,10 +32,10 @@ #include "phase_mechanics_eigendeformation.f90" #include "phase_mechanics_eigendeformation_cleavageopening.f90" #include "phase_mechanics_eigendeformation_slipplaneopening.f90" +#include "phase_mechanics_eigendeformation_thermalexpansion.f90" #include "phase_thermal.f90" #include "phase_thermal_dissipation.f90" #include "phase_thermal_externalheat.f90" -#include "phase_mechanics_eigendeformation_thermalexpansion.f90" #include "phase_damage.f90" #include "phase_damage_isobrittle.f90" #include "phase_damage_isoductile.f90" diff --git a/src/phase.f90 b/src/phase.f90 index f37efd043..d532e4596 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -186,6 +186,11 @@ module phase real(pReal) :: T end function thermal_T + module function thermal_dot_T(ph,me) result(dot_T) + integer, intent(in) :: ph,me + real(pReal) :: dot_T + end function thermal_dot_T + module subroutine constitutive_mech_setF(F,co,ip,el) real(pReal), dimension(3,3), intent(in) :: F diff --git a/src/phase_mechanics.f90 b/src/phase_mechanics.f90 index 91d455279..2b2144b97 100644 --- a/src/phase_mechanics.f90 +++ b/src/phase_mechanics.f90 @@ -137,6 +137,11 @@ submodule(phase) mechanics logical, dimension(:,:), allocatable :: myKinematics end function kinematics_slipplane_opening_init + module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics) + integer, intent(in) :: kinematics_length + logical, dimension(:,:), allocatable :: myKinematics + end function kinematics_thermal_expansion_init + module subroutine plastic_isotropic_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group @@ -365,6 +370,7 @@ module subroutine mech_init(phases) if(maxval(phase_Nkinematics) /= 0) then where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_cleavage_opening_ID where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID + where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_thermal_expansion_ID endif end subroutine mech_init diff --git a/src/phase_mechanics_eigendeformation_thermalexpansion.f90 b/src/phase_mechanics_eigendeformation_thermalexpansion.f90 index b9483977a..7a6a7daa3 100644 --- a/src/phase_mechanics_eigendeformation_thermalexpansion.f90 +++ b/src/phase_mechanics_eigendeformation_thermalexpansion.f90 @@ -3,7 +3,7 @@ !> @brief material subroutine incorporating kinematics resulting from thermal expansion !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:thermal) thermalexpansion +submodule(phase:eigendeformation) thermalexpansion integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance @@ -93,21 +93,21 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) real(pReal) :: T, dot_T - T = current(ph)%T(me) - dot_T = current(ph)%dot_T(me) + T = thermal_T(ph,me) + dot_T = thermal_dot_T(ph,me) associate(prm => param(kinematics_thermal_expansion_instance(ph))) - Li = dot_T * ( & - prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient - + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient - + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient - ) / & - (1.0_pReal & - + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. & - + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. & - + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. & - ) - end associate + Li = dot_T * ( & + prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient + + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient + + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient + ) / & + (1.0_pReal & + + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. & + + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. & + + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. & + ) + end associate dLi_dTstar = 0.0_pReal end subroutine thermalexpansion_LiAndItsTangent diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 96ac291ae..621bb4f4a 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -32,10 +32,7 @@ submodule(phase) thermal logical, dimension(:,:), allocatable :: mySources end function externalheat_init - module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics - end function kinematics_thermal_expansion_init + module subroutine externalheat_dotState(ph, me) @@ -121,11 +118,6 @@ module subroutine thermal_init(phases) maxval(thermalState(ph)%p%sizeDotState)) enddo PhaseLoop2 -!-------------------------------------------------------------------------------------------------- -!initialize kinematic mechanisms - if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) & - phase_kinematics = KINEMATICS_thermal_expansion_ID - end subroutine thermal_init @@ -257,6 +249,20 @@ module function thermal_T(ph,me) result(T) end function thermal_T +!---------------------------------------------------------------------------------------------- +!< @brief Get rate of temperature (for use by non-thermal physics) +!---------------------------------------------------------------------------------------------- +module function thermal_dot_T(ph,me) result(dot_T) + + integer, intent(in) :: ph, me + real(pReal) :: dot_T + + + dot_T = current(ph)%dot_T(me) + +end function thermal_dot_T + + !---------------------------------------------------------------------------------------------- !< @brief Set temperature !----------------------------------------------------------------------------------------------