Heat Transfer Theory

Assumptions

1. Equations

heat equation
\[\rho C_p \frac{\partial T}{\partial t} - \nabla \cdot \left( k \nabla T \right) = Q, \quad \text{ in } \Omega\]

which is completed with boundary conditions and initial value

\[\text{at } t=0, \quad T(x,0) = T_0(x)\]

1.1. Convective heat transfer

convective heat equation
\[\rho C_p \left( \frac{\partial T}{\partial t} + \boldsymbol{u} \cdot \nabla T \right) - \nabla \cdot \left( k \nabla T \right) = Q, \quad \text{ in } \Omega\]

1.2. Steady case

steady heat equation
\[ - \nabla \cdot \left( k \nabla T \right) = Q, \quad \text{ in } \Omega\]
steady convective heat equation
\[\rho C_p \boldsymbol{u} \cdot \nabla T - \nabla \cdot \left( k \nabla T \right) = Q, \quad \text{ in } \Omega\]

1.3. Multi-materials

Given a domain \(\Omega \subset \mathbb{R}^d, d=1,2,3\), \(\Omega\) is partitioned into \(N_r\) regions \(\Omega_i,i=1,\ldots,N_r\) corresponding to different materials (solid or fluid). We consider \(\rho_i\), \(C_{p,i}\) and \(k_i\) the material properties defined in each regions \(\Omega_i\). We define also \(\boldsymbol{n}_i\) the outward unit normal vector associated to the boundary \(\partial \Omega_i\).

\[\begin{eqnarray} \rho_i C_{p,i} \frac{\partial T}{\partial t} - \nabla \cdot \left( k_i \nabla T \right) &=& Q, \quad &\text{ in }& \Omega_i,i=1,\ldots,N_r \\ T_{|_{\Omega_i}} &=& T_{|_{\Omega_j}}, \quad &\text{ on }& \partial \Omega_i \cap \Omega_j = \Gamma_{ ij}, \forall i \neq j \\ -k_i \nabla T \cdot \boldsymbol{n}_i &=& k_j \nabla T \cdot \boldsymbol{n}_j, \quad &\text{ on }& \partial \Omega_i \cap \Omega_j = \Gamma_{ ij}, \forall i \neq j \\ \end{eqnarray}\]

We assume the operator \(\mathcal{L}\) tel que \(\mathcal{L} T = \rho_i C_{p,i} \frac{\partial T}{\partial t} - \nabla \cdot \left( k_i \nabla T \right)\) is elliptical.

We multiply \(\mathcal{L} u = Q\) by a function test \(v \in \mathbf{V}\) and integrates by part on \(\Omega_i\). Which give:

\[\rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{\partial T}{\partial t} v - \int_{\Omega_i} \nabla \cdot \left[ k_i \nabla T \right] v = \int_{\Omega_i} Qv, \quad \forall v \in H^1_0(\Omega) \quad for i = 1, \cdots , N_r\]

By the formula of Green, we get

\[\rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{\partial T}{\partial t} v + \int_{\Omega_i} k_i(y) \nabla T \cdot \nabla v- \int_{\partial \Omega} k_i \nabla T \cdot \boldsymbol{n}_i v = \int_{\Omega_i} Qv \quad \forall v \in \mathbf{V}\]

Additivity of the integral, we have

\[\sum_{ i=1}^{N_r} \left( \rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{\partial T}{\partial t} v + \int_{\Omega_i} k_i \nabla T \cdot \nabla v- \int_{\partial \Omega_i} k_i \nabla T \cdot \boldsymbol{n}_i v \right) = \sum_{ i=1}^{N_r} \left( \int_{\Omega_i} Qv \right) \forall v \in \mathbf{V}\]

Note that

\[\bigcup_{ i=1}^{ N } \partial \Omega_i = \bigcup_{ i,j} \Gamma_{ ij} \cup \partial \Omega\]

Use the conditions in the interfaces, we get

\[\sum_{ i=1}^{N_r} \left( \rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{\partial T}{\partial t} v + \int_{\Omega_i} k_i \nabla T \cdot \nabla v- \int_{\partial \Omega} k_i \nabla T \cdot \boldsymbol{n} v \right) = \sum_{ i=1}^{N_r} \left( \int_{\Omega_i} Qv \right) \forall v \in \mathbf{V}\]

Using the implicit Euler method for the time term:

\[\frac{\partial T}{\partial t} (t^{ k+1}) \approx \frac{ T (t^{ k+1}) - T(t^k)}{ dt} \quad \forall t^k \in \mathbb{ R^+} \text{ et } k \in \mathbb{N}\]

Denoting \(T^k = T(t^k)\), we write the formula in \(t^{ k+1}\), we obtain:

\[\sum_{ i=1}^{N_r} \left( \rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{ T^{k+1}}{dt} v + \int_{\Omega_i} k_i \nabla T^{k+1} \cdot \nabla v - \int_{\partial \Omega} k_i \nabla T^{k+1} \cdot \boldsymbol{n} v \right) = \sum_{ i=1}^{N_r} \left( \int_{\Omega_i} \frac{T^{k}}{dt} v + \int_{\Omega_i} Qv \right) \quad \forall v \in \mathbf{V}\]

So, the weak wording becomes:

The weak formulation
\[\text{ On cherche } T \in \mathbf{H} \text{ telle que:} \\ a(T^{k+1}, v) = l(v) \quad \forall v \in \mathbf{V} \\ \text{ and} \quad a(T^{k+1}, v) = \sum_{ i=1}^{N_r} \left( \rho C_{p,i} \displaystyle \int_{\Omega_i} \frac{ T^{k+1}}{dt} v + \int_{\Omega_i} k_i \nabla T^{k+1} \cdot \nabla v - \int_{\partial \Omega} k_i \nabla T^{k+1} \cdot \boldsymbol{n} v \right) \\ l(v) = \sum_{ i=1}^{N_r} \left( \int_{\Omega_i} \frac{T^{k}}{dt} v + \int_{\Omega_i} Qv \right)\]

So we have \(a(u_{k+1},v)\) a continuous bilinear form coercive in \(v \in \mathbf{V}\) and \(l(\phi)\) a continuous linear form . We are in a Hilbert space, so we have all the conditions for the application of the Lax-Milgram theorem. So this problem is well posed.

Correct approximation:

We use the Galerkin approximation method:

Let \(\{ \mathcal{T}_h \}\) a family of meshes of \(:\Omega\).

Let \(\{ \mathcal{K}, P, \sum \}\) a finite element of Lagrange of reference of the degree \(k \geq 1\).

Let \(P^k_{c,h}\) the conforming approximation space defined by

\[P^k_{ c,h} = \{ v \in C^0(\Omega), \forall \mathcal{K} \in \mathcal{T}_h, v|_{\mathcal{K}} \in \mathbb{P}_k(\mathcal{K}) \}\]

To obtain a conformal approximation in V, we add the boundary conditions

\[V_h = P^k_{c,h} \cap V\]

Discrete problem is written:

Problème discrète
\[\text{ Find } T_h \in V_h \text{ such that} \\ a(T_h, v_h) = l(v_h) \quad \forall v_h \in V_h\]

Let \(\{ \varphi_1, \varphi_2, ..., \varphi_N \}\) the base of \(V_h\). An element \(T_h \in V_h\) is written as

\[T_h = \sum^{N}_{l=1} T_l \varphi_l\]

Using \(v\) as a basic function of \(V_h\), our problem becomes

\[\sum_{ i=1}^{N_r} \left( \rho C_{p,i} \displaystyle \int_{\Omega_{i}} \sum_{ l=1}^N T^{k+1}_l \frac{ \varphi_l }{dt} \varphi_j + \int_{\Omega_i} k_i \sum_{ l=1}^N T^{k+1}_l \nabla \varphi_l \cdot \nabla \varphi_j - \int_{\partial \Omega} k_i \sum_{ l=1}^N T^{k+1}_l \nabla \varphi_l \cdot \boldsymbol{n} \varphi_j \right) = \sum_{ i=1}^{N_r} \left( \int_{\Omega_i} \sum_{ l=1}^N T^{k}_l \frac{ \varphi_l }{dt} \varphi_j + \int_{\Omega_i} Q \varphi_j \right)\]

The variational problem of approximation is then equivalent to a linear system

Algebraic problem
\[\text{Determine } T_l \text{ satisfying} \\ \sum_{ l=1}^N a(\varphi_l, \varphi_j) T^{k+1}_l = l(\varphi_j) \forall j = 1, \cdots , N\]

Introduce

\[A = (a(\varphi_i , \varphi_j)), \quad 1 \leq i,j \leq N , \\ U^{k+1} = (T_1^{k+1}, T_2^{k+1}, ..., T_N^{k+1}) \in \mathbb{R}^{N}, \\ F = (l(\varphi_1), l(\varphi_2), ..., l(\varphi_N)) \in \mathbb{R}^{N}\]

We write the system in matrix form

\[AU = F\]

1.4. Variational formulation and discretization of the heat equation with radiative boundary conditions on several surfaces

Radiative heat transfer is not yet available in the toolboxes. An application implementing radiative heat is currently available in feelpp/doc/manual/heat.

Let \(\partial\Omega_D\) and \(\partial \Omega_N\) be the portions of the boundary where Dirichlet and Neumann boundary conditions are applied, respectively. Let us write the variational formulation of the heat equation: find \(T \in H^1((0,T);H^1(\Omega))\) such that, for all \(\phi \in H^1_{0,\partial \Omega_D}(\Omega)\)

\[\int_\Omega \rho c_p \partial_t T \phi + \int_\Omega k \nabla T \cdot \nabla \phi + \int_{\partial \Omega_R} k \nabla T \cdot \vec{n} \phi = \int_\Omega S \phi - \int_{\partial \Omega_N} k \nabla T \cdot \vec{n} \phi.\]

When the radiative boundary \(\partial \Omega_R\) is composed of several subsurfaces that can exchange heat through radiation, the associated radiative boundary condition is complex. In fact, each surface receives heat contributions from the other ones, proportionally to the values of the corresponding view factors. From Modest’s book [Radiative_heat_transfer], equation (5.28), the radiative heat flux at point \(x\) of \(\partial \Omega_R\) is

\[\frac{q(x)}{\epsilon(x)} - \int_{\partial \Omega_R} (\frac{1}{\epsilon(x')}-1)q(x')dF_{x-x'} = E_b(x) - \int_{\partial \Omega_R} E_b(x') dF_{x-x'}.\]

In this equation, \(q=\nabla T \cdot \vec{n}, E_b(x)=\sigma T^4\), \(\epsilon(x)\) is the emittance and \(dF_{x-x'}\) is the view factor between the infinitesimal areas surrounding points \(x,x'\).

1.4.1. Picard’s iteration

Let us now propose a variational formulation for this equation. Let \(\psi \in \mathbb{P}^0_{d}(\partial \Omega_R)\) be discontinuous, piecewise constant basis functions. Functions \(\psi\) are elementwise discontinuous; however, one could choose to work, as a first approximation, with \(\Psi \in P\), where \(P\) a space of discontinuous, piecewise constant basis functions, where discontinuities are not elementwise, but for example discontinuous in correspondence of different radiating surfaces. In the following, we will use \(\psi\) to denote test functions and \(N_h\) to denote the cardinality of the test space. We have

\[\begin{multline} \int_{\partial \Omega_R} \frac{q(x)}{\epsilon(x)} \psi_j - \int_{\partial \Omega_R}\int_{\partial \Omega_R} (\frac{1}{\epsilon(x')}-1)q(x')dF_{x-x'} \psi_j \\ = \int_{\partial \Omega_R} E_b(x) \psi_j- \int_{\partial \Omega_R} \int_{\partial \Omega_R} E_b(x') dF_{x-x'} \psi_j. \end{multline}\]

By decomposing \(q(x)=\sum_{i=0}^{N_h} q_i \psi_i(x)\), the first term of the left-hand side gives rise to a mass matrix \(M_{ij} = \int_{\partial \Omega_R} \epsilon_i^{-1} \psi_i\psi_j\). The non-linearity in \(q\) is handled iteratively: the second term on the left-hand side is treated explicitly and moved to the right-hand side. It can be decomposed as a function of coefficients \(q_i\) as \(N_{j}(q) = \int_{\partial \Omega_R} \Big( \sum_{A_k} \frac{1}{A_k} \int_{A_k} (\frac{1}{\epsilon_k}-1) q_k F_{ik} \, dA \Big)\psi_j\). The right-hand side is of the form \(D_j(T) = \int_{\partial \Omega_R} ( \sigma T^4 + \sum_{A_k} \frac{1}{A_k} \int_{A_k} T^4_k F_{jk} \, dA) \psi_j\).

Due to the presence of the fourth power of the temperature and the non-linearity of the equation with respect to temperature and heat flux, two nested iterative loops are proposed. In the following algorithm, \(n\) denotes time indices, \(k\) denotes the indices of the temperature loop and \(l\) denotes the indices of the flux loop.

Double iterative loop
\begin{algorithm}
    \caption{Solution of the heat equation with radiative BC on the timestep $\Delta t$.}
    \begin{algorithmic}
    \STATE $T^{n}=T^{n,0}$, $q^{n}=q^{n,0,0}$
    \STATE $T^{n+1,0} = Heat_{\Delta t}(T^n,q^{n})$;
    \STATE $q^{n+1,0,0} = q^{n}$
    \WHILE{$||T^{n+1,k+1} -T^{n+1,k}||/||T^{n+1,0}|| > \tau_T$}
        \STATE $T^{n+1,k} \leftarrow T^{n+1,k+1}$
        \STATE $q^{n+1,k,0} = M_{ij}^{-1} (D_j(T^{n+1,k})-N_{ij}(q^{n+1,k+1,0}))$
        \WHILE{$||q^{n+1,k,l} -q^{n+1,k,l-1}||/||q^{n+1,k,0}|| > \tau_q$}
            \STATE $q^{n+1,k,l-1} \leftarrow q^{n+1,k,l}$
            \STATE $q^{n+1,k,l} = M_{ij}^{-1} (D_j(T^{n+1,k})-N_{j}(q^{n+1,k,l-1}))$
        \ENDWHILE
        \STATE $q^{n+1,k+1,0} \leftarrow q^{n+1,k,l}$
        \STATE $T^{n+1,k+1} = Heat_{\Delta t}(T^{n+1,k},q^{n+1,k+1,0})$
    \ENDWHILE
    \STATE $T^{n+1} \leftarrow T^{n+1,k+1}$
    \STATE $q^{n+1} \leftarrow q^{n+1,k+1,0}$
    \end{algorithmic}
\end{algorithm}

1.4.2. Newton non-linear iteration

The previous integral equation can be included in the variational formulation of the PDE as a boundary condition on the heat flux \(q(x)=\Big(-k \nabla T \cdot n \Big)(x)\):

\[q(x) = \epsilon(x) \Bigg( E_b(x) - \int_{\partial \Omega_R} E_b(x') dF_{x-x'} + \int_{\partial \Omega_R} \Big( \frac{1}{\epsilon(x')}-1\Big)q(x') dF_{x-x'} \Bigg).\]

The linearization of the boundary condition generates the following terms to be added to the Jacobian matrix \(DF(T)\), at iteration \(k\),

\[DF(T_k)\alpha = \int_{\partial \Omega_R} 4\sigma \epsilon T_k^3 \alpha \varphi.\]

The following term is added to the residual \(F(T)\) at iteration \(k\), where we denoted by \(A\) and \(A'\) the areas surrounding points \(x\) and \(x'\).

\[\begin{multline} F(T_k) = \int_A \epsilon(x) \sigma T^4_k \varphi - \int_A \epsilon(x) \Big[\int_{A'} \sigma T^4_k(x') dF_{dA-dA'}\Big] \varphi(x) +\\ + \int_A \epsilon(x) \Big[\int_{A'} \Big( \frac{1}{\epsilon(x')} -1\Big)(-k\nabla T_k \cdot \vec{n}) dF_{dA-dA'} \Big]\varphi(x). \end{multline}\]

In practice, the surface of the cavity is divided into smaller surfaces and the view factors \(F_{ij}\) are often already computed and stored for all pairs \(ij\) of (finite) subsurfaces. In that case, the differential \(dF_{dA-dA'}\) can be substituted by the ratio \(F_{ji}/|A_i|\). The meaning of \(dF_{dA-dA'}\) is ``the ratio (diffuse energy leaving $dA$ and intercepted by \(dA'\))/(total energy leaving \(dA\)) '', and the ratio \(F_{ji}/|A_i|\) corresponds to the ratio of (energy leaving surface $j$ and intercepted by surface \(i\))/(total energy leaving surface \(j\)) divided by the area of the intercepting surface. The idea behind the substitution is that integrals \( \Big(\int_{A'} \sigma T^4_k(x') dF_{dA-dA'}\Big) \) and \(\Big(\int_{A'} \Big( \frac{1}{\epsilon(x')} -1\Big)(-k\nabla T_k \cdot \vec{n}) dF_{dA-dA'} \Big)\) produce values for all \(x\in A\), by attributing them a portion of incident radiation which is proportional to the view factor. By dividing \(F_{ji}\) by \(|A_i|\) we construct a density for the view factor, constant over the whole receiving surface, to average the incoming radiation. We hence make the approximation that \( dF_{dA-dA'} = \cos(\theta_i)\cos(\theta_j)/(\pi S^2)dA' \approx F_{ji}/|A_i| dA' = F_{ij}/|A_j| dA' \).

2. Bibliography

[Radiative_heat_transfer] Modest, M.F., Radiative Heat Transfer, Elsevier Science (2013) doi.org/10.1016/C2010-0-65874-3