3.4.8. Conjugate Heat Transfer

3.4.8.1. General Formulation

A conjugate heat transfer problem is one in which conductive heat transfer in a solid region is coupled to the convective heat transfer in a neighboring fluid. In its most general form, the coupling at the boundary is governed by the conservation of energy, such that heat flux out of the solid is equal to heat flux into the fluid:

(3.1108) ${\bf q}_s \cdot {\bf n} = {\bf q}_f \cdot {\bf n}$

where ${\bf q}_s$ and ${\bf q}_f$ are the heat flux in the solid and fluid, respectively, and ${\bf n}$ is the surface normal directed {em into} the solid and {em out of} the fluid.

The exact form in which equation ((3.1108)) is implemented depends on whether the fluid flow is laminar or turbulent, since different expressions must be used in these cases for ${\bf q}_f$ . The heat flux in the solid is always due to conduction alone, but there are several possible choices that could be made for the discretization of this flux in space and time.

3.4.8.2. Time Integration

In Fuego, conjugate heat transfer is implemented through loose coupling between the fluid and solid regions, meaning that at each time step, each region is solved separately by treating information from the neighboring region as “given”, and no extra iterations are done between regions to ensure convergence at a single time step. The specific algorithm used can be described as a temperature-forward, flux-back scheme. At a given time step $n$ , the fluid equations are solved using the current solid temperature as a Dirichlet boundary condition; the temperature field of the fluid is thus updated to state $n+1$ everywhere except on the boundary. Then, the heat flux in the fluid at step $n+1$ is computed and transferred to the solid. Finally, the solid region is solved, updating to state $n+1$ using the information from the fluid as a flux boundary condition.

Rather than applying the fluid heat flux directly to the solid, we choose to write the solid boundary condition in the form of a convective heat flux boundary condition:

(3.1109) ${\bf q}_s({\bf x}) \cdot {\bf n} = h({\bf x}) \left( T_\infty({\bf x}) - T_s({\bf x}) \right)$

where $h$ is a convection coefficient, $T_\infty$ is the fluid temperature away from the wall, and $T_s$ is the solid surface temperature. Both $h$ and $T_\infty$ are computed from the fluid temperature field in a way that will be specified, while $T_s$ is left free in the solution of the solid region temperature. This formulation can be shown to be more stable than the alternative of simply transferring the heat flux in the fluid and applying it as a pure Neumann boundary condition to the solid.

Using superscripts to denote time step, the loosely coupled integration scheme can thus be written as:

(3.1110) $T_f^{n+1} = T_s^{n} \mspace{216mu} \text{on $\Gamma_{fs}$} \\$

(3.1111) ${\bf q}_s^{n+1}({\bf x}) \cdot {\bf n} = h^{n+1}({\bf x}) \left( T^{n+1}_\infty({\bf x}) - T^{n+1}_s({\bf x}) \right) \quad \text{on $\Gamma_{fs}$}$

where $\Gamma_{fs}$ is the fluid-solid interface.

3.4.8.3. Discretization of Conduction Region Boundary Condition

The quantity that is needed for a flux boundary condition condition in our CVFEM formulation is the heat flux integrated over the interface surface area associated with each node on the surface. Equation ((3.1111)) is applied to the conduction region at each surface node by assuming that $h$ , $T_s$ and $T_\infty$ can be treated as constants on that node’s sub-control surfaces:

(3.1112) $Q^{n+1}_{s,I} = \int_{SCS_I} {\bf q}_s \cdot {\bf n} dA = h^{n+1}_I A_I \left( T^{n+1}_{\infty,I} - T^{n+1}_{s,I} \right)$

where $A_I$ is the surface area associated with node $I$ . The nodal data $h^{n+1}_I$ and $T^{n+1}_{\infty,I}$ are computed from the fluid solution at time step $n+1$ (see section Computation of Convection Temperature and Coefficient), while $T^{n+1}_{s,I}$ is a degree of freedom solved during the conduction region solution.

3.4.8.4. Computation of Convection Temperature and Coefficient

On the fluid side, the total heat transfer associated with a node on the fluid-solid interface is the integral of the heat flux over that node’s sub-control surfaces on the interface:

(3.1113) $Q^{n+1}_{f,I} = \int_{SCS_I} {\bf q}^{n+1}_f \cdot {\bf n} dA.$

Consider the case in which fluid and solid surfaces meshes conform exactly at the interface. Then, every fluid node can be associated with a corresponding solid node, and using Equations ((3.1108)) and ((3.1112)) we have:

(3.1114) $Q^{n+1}_{f,I} &= Q^{n+1}_{s,I} \\ &= h^{n+1}_I A_I \left( T^{n+1}_{\infty,I} - T^{n+1}_{s,I} \right) \\ &\approx h^{n+1}_I A_I \left( T^{n+1}_{\infty,I} - T^{n+1}_{f,I} \right)$

where the last line (where $T^{n+1}_{f,I}$ is substituted for $T^{n+1}_{s,I}$ ) follows approximately from ((3.1110)); this approximation is of the same order accuracy as the time integration scheme, and for steady state it is exact. In cases in which the surface meshes do not conform exactly, the nodal values of $h_I$ and $T_{\infty,I}$ are passed through an interpolation transfer, introducing a small amount of error.

The total heat transfer $Q_{f,I}$ must be decomposed into two components: $Q_{W,I}$ representing the variables of the fluid at nodes on the surface (“wall”), and $Q_{\infty,I}$ representing variables at nodes away from the surface:

(3.1115) $Q_{f,I} = Q_{W,I} + Q_{\infty,I}$

The way in which this decomposition is done depends on whether the flow is laminar or turbulent, as will be discussed. Comparing this decomposition with (Computation of Convection Temperature and Coefficient), it is clear that:

(3.1116) $Q_{W,I} &= -h_I A_I T_{f,I} \\ Q_{\infty,I} &= h_I A_I T_{\infty,I}$

Rearranging:

(3.1117) $h_I &= -\frac{Q_{W,I}}{T_{f,I} A_I} \\ T_{\infty,I} &= \frac{Q_{\infty,I}}{h_I A_I}$

Finally, we must define the decomposition of $Q_{f,I}$ for laminar and turbulent flow. It is possible when using this approach to end up with negative values for $h_I$ , which appear non-physical to the analyst and are detrimental to the numerical stability of the conduction solve since they reduce diagonal dominance of the linear system. Since the choice of these parameters is arbitrary as long as they reproduce the correct energy flux, when this occurs we reverse the sign of $h_I$ and re-compute $T_{\infty,I}$ as

(3.1118) $T_{\infty,I} = \frac{Q_{f,I}}{h_I} + T_{f,I}$

3.4.8.4.1. Resolution of Boundary Layer

The fluid velocity at the solid surface is zero for laminar flows and turbulent flow models in which the boundary layer is resolved, so all heat transfer in the fluid near walls is due to conduction:

(3.1119) ${\bf q}_f({\bf x}) = -\kappa_f({\bf x}) {\bf \nabla} T({\bf x})$

where $\kappa_f$ is the thermal conductivity of the fluid. Substituting this into ((3.1113)) and using the finite element interpolation for $T({\bf x})$ gives:

(3.1120) $Q_{W,I} = -\int_{SCS_I} \kappa_f \sum_J ({\bf n} \cdot {\bf \nabla} N_J) T_J dA$

where $N_J$ and $T_J$ are respectively the FEM shape function and temperature degree of freedom associated with node $J$ .

The most obvious way of decomposing $Q_{f,I}$ is by breaking the summation into two summations, one over boundary nodes, one over off-boundary nodes:

(3.1121) $Q_{f,I} &= -\int_{SCS_I} \kappa_f \sum_{J \in B} ({\bf n} \cdot {\bf \nabla} N_J) T_J dA \\ Q_{\infty,I} &= -\int_{SCS_I} \kappa_f \sum_{J \notin B} ({\bf n} \cdot {\bf \nabla} N_J) T_J dA$

where $B$ is the set of nodes on the wall.

These quantities, when substituted into (Computation of Convection Temperature and Coefficient), give the computed values of $h_I$ and $T_{\infty,I}$ .

3.4.8.4.2. Turbulent flow modeling with wall functions

In turbulent flow where the boundary layer is not resolved, wall boundary conditions are applied by assuming that the first layer of nodes in the fluid lies not exactly on the solid interface, but slightly away from the wall in the turbulent boundary layer. Various laws of the wall can then be used to relate quantities at these nodes to the wall values.

(3.1122) $Q_{f,I} = c_I A_I (H_I - H_{W,I})$

where $H_I$ is the nodal enthalpy, $H_{W,I}$ is the corresponding enthalpy exactly at the wall, and $c_I$ is a coefficient that depends on the flow variables. The most obvious decomposition is to let $Q_{W,I}=-c_I A_I H_{W,I}$ and $Q_{\infty,I}=c_I A_I H_I$ . However, this most obvious decomposition is incorrect. The difficulty is that enthalpy is measured on a relative scale, rather than an absolute scale like temperature. For example, consider the case where $H_I=0$ . This does not imply that $T_I=0$ ; in Fuego, it usually corresponds to something near standard temperature and pressure. However, the obvious decomposition when substituted into (Computation of Convection Temperature and Coefficient) gives $T_{\infty,I}=0$ , which is clearly the wrong value for the conduction region boundary condition.

Thus, we should choose a decomposition that has $Q_{\infty,I}=0$ only if $T_{\infty,I}$ should be zero. The correct choice is:

(3.1123) $Q_{W,I} &= -c_I T_{W,I} \left (\frac{H_I - H_{W,I}}{T_I - T_{W,I}} \right) \\ Q_{\infty,I} &= c_I T_I \left( \frac{H_I - H_{W,I}}{T_I - T_{W,I}} \right)$

where $T_{W,I}$ is the wall temperature (which for conjugate heat transfer has been obtained from the solid at the previous time step), and $T_I$ is the temperature value at node $I$ (slightly away from the wall). These expressions are undefined if $T_I = T_{W,I}$ ; in that case, the fraction $\Delta H / \Delta T$ is approximated using the limiting value given by the specific heat $c_p$ .