using new structure for thermal

This commit is contained in:
Martin Diehl 2021-01-24 13:26:01 +01:00
parent a933fe57e7
commit e22d76be9e
6 changed files with 142 additions and 105 deletions

View File

@ -200,10 +200,10 @@ module constitutive
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
end subroutine constitutive_mech_setF end subroutine constitutive_mech_setF
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) :: co, ce integer, intent(in) :: ce, co
end subroutine constitutive_thermal_setT end subroutine constitutive_thermal_setField
module subroutine constitutive_damage_set_phi(phi,co,ce) module subroutine constitutive_damage_set_phi(phi,co,ce)
real(pReal), intent(in) :: phi real(pReal), intent(in) :: phi
@ -353,7 +353,7 @@ module constitutive
constitutive_restartWrite, & constitutive_restartWrite, &
constitutive_restartRead, & constitutive_restartRead, &
integrateDamageState, & integrateDamageState, &
constitutive_thermal_setT, & constitutive_thermal_setField, &
constitutive_damage_set_phi, & constitutive_damage_set_phi, &
constitutive_damage_get_phi, & constitutive_damage_get_phi, &
constitutive_mech_getP, & constitutive_mech_getP, &

View File

@ -272,15 +272,16 @@ end function thermal_T
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Set temperature !< @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 integer, intent(in) :: ce, co
current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T 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

View File

@ -191,6 +191,9 @@ function grid_thermal_spectral_solution(timeinc) result(solution)
call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), & call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), &
(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, &
1,ce) 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) homogenization_T(ce) = T_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
@ -233,6 +236,9 @@ subroutine grid_thermal_spectral_forward(cutBack)
(T_current(i,j,k) - & (T_current(i,j,k) - &
T_lastInc(i,j,k))/params%timeinc, & T_lastInc(i,j,k))/params%timeinc, &
1,ce) 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) homogenization_T(ce) = T_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
else else
@ -283,8 +289,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = ce + 1 ce = ce + 1
call thermal_conduction_getSource(Tdot, 1,ce) call thermal_conduction_getSource(Tdot, 1,ce)
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) &
+ thermal_conduction_getMassDensity (1,ce)* & + thermal_conduction_getMassDensity (ce)* &
thermal_conduction_getSpecificHeat(1,ce)*(T_lastInc(i,j,k) - & thermal_conduction_getSpecificHeat(ce)*(T_lastInc(i,j,k) - &
T_current(i,j,k))& T_current(i,j,k))&
+ mu_ref*T_current(i,j,k) + mu_ref*T_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
@ -315,7 +321,7 @@ subroutine updateReference
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
K_ref = K_ref + thermal_conduction_getConductivity(1,ce) 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 enddo; enddo; enddo
K_ref = K_ref*wgt K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)

View File

@ -75,8 +75,7 @@ module homogenization
el !< element number el !< element number
end subroutine mech_partition end subroutine mech_partition
module subroutine thermal_partition(T,ce) module subroutine thermal_partition(ce)
real(pReal), intent(in) :: T
integer, intent(in) :: ce integer, intent(in) :: ce
end subroutine thermal_partition end subroutine thermal_partition
@ -112,11 +111,40 @@ module homogenization
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
end function mech_updateState 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 end interface
public :: & public :: &
homogenization_init, & homogenization_init, &
materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent, &
thermal_conduction_getSpecificHeat, &
thermal_conduction_getConductivity, &
thermal_conduction_getMassDensity, &
homogenization_thermal_setfield, &
homogenization_forward, & homogenization_forward, &
homogenization_results, & homogenization_results, &
homogenization_restartRead, & homogenization_restartRead, &
@ -276,7 +304,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
do ip = FEsolving_execIP(1),FEsolving_execIP(2) do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip ce = (el-1)*discretization_nIPs + ip
call thermal_partition(homogenization_T(ce),ce) call thermal_partition(ce)
do co = 1, homogenization_Nconstituents(ho) do co = 1, homogenization_Nconstituents(ho)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then

View File

@ -3,6 +3,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_thermal submodule(homogenization) homogenization_thermal
use lattice
type :: tDataContainer type :: tDataContainer
real(pReal), dimension(:), allocatable :: T, dot_T real(pReal), dimension(:), allocatable :: T, dot_T
@ -67,15 +68,18 @@ end subroutine thermal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Partition temperature onto the individual constituents. !> @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 integer, intent(in) :: ce
real(pReal) :: T, dot_T
integer :: co 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)) 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 enddo
end subroutine thermal_partition end subroutine thermal_partition
@ -93,4 +97,89 @@ module subroutine thermal_homogenize(ip,el)
end subroutine thermal_homogenize 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 end submodule homogenization_thermal

View File

@ -26,9 +26,6 @@ module thermal_conduction
public :: & public :: &
thermal_conduction_init, & thermal_conduction_init, &
thermal_conduction_getSource, & thermal_conduction_getSource, &
thermal_conduction_getConductivity, &
thermal_conduction_getSpecificHeat, &
thermal_conduction_getMassDensity, &
thermal_conduction_putTemperatureAndItsRate, & thermal_conduction_putTemperatureAndItsRate, &
thermal_conduction_results thermal_conduction_results
@ -110,90 +107,6 @@ subroutine thermal_conduction_getSource(Tdot, ip,el)
end subroutine thermal_conduction_getSource 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 !> @brief updates thermal state with solution from heat conduction PDE
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------