diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 6cf7b5e46..149d09a49 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -200,11 +200,11 @@ module constitutive integer, intent(in) :: co, ip, el end subroutine constitutive_mech_setF - module subroutine constitutive_thermal_setT(T,co,ce) - real(pReal), intent(in) :: T - integer, intent(in) :: co, ce - end subroutine constitutive_thermal_setT - + module subroutine constitutive_thermal_setField(T,dot_T, co,ce) + real(pReal), intent(in) :: T, dot_T + integer, intent(in) :: ce, co + end subroutine constitutive_thermal_setField + module subroutine constitutive_damage_set_phi(phi,co,ce) real(pReal), intent(in) :: phi integer, intent(in) :: co, ce @@ -353,7 +353,7 @@ module constitutive constitutive_restartWrite, & constitutive_restartRead, & integrateDamageState, & - constitutive_thermal_setT, & + constitutive_thermal_setField, & constitutive_damage_set_phi, & constitutive_damage_get_phi, & constitutive_mech_getP, & diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index a82262e81..57b8a3117 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -272,15 +272,16 @@ end function thermal_T !---------------------------------------------------------------------------------------------- !< @brief Set temperature !---------------------------------------------------------------------------------------------- -module subroutine constitutive_thermal_setT(T,co,ce) +module subroutine constitutive_thermal_setField(T,dot_T, co,ce) - real(pReal), intent(in) :: T + real(pReal), intent(in) :: T, dot_T integer, intent(in) :: ce, co current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T + current(material_phaseAt2(co,ce))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T -end subroutine constitutive_thermal_setT +end subroutine constitutive_thermal_setField diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index d164f11e7..44d9e22b1 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -191,6 +191,9 @@ function grid_thermal_spectral_solution(timeinc) result(solution) call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), & (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & 1,ce) + call homogenization_thermal_setField(T_current(i,j,k), & + (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & + ce) homogenization_T(ce) = T_current(i,j,k) enddo; enddo; enddo @@ -233,6 +236,9 @@ subroutine grid_thermal_spectral_forward(cutBack) (T_current(i,j,k) - & T_lastInc(i,j,k))/params%timeinc, & 1,ce) + call homogenization_thermal_setField(T_current(i,j,k), & + (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & + ce) homogenization_T(ce) = T_current(i,j,k) enddo; enddo; enddo else @@ -283,8 +289,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = ce + 1 call thermal_conduction_getSource(Tdot, 1,ce) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & - + thermal_conduction_getMassDensity (1,ce)* & - thermal_conduction_getSpecificHeat(1,ce)*(T_lastInc(i,j,k) - & + + thermal_conduction_getMassDensity (ce)* & + thermal_conduction_getSpecificHeat(ce)*(T_lastInc(i,j,k) - & T_current(i,j,k))& + mu_ref*T_current(i,j,k) enddo; enddo; enddo @@ -315,7 +321,7 @@ subroutine updateReference do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 K_ref = K_ref + thermal_conduction_getConductivity(1,ce) - mu_ref = mu_ref + thermal_conduction_getMassDensity(1,ce)* thermal_conduction_getSpecificHeat(1,ce) + mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce) enddo; enddo; enddo K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 86caafef3..f316d3969 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -75,8 +75,7 @@ module homogenization el !< element number end subroutine mech_partition - module subroutine thermal_partition(T,ce) - real(pReal), intent(in) :: T + module subroutine thermal_partition(ce) integer, intent(in) :: ce end subroutine thermal_partition @@ -112,11 +111,40 @@ module homogenization logical, dimension(2) :: doneAndHappy end function mech_updateState + + module function thermal_conduction_getConductivity(ip,el) result(K) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: K + + end function thermal_conduction_getConductivity + + module function thermal_conduction_getSpecificHeat(ce) result(c_P) + integer, intent(in) :: ce + real(pReal) :: c_P + end function thermal_conduction_getSpecificHeat + + module function thermal_conduction_getMassDensity(ce) result(rho) + integer, intent(in) :: ce + real(pReal) :: rho + end function thermal_conduction_getMassDensity + + module subroutine homogenization_thermal_setField(T,dot_T, ce) + integer, intent(in) :: ce + real(pReal), intent(in) :: T, dot_T + end subroutine homogenization_thermal_setField + end interface public :: & homogenization_init, & materialpoint_stressAndItsTangent, & + thermal_conduction_getSpecificHeat, & + thermal_conduction_getConductivity, & + thermal_conduction_getMassDensity, & + homogenization_thermal_setfield, & homogenization_forward, & homogenization_results, & homogenization_restartRead, & @@ -276,7 +304,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE ho = material_homogenizationAt(el) do ip = FEsolving_execIP(1),FEsolving_execIP(2) ce = (el-1)*discretization_nIPs + ip - call thermal_partition(homogenization_T(ce),ce) + call thermal_partition(ce) do co = 1, homogenization_Nconstituents(ho) ph = material_phaseAt(co,el) if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index ce55d1573..80137f13a 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -3,6 +3,7 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_thermal + use lattice type :: tDataContainer real(pReal), dimension(:), allocatable :: T, dot_T @@ -67,15 +68,18 @@ end subroutine thermal_init !-------------------------------------------------------------------------------------------------- !> @brief Partition temperature onto the individual constituents. !-------------------------------------------------------------------------------------------------- -module subroutine thermal_partition(T,ce) +module subroutine thermal_partition(ce) - real(pReal), intent(in) :: T integer, intent(in) :: ce + real(pReal) :: T, dot_T integer :: co + + T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) + dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - call constitutive_thermal_setT(T,co,ce) + call constitutive_thermal_setField(T,dot_T,co,ce) enddo end subroutine thermal_partition @@ -93,4 +97,89 @@ module subroutine thermal_homogenize(ip,el) end subroutine thermal_homogenize +!-------------------------------------------------------------------------------------------------- +!> @brief return homogenized thermal conductivity in reference configuration +!-------------------------------------------------------------------------------------------------- +module function thermal_conduction_getConductivity(ip,el) result(K) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: K + + integer :: & + co + + + K = 0.0_pReal + + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + K = K + crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el))) + enddo + + K = K / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + +end function thermal_conduction_getConductivity + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized specific heat capacity +!-------------------------------------------------------------------------------------------------- +module function thermal_conduction_getSpecificHeat(ce) result(c_P) + + integer, intent(in) :: ce + real(pReal) :: c_P + + integer :: co + + + c_P = 0.0_pReal + + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + c_P = c_P + lattice_c_p(material_phaseAt2(co,ce)) + enddo + + c_P = c_P / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + +end function thermal_conduction_getSpecificHeat + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized mass density +!-------------------------------------------------------------------------------------------------- +module function thermal_conduction_getMassDensity(ce) result(rho) + + integer, intent(in) :: ce + real(pReal) :: rho + + integer :: co + + + rho = 0.0_pReal + + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + rho = rho + lattice_rho(material_phaseAt2(co,ce)) + enddo + + rho = rho / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + +end function thermal_conduction_getMassDensity + + + +!-------------------------------------------------------------------------------------------------- +!> @brief Set thermal field and its rate (T and dot_T) +!-------------------------------------------------------------------------------------------------- +module subroutine homogenization_thermal_setField(T,dot_T, ce) + + integer, intent(in) :: ce + real(pReal), intent(in) :: T, dot_T + + + current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) = T + current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) = dot_T + + +end subroutine homogenization_thermal_setField + end submodule homogenization_thermal diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 1711c5512..ccca75026 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -26,9 +26,6 @@ module thermal_conduction public :: & thermal_conduction_init, & thermal_conduction_getSource, & - thermal_conduction_getConductivity, & - thermal_conduction_getSpecificHeat, & - thermal_conduction_getMassDensity, & thermal_conduction_putTemperatureAndItsRate, & thermal_conduction_results @@ -110,90 +107,6 @@ subroutine thermal_conduction_getSource(Tdot, ip,el) end subroutine thermal_conduction_getSource -!-------------------------------------------------------------------------------------------------- -!> @brief return homogenized thermal conductivity in reference configuration -!-------------------------------------------------------------------------------------------------- -function thermal_conduction_getConductivity(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - thermal_conduction_getConductivity - - integer :: & - co - - - thermal_conduction_getConductivity = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - thermal_conduction_getConductivity = thermal_conduction_getConductivity + & - crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el))) - enddo - - thermal_conduction_getConductivity = thermal_conduction_getConductivity & - / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function thermal_conduction_getConductivity - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized specific heat capacity -!-------------------------------------------------------------------------------------------------- -function thermal_conduction_getSpecificHeat(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal) :: & - thermal_conduction_getSpecificHeat - - integer :: & - co - - - thermal_conduction_getSpecificHeat = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & - + lattice_c_p(material_phaseAt(co,el)) - enddo - - thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & - / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function thermal_conduction_getSpecificHeat - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized mass density -!-------------------------------------------------------------------------------------------------- -function thermal_conduction_getMassDensity(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal) :: & - thermal_conduction_getMassDensity - - integer :: & - co - - - thermal_conduction_getMassDensity = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & - + lattice_rho(material_phaseAt(co,el)) - enddo - - thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & - / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function thermal_conduction_getMassDensity - - !-------------------------------------------------------------------------------------------------- !> @brief updates thermal state with solution from heat conduction PDE !--------------------------------------------------------------------------------------------------