2.2.4. Heat Transfer Solver (HEATTRANS module)

2.2.4.1. Governing equations

Tdyn solves the transient Heat Transfer Equations in a given domain \(\Omega := \left\{ \Omega_F \cup \Omega_S \right\}\) (being \(\Omega_F\) the fluid domain and \(\Omega_S\) the solid domain) and time interval \((0, t)\):

(2.335)\[\begin{aligned} & \rho C_P \left( \frac{\partial \theta}{\partial t} + \left( \boldsymbol{u} \cdot \nabla \right) \theta \right) - \nabla \cdot \left( \kappa _F \nabla \theta \right) = \rho _F q_F & \text{in } \Omega_F \times (0,t) \end{aligned}\]

(2.336)\[\begin{aligned} & \rho _S C \frac{\partial \theta}{\partial t} - \nabla\cdot (\kappa _S \nabla\theta) = \rho _S q_S & \text{in } \Omega_S \times (0,t) \end{aligned}\]

where \(\theta = \theta (\boldsymbol{x}, t)\) denotes the temperature field, \(\rho _S\) the (constant) solid density, \(k_S\) the solid thermal conduction tensor, \(k_F\) the fluid thermal conduction constant, \(C_P\) the fluid specific heat constant, \(C\) the solid specific heat constant and \(q_S\) and \(q_F\) the solid and fluid, volumetric heat source distributions respectively. The above equations need to be combined with the following boundary conditions:

(2.337)\[\begin{aligned} & \theta = \theta_c & \text{in } \Gamma _\theta \times (0,t) \end{aligned}\]

(2.338)\[\begin{aligned} & n\kappa _F \nabla\theta = f_F & \text{in } \Gamma _F \times (0,T) \end{aligned}\]

(2.339)\[\begin{aligned} & n\kappa _S \nabla\theta = f_S & \text{in } \Gamma _S \times (0,T) \end{aligned}\]

(2.340)\[\begin{aligned} & \theta(x,0) = \theta_0 (x) & \text{in } \Omega \times \left\{ 0 \right\} \end{aligned}\]

In the above equations, \(\Gamma := \partial \Omega\) denotes the boundary of the domain \(\Omega\), with \(\boldsymbol{n}\) the normal unit vector, \(\theta_c\) is the temperature field on \(\Gamma_{\theta}\) (the part of the boundary of Dirichlet type, or prescribed temperature type), \(f_F\) the prescribed heat flux on \(\Gamma_F\) (prescribed heat flux fluid boundary), \(f_S\) the prescribed heat flux on \(\Gamma_S\) (prescribed heat flux solid boundary) and \(\theta_0\) the initial temperature field.

The spatial discretization of the heat transfer equations is done by means of the finite element method, while for the time discretization implicit first and second order schemes have been implemented.

2.2.4.2. Stabilized heat transfer equations

The Finite Increment Calculus theory presented above is also valid for the heat transfer balance equations Eq. (4-1). The stabilized form of the governing differential equations in 3 dimensions is written as follows:

(2.341)\[\begin{aligned} & r_F - \frac{h_{Fj}}{2} \frac{\partial r_F}{\partial x_j} = 0 & \text{in } \Omega_F, j = 1,2,3 \end{aligned}\]

(2.342)\[\begin{aligned} & r_S - \frac{h_{Sj}}{2} \frac{\partial r_S}{\partial x_j} = 0 & \text{in }\Omega_S, j = 1,2,3 \end{aligned}\]

Where

(2.343)\[\begin{aligned} & r_F = \rho C_P \frac{\partial \theta}{\partial t} - \kappa _F \Delta \theta - \rho q_F \end{aligned}\]

(2.344)\[\begin{aligned} & r_S = \rho _S C \frac{\partial \theta}{\partial t} - \nabla \cdot \left( \kappa _S \nabla\theta \right) - \rho _S q_S \end{aligned}\]

and the new boundary conditions are:

(2.345)\[\begin{aligned} & \theta = \theta_c & \text{in } \Gamma _\theta \times (0,t) \end{aligned}\]

(2.346)\[\begin{aligned} & n\kappa _F \nabla\theta + n \frac{h_{Fj}}{2} r_F = f_F & \text{in } \Gamma _F \times (0,T) \end{aligned}\]

(2.347)\[\begin{aligned} & n\kappa _S \nabla\theta + n \frac{h_{Sj}}{2} r_S = f_S & \text{in } \Gamma _S \times (0,T) \end{aligned}\]

(2.348)\[\begin{aligned} & \theta(x,0) = \theta_0 (x) & \text{in } \Omega \times \left\{ 0 \right\} \end{aligned}\]

The stabilized heat transfer balance equations are discretized in time in the standard way, to give the following equations:

(2.349)\[\begin{aligned} & \theta^{j+1,n+1}= \theta^n - \frac{\Delta t}{\rho C_P} \left[ \rho C_P \left( \boldsymbol{u} \cdot \nabla \right) \theta - \nabla \cdot \left( \kappa _F \nabla \theta \right) - \rho q_F \right]^{j,n+\theta} \end{aligned}\]

(2.350)\[\begin{aligned} & \theta^{j+1,n+1}= \theta^n - \frac{\Delta t}{\rho C} \left[ \nabla \cdot \left( \kappa _S \nabla \theta \right) - \rho q_S \right]^{j,n+\theta} \end{aligned}\]

2.2.4.3. Radiation models

Heat transfer by radiation can be taken into account within Tdyn by using one of the following radiation models.

2.2.4.3.1. P-1 radiation model

The energy conservation equation in terms of the temperature \(T\), can be written as:

(2.351)\[\begin{aligned} & \rho c \left[ \frac{\partial T}{\partial t} +u_j \frac{\partial T}{\partial x_j} \right] - \frac{\partial}{\partial x_j} \left( k \cdot \frac{\partial T}{\partial x_j} \right) + r T + p \frac{\partial u_j}{\partial x_j} - Φ = s_T \end{aligned}\]

Being \(c\) the heat capacity, \(k\) the thermal conductivity, \(r\) the heat reaction term, \(s\) the heat source term, and \(Φ\) the viscous dissipation function given by:

(2.352)\[\begin{aligned} & Φ = \boldsymbol{\tau} : \nabla u \end{aligned}\]

Where \(\tau\) is the viscous stress tensor of the fluid. The FEM formulation for the temperature equation can be written as:

(2.353)\[\begin{split}\begin{aligned} & \int_\Omega W \rho c \frac{\boldsymbol{T}^{n+\theta,k} - T^n}{\theta ∆t} d\Omega + \int_\Omega W \rho c \left[ u_j^{n+\theta,k-1} \frac{\partial}{\partial x_j} T^{n+\theta,k} \right] d\Omega + \\ & + \int_\Omega \tau \rho c \left[ u_j^{n+\theta,k-1} \frac{\partial W}{\partial x_j} \right] \left[ u_j^{n+\theta,k-1} \frac{\partial}{\partial x_j} T^{n+\theta,k} \right] d\Omega + \\ & + \int_\Omega k \frac{\partial W}{\partial x_j} \frac{\partial}{\partial x_j} T^{n+\theta,k} d\Omega + \int_\Omega W r T^{n+\theta,k} d\Omega = \\ & = -\int_\Omega W s_{T+Φ+p} d\Omega + \int_{\Gamma}n_j W q_j^* d\Gamma + \\ & + \int_\Omega \tau \rho c \left[ u_j^{n+\theta,k-1} \frac{\partial W}{\partial x_j} \right] ς_i^{n+\theta,k-1} d\Omega \end{aligned}\end{split}\]

(2.354)\[\begin{aligned} & \int_\Omega W ς_i^{n+1,k} d\Omega = \int_\Omega W \left[ u_j^{n+\theta,k} \frac{\partial}{\partial x_j} T^{n+\theta,k} \right] d\Omega \end{aligned}\]

Where \(s_{T + Φ + p}\) is the total heat/sink source term, including viscous dissipation and rate of work for volume change, and \(q^*\) is the heat flux at the walls.

The so called P-N approximation reduces the integral equations of radiative transfer in media to differential equations by approximating the transfer relations with a finite set of moment equations. The starting point in the following equation for the variation of intensity of radiation \(i^{'} (W/m2 \cdot sr)\) at a position \(\boldsymbol{x}\), along the \(\boldsymbol{s}\) direction [Siegel_1992]:

(2.355)\[\begin{aligned} & \frac{di_{\lambda}^{'}}{ds} = a_{\lambda} i_{\lambda b}^{'} - \left( a_{\lambda} + \sigma_{s \lambda} \right) i_\lambda^{'} + \frac{\sigma_{s \lambda}}{4 \pi} \int_0^{4 \pi} i_{\lambda}^{'} \cdot ϕ \cdot d\omega \end{aligned}\]

Where \(a_{\lambda}\) is the absorption coefficient and \(\sigma_{s \lambda}\) is the dispersion coefficient for the given wave-length \(\lambda\), and the sub-index \(b\) refers to the black body. If the medium is assumed gray, with uniform scattering and absorption coefficients, and to scatter isotropically, above equation can be integrated over all \(\lambda\), to give:

(2.356)\[\begin{aligned} & \frac{di^{'}}{ds} = a i_b^{'} - \left( a + \sigma_s \right) i^{'} + \frac{\sigma_s}{4 \pi} \int_0^{4 \pi} i^{'} \cdot d\omega \end{aligned}\]

Where the first term at the right hand side represents the gain of intensity of radiation by emission, the second term the losses by absorption and scattering, and the third term the gain by scattering into \(\boldsymbol{s}\) direction. To develop the P-N method, the intensity \(i^{'}\) is expanded in an orthogonal series of spherical harmonics. The \(P-1\) approximation is obtained by retaining the first two terms of the series, and gives [Siegel_1992]:

(2.357)\[\begin{aligned} & \sum_{i=1}^3 \frac{\partial i^{'(i)}}{\partial x_i} = a \left( 4 \pi i_b^{'} - i^{'(0)} \right) \end{aligned}\]

In the above equation \(i^{'(0)}\) and \(i^{'(i)}\) represents the zeroth and first moments of the radiation intensity:

(2.358)\[\begin{aligned} & i^{'(0)} = G = \int_0^{4 \pi} i^{'} \cdot d\omega \end{aligned}\]

(2.359)\[\begin{split}\begin{aligned} & i^{'(i)} = q_{r,i} = \int_0^{4 \pi} i^{'} \cdot s_i \cdot d\omega \\ & \boldsymbol{q}_r = \int_0^{4 \pi} i^{'} \cdot \boldsymbol{s} \cdot d\omega \end{aligned}\end{split}\]

Where \(\boldsymbol{q}_r\) is the radiative heat flux vector, and \(G\) is the incident radiation. Therefore, the \(P-1\) approximation yields:

(2.360)\[\begin{aligned} & -\nabla \cdot \boldsymbol{q}_r = a \left( G - 4 \pi i_b^{'} \right) \end{aligned}\]

The first term in the above equation, represents the net radiative energy supplied per unit volume. Furthermore, \(P-1\) model gives a following additional relation between \(G\) and \(q_r\):

(2.361)\[\begin{aligned} & \boldsymbol{q}_r = - \frac{1}{3 \left( a + \sigma_s \right)} \nabla G = - \gamma \nabla G \end{aligned}\]

Combining above relations, the resulting transport equation for \(G\) is:

(2.362)\[\begin{aligned} & \nabla \cdot \left( \gamma \nabla G \right) - a G = - 4 e \sigma T^4 \end{aligned}\]

Where \(\sigma = 5.672 \cdot 10-8 W/m2\cdot K4\) is the Stefan-Boltzmann constant, \(\epsilon\) is the emissivity of the material, and the black body radiant energy has been evaluated using:

(2.363)\[\begin{aligned} & i_b^{'} = \frac{ \sigma T^4 }{\pi} \end{aligned}\]

The FEM formulation for the incident radiation equation can be written as:

(2.364)\[\begin{aligned} & \int_\Omega \gamma \frac{\partial W}{\partial x_j} \frac{\partial}{\partial x_j} \boldsymbol{G}^k d\Omega + \int_\Omega W a \boldsymbol{G}^k d\Omega = \int_\Omega W 4 \epsilon \sigma \left( T^{n+\theta,k} \right)^4 d\Omega + \int_{\Gamma} W q_{r,w} d\Gamma \end{aligned}\]

Where \(q_(r,w)\) is the flux of the incident radiation al the walls. From the solution of the previous equation, the heat source term that must be added to the energy equation (temperature) is calculated as follows:

(2.365)\[\begin{aligned} & s_T = - \nabla \cdot \boldsymbol{q}_r = a G - 4 \epsilon \sigma T^4 \end{aligned}\]

In the implemented model the emissivity is assumed to be equal to the absorptivity \(a\), and therefore the above equations can be written as:

(2.366)\[\begin{aligned} & \int_\Omega \gamma \frac{\partial W}{\partial x_j} \frac{\partial}{\partial x_j} \boldsymbol{G}^k d\Omega + \int_\Omega W a \boldsymbol{G}^k d\Omega = \int_\Omega W 4 a \sigma \left( T^{n+\theta,k} \right)^4 d\Omega+ \int_{\Gamma} W g_w d\Gamma \end{aligned}\]

(2.367)\[\begin{aligned} & s_T = - \nabla \cdot \boldsymbol{q}_r = a \left( G - 4 \sigma T^4 \right) \end{aligned}\]

The boundary condition for the incident radiation equation at the walls is given by (assuming an inwards normal vector \(\boldsymbol{n}\)):

(2.368)\[\begin{aligned} & g_w = \gamma \nabla G \cdot \boldsymbol{n} = \gamma \frac{\partial G}{\partial n} \end{aligned}\]

Where \(g_w\) can be calculated, assuming that walls are diffuse gray surfaces, as [21]:

(2.369)\[\begin{aligned} & g_w = \frac{\epsilon_w}{2 ( 2 -\epsilon_w )} \left(4 \sigma T_w^4 -G_w \right) \end{aligned}\]

Where \(\epsilon_w\) is the wall emissivity, and \(T_w\) is the wall temperature.

Remark: It is usually assumed that the emissivity at the inlet and outlets is 1 (black body).

The final FEM formulation for the incident radiation equation can be written as:

(2.370)\[\begin{split}\begin{aligned} & \int_\Omega \Gamma \frac{\partial W}{\partial x_j} \frac{\partial}{\partial x_j} \boldsymbol{G}^k d\Omega + \int_\Omega W a \boldsymbol{G}^k d\Omega + \frac{\epsilon _w}{2(2-\epsilon _w )} \int_{\Gamma} W \boldsymbol{G}^k d\Gamma = \\ & = \int_\Omega W 4 a \sigma \left( T^{n+\theta,k} \right)^4 d\Omega + \frac{\epsilon _w}{2(2-\epsilon _w )} \int_{\Gamma} W 4 \sigma T_w^4 d\Gamma \end{aligned}\end{split}\]

2.2.4.3.2. Surface to surface (S2S) radiation model

The surface-to-surface \((S2S)\) radiation model can be used to account for the radiation exchange in an enclosure of gray-diffuse surfaces. This energy exchange depends on:

• Surfaces size

• Separation distance

• Surfaces orientation

All these factors are casted together into a geometric function called the view factor \(F_{ij}\) that can be computed as follows:

(2.371)\[\begin{aligned} & F_{ij} = \frac{1}{A_i} \int_{A_i} \int_{A_j} \frac{cos \theta_i \cdot cos \theta_j}{ \pi r^2} \delta_{ij} dA_i dA_j \end{aligned}\]

where \(\delta_{ii} = 0\), \(\delta_{ij} = 0\) when \(i, j\) are blocking surfaces and \(\delta_{ij} = 1\) when \(i,j\) are non-blocking surfaces.

The S2S model main assumptions are:

• All surfaces are gray and diffuse. According to this, if a certain amount of radiation is incident on a surface, then a fraction is reflected, a fraction is absorbed, and a fraction is transmitted.

• Absorption, emission and scattering of radiation by the medium can be ignored (non-participating media)

• Emissivity and absorptivity are independent of the radiation wave length

• Emissivity is taken to be identical to absorptivity (\(\epsilon ≡ \alpha\))

• Reflectivity is independent of the outgoing or the incoming directions

In most applications of the \(S2S\) model the surfaces are opaque to thermal radiation (in the infrared spectrum). In such a situation the following relations are fulfilled:

• \(\tau = 0\)

• \(\alpha = \epsilon\)

• \(\rho = 1 - \epsilon\)

and the only independent parameter characteristic of the material is the emissivity \(\epsilon\). When the \(S2S\) model is used, it is also possible to define a “partial enclosure” that allows disabling the view factor calculation for walls with negligible emission/absorption or walls that have uniform temperature.

The energy flux leaving a given surface is composed of directly emitted and reflected energy. Hence, the total energy flux leaving the surface \(k (q_{out,k})\) has two contributions, one depending on its own temperature and emissivity, and another one depending on the energy flux incident on the surface from the surroundings.

(2.372)\[\begin{aligned} & q_{out,k} = \epsilon_k \cdot \sigma \cdot T_k^4 + \rho_k \cdot q_{in,k} \end{aligned}\]

The energy flux incident from the surroundings is given by:

(2.373)\[\begin{aligned} & q_{in,k} = \sum_{j=1}^N F_{kj} \cdot q_{out,j} \end{aligned}\]

where \(N\) is the number of surface clusters (i.e. element patches) considered. Hence, the total energy flux leaving a given surface can be expressed as follows:

(2.374)\[\begin{aligned} & q_{out,k} = \epsilon_k \cdot \sigma \cdot T_k^4 + \rho_k \sum_{j=1}^N F_{kj} \cdot q_{out,j} \end{aligned}\]

(2.375)\[\begin{aligned} & q_{out,k} -\rho_k \sum_{j=1}^N F_{kj} \cdot q_{out,j} = \epsilon_k \cdot \sigma \cdot T_k^4 \end{aligned}\]

This can be rewritten in matrix form as:

(2.376)\[\begin{aligned} & \boldsymbol{K} \cdot \boldsymbol{J} = \boldsymbol{E} \end{aligned}\]

where \(K\) is a \(NxN\) matrix, \(J\) is the radiosity vector, and \(E\) is the emissive power. In the above equations, \(T_k\) and the material dependent property \(\epsilon_k\) are calculated for each patch by averaging the element’s values as follows:

(2.377)\[\begin{split}\begin{aligned} & T_k = \frac{ \sum_e A_e \cdot T_e }{ \sum_e A_e } \\ & \epsilon_k = \frac{ \sum_e A_e \cdot \epsilon_e }{ \sum_e A_e } \end{aligned}\end{split}\]

Finally, equation 5 can be solved to obtain the outgoing fluxes \(q_{out,k}\). And the net heat fluxes can be obtained as follows:

(2.378)\[\begin{aligned} & q_{net,k} = q_{out,k} - q_{in,k} = q_{out,k} - \sum_{j=1}^N F_{kj} \cdot q_{out,j} \end{aligned}\]

where Eq. (3) has been used to express the incident fluxes \(q_{in,k}\). The solution process outlined above provides the net heat flux at the element/patch level. This can be applied as boundary condition of the boundary value problem to be solved by the finite elements method after proper transformation of element/patch values to nodal values.