Radiation coupled with convection

Figures 1. a) Triangular prism. b) Heat flow rate, \(\varphi\), balance on a surface: \(cv\) - convection; \(cd\) - conduction; \(SW\) - short wave radiation; \(LW\) - long wave radiation.

Problem

Let’s consider a triangular prism of infinite length (Figure 1a) with the characteristics (Table 1):

• surface 0: transparent in short-wave, opaque and absorbant in long-wave;

• surfaces 1 and 2: opaque and absorbant in both short-wave and in long-wave.

Table 1. Transmittance \(\tau_{SW} \equiv \tau\), absorptance \(\alpha_{SW} \equiv \alpha\), and emissivity \(\alpha_{LW} \equiv \varepsilon\) of surfaces in short wave (SW) and long wave (LW) radiation. All surfaces are opaque in long-wave, \(\tau_{LW} = 0\)

Surface	𝜏	𝛼	ε
0	1	0	0.8
1 & 2	0	0.8	0.8

The short wave (solar) radiation received by surfaces 1 and 2 are \(E^o_1 = 120 \, \mathrm{W/m^2}\) and \(E^o_2 = 100 \, \mathrm{W/m^2},\) respectively.

The indoor temperature of the prism is controlled at \(T_i = 20\,\mathrm{°C}\). The heat convection coefficient between the walls and the indoor air is \(h_i = 5 \, \mathrm{W/(m^2·K)}\) for all walls.

The temperature of each wall is considered homogeneous. The heat balance of a wall is represented in Figure 1b. The heat flow rate by conduction is negligible for all walls, \(\varphi_{cd} = 0\) (Figure 1b).

In this conditions, find the temperatures of the walls.

Note: Stefan-Bolzmann contant, \(\sigma = 5.670374419...×10^{−8} \mathrm{W⋅m^{−2}⋅K^{−4}}\).

Solution

Problem analysis

Physical model

The indoor air temperature, \(T_i\), is controlled; since it does not change as a result of what happens in the system, it is an input.

The short-wave (solar) irradiances received by the surfaces, \(E^o\), are also inputs.

The temperatures of the surfaces, \(\theta\), are three dependent variables represented by temperature nodes; they are the outputs (dependent variables).

Computational model

The known variables (the inputs of the computational model) are the inputs of the physical model and the outputs of the computational model are the outputs of the physical model. It is a direct problem.

The problem is described by:

1. Vector of outputs

   • \(\theta\) : vector of unknown tempertaures (output of the model).

2. Vectors of variables:

   • \(E^o\): direct short-wave irradiances of the surfaces, W/m².

   • \(E\): total short-wave irradiances of the surfaces, i.e., the sum of direct irradiance and reflected irradiances received from the other surfaces, W/m².

   • \(M^o\): radiant emmitance or total emmisive power of a black body, W/m².

   • \(T_a\): air temperatures beyond the boundary layer of each surface, °C.

3. Matrices of coefficients:

   • Short-wave radiation:

     • \(\alpha\): coefficients of absorption (diagonal matrix), \(0 \le \alpha_i \le 1\),

     • \(\rho\): coefficients of reflexion (diagonal matrix), \(0 \le \rho_i \le 1\),

     • \(\tau\): coefficients of transmittance (diagonal matrix), \(0 \le \tau_i \le 1\),

     • \(F\): view factors, \(0 \le F_{i,j} \le 1\).

   • Long-wave radiation:

     • \(\varepsilon\): hemispherical emissivity , \(0 \le \varepsilon_i \le 1\),

     • \(F\): view factors (same as for short-wave radiation).

   • Convection

     • \(h\): coefficients of heat convection (diagonal matrix), W/(m²·K).

Note:

According to Kirchhoff’s law of thermal radiation, the absorptivity and the emmisivity are equal for the same wavelength, \(\alpha_{\lambda} = \varepsilon_{\lambda}\). Therefore, in order to simplify the notations, absorptivity \(\alpha\) is used for short-wave radiation and emissivity \(\varepsilon\) is used for long-wave radiation.

import numpy as np

np.set_printoptions(precision=1)

# Stefan-Bolzmann constant
σ = 5.670e-8    # W/(m²·K⁴)

# Vectors of inputs
# =================
# Short-wave irradiances of the surfaces
E0 = np.zeros(3)
E0[1], E0[2] = 120, 100     # W/m²

# Air temperature beyond the boundary layer
Ta = 20                     # °C, indoor air temperature
Ta = Ta * np.ones(3)        # the same for all 3 surfaces

# Matrices of coefficients
# ========================
# Transmittance - short-wave (LW)
τ = np.zeros(3)
τ[0] = 1
τ = np.diag(τ)

# Absorptance - short-wave (SW)
α = np.zeros(3)
α[1] = α[2] = 0.8
α = np.diag(α)

# Reflexion - short-wave (SW)
ρ = 1 - (α + τ)

# Emissivity - long-wave (LW)
ε = 0.8 * np.eye(3)

# Heat convection coefficient
h = 5       # W/(m²·K)
h = h * np.eye(3)


# Identity matrix
I = np.eye(3)

Plan for problem solving

The problem to solve is to find the three surface temperatures, \(\theta\). In linear form, it means to solve the equation

\[A_{\theta} \theta = b_{\theta}\]

where:

• \(A_{\theta}\) is a matrix with elements depending on coefficients;

• \(b_{\theta}\) - vector with elements depending on coefficients and inputs;

• \(\theta\) - vector of surface temperatures (size 3).

The system of equations for temperatures will be obtained from the heat balance on the surfaces (Figure 1b):

\[\varphi_{cv} + \varphi_{SW} - \varphi_{LW} + \varphi_{cd} = 0 \]

where:

• \(\varphi_{cv}\) is the flow rate by convection, \(\varphi_{cv} = h (T_a - \theta)\);

• \(\varphi_{SW}\) - short-wave radiation absorbed by the surfaces, \(\varphi_{SW} = A_{E^o}E^o\);

• \(\varphi_{LW}\) - long-wave radiation, \(\varphi_{LW} = A_{M^o}M^o\), where \(M^o_i\) is the vector of emitances \(M^o_i = \sigma T_i^4\);

• \(\varphi_{cd} = 0\) - flow-rate through conduction; zero in this specific problem.

View factors

Since the surfaces are flat, the self-viewing factors are zero: $\(F_{0,0} = F_{1,1} = F_{2,2} \)$

Due to symmetry, \(F_{0,1} = F_{0,2}\). From the summation relation of view factors, \(\sum_{j} F_{0, j} = F_{0, 0} + F_{0, 1} + F_{0, 2} = 1\), it results that \(F_{0, 1} = F_{0, 2} = 1/2\).

From the reciprocity relation, \(F_{0,1} S_0 = F_{1,0} S_1\) and from the relation between surface areas, \(S_0 = \sqrt{2}S_1\), it results that \(F_{1, 0} = \sqrt{2} / 2\).

From the summation relation of view factors, \(\sum_{j} F_{1, j} = F_{1, 0} + F_{1, 1} + F_{1, 2} = 1\), it results that \(F_{1, 2} = 1 - F_{1, 0} = 1 - \sqrt{2} / 2\).

Due to symmetry, \(F_{0,1} = F_{0,2} = 1/2\), \(F_{1,0} = F_{2,0} = \sqrt{2}/2\) and \(F_{1, 2} = F_{2, 1}\) (\(F_{1,2} = 1 - F_{1, 0}\)).

Table 2. View factors.

F	0	1	2
0	0	1/2	1/2
1	$\(\frac{\sqrt{2}}{2}\)$	0	$\(1-\frac{\sqrt{2}}{2}\)$
2	$\(\frac{\sqrt{2}}{2}\)$	$\( 1-\frac{\sqrt{2}}{2}\)$	0

F = np.zeros([3, 3])
F[0, 1] = F[0, 2] = 1 / 2
F[1, 0] = F[2, 0] = np.sqrt(2) / 2
F[1, 2] = F[2, 1] = 1 - F[1, 0]
print("F = ")
print(F)

F = 
[[0.  0.5 0.5]
 [0.7 0.  0.3]
 [0.7 0.3 0. ]]

Short-wave (solar) radiation

The set of equations for short-wave solar radiation is (see Solar radiation abosorbed by the walls in Thermal networks for heat trasfer in buildings):

\[\begin{split} \begin{bmatrix} 1 - \rho_1 F_{1,1} & - \rho_2 F_{1,2} & - \rho_n F_{1,3}\\ - \rho_1 F_{2,1} & 1 - \rho_2 F_{2,2} & - \rho_n F_{2,n3} \\ - \rho_1 F_{n,1} & - \rho_2 F_{n,1} & 1 - \rho_n F_{n,n} \end{bmatrix} \begin{bmatrix} E_1\\ E_2\\ E_3 \end{bmatrix} = \begin{bmatrix} E_{1}^{o}\\ E_{2}^{o}\\ E_{3}^{o} \end{bmatrix} \end{split}\]

or

\[(I - \mathrm{diag}(\rho) \circ F) E = E^o\]

The unknown total short-wave irradiances on walls, in W/m², are then

\[ E = (I - \mathrm{diag}(\rho) \circ F)^{-1} E^o\]

where:

• \(\mathrm{diag}(·)\) is the matrix-to-vector diagonal operator;

• the symbol \(\circ\) represents the Hadamard (or element-wise) product;

\(I =\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \) is the identity matrix;

\(\rho = \begin{bmatrix} \rho_1 & 0 & 0\\ 0 & \rho_2 & 0\\ 0 & 0 & \rho_3 \end{bmatrix}\) - diagonal matrix of reflectances, \(0 \le \rho_i \le 1\);

\(F = \begin{bmatrix} F_{1,1} & F_{1,2} & F_{1,3}\\ F_{2,1} & F_{2,2} & F_{2,3} \\ F_{3,1} & F_{3,2} & F_{3,3} \end{bmatrix}\) - matrix of view factors, \(0 \le F_{i,j} \le 1\);

\(E^o = \begin{bmatrix} E_{1}^{o}\\ E_{2}^{o}\\ E_{3}^{o} \end{bmatrix}\) - vector of direct solar irradiances, W/m²;

\(E = \begin{bmatrix} E_1\\ E_2\\ E_3 \end{bmatrix}\) - vector of unknown total short-wave irradiances, W/m².

E = np.linalg.inv((I - np.diag(ρ) * F)) @ E0
print("E =")
print(E, 'W/m²')

E =
[ 23.4 126.3 107.4] W/m²

Short-wave heat flow rate absorbed by the surfaces is:

\[\varphi_{SW} = \alpha E\]

φSW = np.diag(α) * E
print("φSW = ")
print(φSW, 'W/m²')

φSW = 
[  0.  101.   85.9] W/m²

Long-wave (infrared) thermal radiation

Long-wave absorption coefficients of the surfaces are \(\alpha_{LW,0} = \alpha_{LW,1} = \alpha_{LW,2} = 0.8\).

# α = 0.8 * np.ones(3)
# α = np.diag(α)
# ε = α       # Kirchhoff law

The net flow rate density emmitted by a surface is:

\[\varphi_{LW, i} = \varepsilon_i M^o_i - \alpha_i E_{LW,i}\]

where:

• \(M^o_i\) - emmitance of the surface, \(M^o_i = \sigma T_i^4\);

• \(E_{LW,i}\) - global long wave radiation (illuminace) received by unit area of the surface;

• \(\alpha_i\) - coefficient of absorptivity in long wave radiation;

• \(\varepsilon_i\) - coefficient of hemispherical emmissivity in long wave radiation.

The global long wave radiation received by the surface, \(S_i E_{LW,i}\), is the sum of the fractions \(F_{j,i}\) of the flows emmitted and reflected by other surfaces \(j\) and intercepted by surface \(i\),

\[S_i E_{LW,i} = \sum_j F_{j, i} S_j (\varepsilon_jM^o_j + \rho_j E_{LW, j}) = \sum_j F_{j, i} S_j J_j\]

where:

• \(J_j = \varepsilon_jM^o_j + \rho_j E_{LW, j}\) is the radiosity.

Taking into account the relation of reciprocity, \(F_{i, j} S_i = F_{j, i} S_j\), it results that

\[E_{LW,i} = \sum_j F_{i, j} (\varepsilon_jM^o_j + \rho_j E_{LW, j})\]

For opaque walls, \(\tau_i = 0\). Then \(\rho_i = 1 - \alpha_i\). Since \(\alpha_i = \varepsilon_i\) (Kirchoff’s law of thermal radiation), \(\rho_i = 1 - \varepsilon_i\). Therefore,

\[E_{LW,i} - \sum_j F_{i, j} (1 - \varepsilon_j) E_{LW,j} = \sum_j F_{i,j} \varepsilon_j M^o_j\]

In matrix form, this equation is:

\[ (I - F(I - \varepsilon))E_{LW} = F \varepsilon M^o\]

In our case, the matrix \((I - F(I - \varepsilon))\) is:

print(I - F @ (I - ε))

[[ 1.  -0.1 -0.1]
 [-0.1  1.  -0.1]
 [-0.1 -0.1  1. ]]

and the matrix \(F \, \varepsilon\) is:

print(F @ ε)

[[0.  0.4 0.4]
 [0.6 0.  0.2]
 [0.6 0.2 0. ]]

From this equation, it results the vector of long wave radiation on surfaces:

\[E_{LW} = [I - F(I - \varepsilon)]^{-1} F \varepsilon M^o\]

By substituting \(E_{LW}\) into

\[\varphi_{LW, i} = \varepsilon_i M^o_i - \alpha_i E_{LW,i}\]

we obtain

\[\varphi_{LW} = \varepsilon \{I - [I - F (I - \varepsilon)]^{-1}\ F \varepsilon\} M^o = K M^o\]

In our case, the matrix \(K = \varepsilon \{I - [I - F (I - \varepsilon)]^{-1}\ F \varepsilon\}\) is:

K = ε @ (I - np.linalg.inv(I - F @ (I - ε)) @ F @ ε)
print("K = ")
print(K)

K = 
[[ 0.7 -0.4 -0.4]
 [-0.5  0.7 -0.2]
 [-0.5 -0.2  0.7]]

The equation of \(\varphi_{LW} = K M^o\) is linear in emmitances \(M^o\) but not linear in temperatures, since \(M^o_i = \sigma T_i^4\). In order to linearize the equation \(\varphi_{LW} = K \sigma T^4\), where

\[\begin{split}T^4 = \begin{bmatrix} T_0^4\\ T_1^4\\ T_2^4 \end{bmatrix}\end{split}\]

is the vector of (absolute) temperatures of the surfaces, let’s consider the particular case in which all temperatures have the same value \(\bar{T}\). In this case, the heat flow rate between the surfaces is zero,

\[\bar \varphi_{LW} = K \sigma \bar T^4 = 0\]

Therefore,

\[\varphi_{LW} - \bar \varphi_{LW} = \varphi_{LW} = K \sigma (T^4 - \bar T^4)\]

Taking into account that

\[T_i^4 - \bar T_i^4 = (T_i^2 + \bar T_i^2 )(T_i^2 - \bar T_i^2 ) = (T_i^2 + \bar T_i^2 )(T_i + \bar T_i )(T_i - \bar T_i )\]

the linearized equation for long-wave radiation is then:

\[\varphi_{LW} = K M^o = \sigma K H (T - \bar T) = \sigma K H \theta\]

where the matrix \(H\) is diagonal with elements

\[H_i = (T_i^2 + \bar T_i^2 )(T_i + \bar T_i)\]

Solving this equation requires an initial guess of \(\bar T_i\). If we consider that \(T_i \approx \bar T_i\) and that \(\bar T_i = \bar T = 293.15 \, \mathrm K\) for all \(i\), then

\[T_i^4 - \bar T_i^4 \approx 4 \bar T^3\]

and

\[\varphi_{LW} = K M^o \approx 4 \sigma \bar T^3 K \theta\]

where

\[\begin{split}\theta = \begin{bmatrix} \theta_0\\ \theta_1\\ \theta_2 \end{bmatrix}\end{split}\]

are the surface temperatures expressed in degree Celsius.

For our example, the diagonal matrix \(H \approx 4 \sigma \bar T^3 I\) is:

T0 = 273.15 + 20   # K, temperature
H = 4 * σ * T0**3 * np.eye(3)
print("H =")
print(H)

H =
[[5.7 0.  0. ]
 [0.  5.7 0. ]
 [0.  0.  5.7]]

The matrix \(\sigma \bar T^3 K\) of the linearized expression of heat flow rates

\[\varphi_{LW} = 4 \sigma \bar T^3 K \theta = L \theta\]

is:

L = H @ K
print("L = ")
print(L)

L = 
[[ 4.  -2.  -2. ]
 [-2.8  4.2 -1.4]
 [-2.8 -1.4  4.2]]

Note: For temperatures that are usual in buildings, i.e., \(0°C \leq \bar T - 273.15 \leq 40 °C\), the values of \(4 \sigma \bar T^3\) vary between 4.5 and 7.

T0 = 273.15 + np.array([0, 40])
print(4 * σ * T0**3)

[4.6 7. ]

Coupled short-wave radiation, long-wave radiation and convection

From the heat flow balance for a surface \(i\) shown in Figure 1b, it results that

\[\varphi_{cv} + \varphi_{SW} - \varphi_{LW} + \varphi_{cd} = 0 \]

where:

• \(\varphi_{cv} = h (T_a - \theta)\) is the flow rate by convection;

• \(\varphi_{SW}\) - short-wave radiation absorbed by the surfaces;

• \(\varphi_{LW} = L \theta\) - linearized long-wave radiation;

• \(\varphi_{cd} = 0\) - flow-rate through conduction; it is zero in this specific problem.

By substitution, it results:

\[h (T_a - \theta) + \varphi_{SW} - L \theta = 0 \]

which, by considering that the vector of temperatures is the uknown, becomes:

\[(L + hI) \theta = (\varphi_{SW} + h T_a)\]

which is an equation of the form

\[A_{\theta} \theta = b_{\theta}\]

where:

• \(A_{\theta} = L + hI\) is a matrix with elements depending on coeffcients (\(\alpha\), \(\varepsilon\), \(F\) and \(h\));

• \(b_{\theta} = \varphi_{SW} + h T_a\) - vector with elements depending on coefficients (\(\alpha\), \(\varepsilon\), \(F\) and \(h\)) and inputs (\(E^o\) and \(T_a\));

• \(\theta\) - vector of temperatures (size 3)

By solving for surface temperatures \(\theta\), it results:

\[\theta = (L + hI)^{-1} (\varphi_{SW} + h T_a)\]

θ = np.linalg.inv(L + h * I) @ (φSW + h @ Ta)
print("θ = ")
print(θ, "°C")

θ = 
[26.3 34.9 33.5] °C