Turbulent Heat Transfer Analysis of Silicon Carbide Ceramic Foam as a Solar Volumetric Receiver

Abstract

A volumetric solar receiver receives the concentrated radiation generated by a large number of heliostats. Turbulent heat transfer occurs from the solid matrix to the air as it passes through the porous receiver. Such combined heat transfer within the receiver, including radiation, convection and conduction, is studied using a local thermal non-equilibrium model. Both the Rosseland approximation and the P1 model are applied to consider the radiative heat transfer through the solar receiver. Furthermore, the low Mach approximation is exploited to investigate the compressible flow through the receiver. Analytic solutions are obtained for the developments of air and ceramic temperatures as well as the pressure along the flow direction. Since the corresponding fluid and solid temperature variations generated under the Rosseland approximation agree fairly well with those based on the P1 model, the Rosseland approximation is used for further analysis. The results indicate that the pore diameter must be larger than its critical value to obtain high receiver efficiency. Moreover, it has been found that optimal pore diameter exists for achieving the maximum receiver efficiency under the equal pumping power. The solutions provide effective guidance for a novel volumetric solar receiver design of silicon carbide ceramic foam.

1. Introduction

A solar volumetric receiver is required to have the resistance to temperature as high as 1000 degree Celsius, high porosity for sufficiently large extinction volume such that the concentrated solar radiation penetrates through the receiver, high cell density to achieve large specific surface area and sufficiently high effective thermal conductivity to avoid possible thermal spots. Extruded monoliths with parallel channels (i.e. honeycomb structure) are being used in some solar power plants in Europe, including the solar power tower plant of 1.5 MW built in 2009, in Julich in Germany [1, 2]. However, in such conventional receivers, both thermal spots [3] and flow instabilities [4] have been often reported. In the monolith receiver, locally high solar flux leads to a low mass flow with high temperature, whereas locally low solar flux leads to a high mass flow with low temperature. This causes the absorber material to exceed the design temperature locally, which then leads to its destruction although the average temperature is comparatively low. These difficulties encountered in the receiver must be overcome to run the power plant safely.

In consideration of these requirements, ceramic foams have come to draw attention as a possible candidate to replace the conventional extruded monoliths with parallel channels. Many researchers including Becker et al. [4], Fend et al. [5] and Bai [6] focused on porous ceramic foams as a promising absorber material. Recently, Sano et al. [7] carried out a local non-thermal equilibrium analysis to investigate the receiver efficiency under the equal pumping power. For the first time, the complete set of analytical solutions based on the two-energy equation model of porous media was presented, so as to fully account for the combined effects of tortuosity; thermal dispersion and compressibility on the convective, conductive and radiative heat transfer within a ceramic foam receiver. In their analysis, however, the Rosseland approximation was applied to account for the radiative heat transfer through the solar receiver. It is well known that the Rosseland approximation ceases to be valid near boundaries. Although no wall boundaries exist for the case of the one-dimensional analysis of the solar volumetric receiver, the validity of applying the Rosseland approximation near the inlet boundary of the receiver has not been investigated yet. Furthermore, the effects of turbulence mixing on the heat transfer were not considered.

In this study, the validity of the Rosseland approximation [7] will be examined by comparing the results based on the Rosseland approximation and the results obtained from solving the irradiation transport equation based on the P1 model. The set of the equations will be reduced to a fifth-order ordinary differential equation for the air temperature. Once the air temperature distribution is determined, the pressure distribution along the flow direction can readily be estimated from the momentum equation with the low Mach approximation. Thus, the receiver efficiency, namely, the ratio of the air enthalpy flux increase to the concentrated solar heat flux, can be compared under the equal pumping power, so as to investigate the optimal operating conditions. Some analytical and numerical investigations [3, 4, 5, 6, 7, 8] have been reported elsewhere. However, none of them appeared to elucidate well the combined effects of turbulence, compressibility, radiation, convection and conduction within the volumetric receiver on the developments of air and ceramic temperatures as well as the pressure along the flow direction. This study appears to be the first to provide the complete set of analytical solutions based on the two-energy equation model of porous media [9], fully accounting for the combined effects of turbulence, tortuosity, thermal dispersion, compressibility and radiative heat transfer within a ceramic foam receiver.

2. Volume averaged governing equations

As illustrated in Figure 1, the structure of silicon carbide ceramic foam volumetric receiver may be treated as homogeneous porous medium. Since the dependence of the Darcian velocity on the transverse direction can only be observed in a small region very close to the walls of the passage, we may neglect the boundary effects (i.e. Brinkman term).

Figure 1.

Volumetric receiver.

Based on a theoretical derivation of Darcy’s law, Neuman [10] pointed out that the application of Darcy’s law to compressible fluids is justified as long as Knudsen numbers are sufficiently small to ensure the no-slip conditions at the solid–gas interface. This is usually the case for the volumetric receivers. Thus, allowing the density to vary through the receiver, the following Forchheimer extended Darcy law should hold:

$$-\frac{\partial {⟨p⟩}^f}{\partial x_i}=\frac{{⟨\mu ⟩}^f}{K}⟨u_i⟩+b{⟨\rho ⟩}^f\sqrt{⟨u_j⟩⟨u_j⟩}⟨u_i⟩$$E1

where K and b are the permeability and the inertial coefficients, respectively. Furthermore, by virtue of the volume averaging procedure [11, 12, 13], the microscopic energy equations of the compressible fluid flow phase and the solid phase may be integrated over an elemental control volume V, so as to derive the corresponding macroscopic energy equations. Since the porous medium is considered to be homogeneous, the integration of the two distinct energy equations gives:

For the air:

$$\begin{array}{l}\epsilon \frac{\partial }{\partial t}{\rho }_f{\left(h_{stag}-\frac{p}{\rho }\right)}^f+\epsilon \frac{\partial }{\partial x_j}{\rho }^f{u}_j{}^f{h}_{st}{}^f\hfill \\ =\frac{\partial }{\partial x_j}\left(\epsilon k_f{}^f\frac{\partial T^f}{\partial x_j}+\frac{1}{V}\underset{A_{\mathrm{int}}}{{\displaystyle \int }}{k}_{f}T{n}_{j}dA-\epsilon {\rho }^f{\tilde{h}}_{stag}{\tilde{u}}_j{}^f+\epsilon u_i{\tau }_{ij}{}^f\right)\hfill \\ +\frac{1}{V}\underset{A_{\mathrm{int}}}{{\displaystyle \int }}{k}_f\frac{\partial T}{\partial x_j}{n}_{j}dA\hfill \end{array}$$E2

For the solid matrix:

$$\left(1-\epsilon \right){\rho }_s{c}_s\frac{\partial {⟨T⟩}^s}{\partial t}=\frac{\partial }{\partial x_j}\left(\left(1-\epsilon \right)k_s\frac{\partial {⟨T⟩}^s}{\partial x_j}-\frac{k_s}{V}\underset{A_{\mathrm{int}}}{\int }{\mathit{Tn}}_j\mathit{dA}‐q_R\right)-\frac{1}{V}\underset{A_{\mathrm{int}}}{\int }{k}_f\frac{\partial T}{\partial x_j}{n}_j\mathit{dA}$$E3

where the intrinsic volume average of a certain local variable $\varphi$ in the fluid phase and solid matrix phase can be defined as.

$${⟨\varphi ⟩}^f\equiv \frac{1}{V_f}{\int }_{V_f}\mathit{\varphi dV},{⟨\varphi ⟩}^m\equiv \frac{1}{V_m}{\int }_{V_m}\mathit{\varphi dV}$$E4

Note that subscripts f and m refer to the fluid phase and solid matrix phase, respectively. The decomposition of the local variable $\varphi$ can be expressed in terms of its intrinsic average and the spatial deviation from it:

$$\varphi ={⟨\varphi ⟩}^f+\tilde{\varphi }$$E5

Moreover, $q_R$ is the radiative heat flux, $A_{\mathrm{int}}$ is the interfacial surface area between the fluid and solid matrix phases, while $n_j$ is the normal unit vector from the fluid phase to the solid matrix phase.

In order to simplify the foregoing set of the equations, the low Mach approximation is applied due to the relatively low Mach number when the air flows through a porous medium. Thus, the dynamic pressure change is sufficiently small as compared to the absolute pressure prevailing over the system, such that the stagnant enthalpy is approximated by $h_{\mathit{\text{stag}}}=h+u_k{u}_k/2\cong h$. Combining the foregoing two energy equations namely Eqs. (2) and (3), and, then, noting the continuity of temperature and heat flux at the interface, we obtain the one-equation model for the steady state as follows:

$$\begin{array}{c}\epsilon \frac{\partial }{\partial x_j}{\rho }_f{}^f{u}_j{}^f{h}^f=\frac{\partial }{\partial x_j}(\epsilon k_f{}^f\frac{\partial T^f}{\partial x_j}+\left(1-\epsilon \right)k_s\frac{\partial T^s}{\partial x_j}+\frac{1}{V}\underset{A_{\mathrm{int}}}{{\displaystyle \int }}\left(k_f-k_s\right)T{n}_{j}dA\\ -\epsilon {\rho }_f{}^f\tilde{h}{\tilde{u}}_j{}^f-q_R)\end{array}$$E6

For the time being, let us assume $\partial {⟨T⟩}^f/\partial x_j\cong \partial {⟨T⟩}^s/\partial x_j\cong \partial ⟨T⟩/\partial x_j$ (this assumption will be relaxed shortly). Then, the equation reduces to

$$\begin{array}{c}\epsilon \frac{\partial }{\partial x_j}{\rho }_f{}^f{u}_j{}^f{h}^f=\frac{\partial }{\partial x_j}(\left(\epsilon k_f{}^f+\left(1-\epsilon \right)k_s\right)\frac{\partial T}{\partial x_j}+\frac{1}{V}\underset{A_{\mathrm{int}}}{{\displaystyle \int }}\left(k_f-k_s\right)T{n}_{j}dA\\ -\epsilon {\rho }_f{}^f\tilde{h}{\tilde{u}}_j{}^f-q_R\end{array}$$E7

where

$$⟨\varphi ⟩\equiv \frac{1}{V}{\int }_V\mathit{\varphi dV}$$E8

is the Darcian average of the variable $\varphi$ such that $⟨u_j⟩=\epsilon {⟨u_j⟩}^f$ is the Darcian velocity vector. From the foregoing equation, that is, Eq. (6), the macroscopic heat flux vector $q_i=\left(q_x,q_y,q_z\right)$ and its corresponding stagnant thermal conductivity $k_{\mathit{\text{stag}}}$ may be defined as follows:

$$\begin{array}{c}{q}_i=-k_{stag}\frac{\partial T}{\partial x_i}+q_R+\epsilon {\rho }_f{}^f\tilde{h}{\tilde{u}}_i{}^f\\ =-\left(\epsilon k_f{}^f+\left(1-\epsilon \right)k_s\right)\frac{\partial T}{\partial x_i}-\frac{1}{V}\underset{A_{\mathrm{int}}}{{\displaystyle \int }}\left(k_f-k_s\right)T{n}_{i}dA+q_R+\epsilon {\rho }_f{}^f\tilde{h}{\tilde{u}}_i{}^f\end{array}$$E9

or

$$k_{\mathit{\text{stag}}}\frac{\partial ⟨T⟩}{\partial x_i}=\left(\epsilon {⟨k_f⟩}^f+\left(1-\epsilon \right)k_s\right)\frac{\partial ⟨T⟩}{\partial x_i}+\frac{1}{V}\underset{A_{\mathrm{int}}}{\int }\left(k_f-k_s\right){\mathit{Tn}}_i\mathit{dA}$$E10

The term $\epsilon {⟨{\rho }_f⟩}^f{⟨{\tilde{h}\tilde{u}}_i⟩}^f$ in Eq. (9) describes the thermal dispersion heat flux vector, which serves an additional heat flux resulting from the hydrodynamic mixing of fluid particles passing through pores. On the other hand, the second term on the right-hand side term in Eq. (10) is associated with the surface integral, and it describes the effects of the tortuosity on the macroscopic heat flux, which adjusts the level of the stagnant thermal conductivity from its upper bound $\left(\epsilon k_f+\left(1-\epsilon \right)k_s\right)$ to a correct one. Yang and Nakayama [9] introduced the effective porosity ${\epsilon }^{\ast }$, which is defined as

$${\epsilon }^{\ast }=\frac{k_s-k_{\mathit{\text{stag}}}}{k_s-k_f}=\epsilon +\frac{\epsilon k_f+\left(1-\epsilon \right)k_s-k_{\mathit{\text{stag}}}}{k_s-k_f}$$E11

such that

$$\left({\epsilon }^{\ast }-\epsilon \right)\frac{\partial ⟨T⟩}{\partial x_i}=\frac{1}{V}{\int }_{A_{\mathrm{int}}}{\mathit{Tn}}_i\mathit{dA}$$E12

Using the effective porosity ${\epsilon }^{\ast }$ and the equation of state ${⟨p⟩}^f={⟨{\rho }_f⟩}^{f}R{⟨T⟩}^f$ and ${⟨h⟩}^f=c_p{⟨T⟩}^f$, the volume average energy equations Eqs. (2) and (3) may be concisely rewritten for the steady state for air as:

$$\epsilon c_p\frac{\partial }{\partial x_j}{⟨{\rho }_f⟩}^f{⟨u_j⟩}^f{⟨T⟩}^f=\frac{\partial }{\partial x_j}\left({\epsilon }^{\ast }{⟨k_f⟩}^f\frac{\partial {⟨T⟩}^f}{\partial x_j}+\epsilon k_{{\mathit{dis}}_{\mathit{jk}}}\frac{\partial {⟨T⟩}^f}{\partial x_k}\right)-h_v\left({⟨T⟩}^f-{⟨T⟩}^s\right)$$E13

for the solid matrix phase as:

$$\frac{\partial }{\partial x_j}\left(\left(1-{\epsilon }^{\ast }\right)k_s\frac{\partial {⟨T⟩}^s}{\partial x_j}+q_R\right)-h_v\left({⟨T⟩}^s-{⟨T⟩}^f\right)=0$$E14

Note that the assumption of equal temperature gradients, $\partial {⟨T⟩}^f/\partial x_j\cong \partial {⟨T⟩}^s/\partial x_j\cong \partial ⟨T⟩/\partial x_j$, has been discarded. This practice has been proven to be quite effective in a series of computations (e.g. [8, 9]). According to the gradient diffusion hypothesis [14], the thermal dispersion term is usually expressed as:

$${⟨{\rho }_f⟩}^f{⟨{\tilde{h}\tilde{u}}_j⟩}^f={⟨{\rho }_f⟩}^f{c}_p{⟨{\tilde{T}\tilde{u}}_j⟩}^f=-k_{\mathit{disjk}}\frac{\partial {⟨T⟩}^f}{\partial x_k}$$E15

while the interfacial heat transfer between the solid and fluid phases is modeled using Newton’s cooling law:

$$\frac{1}{V}{\int }_{A_{\mathrm{int}}}{k}_f\frac{\partial T}{\partial x_j}{n}_j\mathit{dA}=h_v\left({⟨T⟩}^s-{⟨T⟩}^f\right)$$E16

where $h_v$ is the volumetric heat transfer coefficient. The Maxwell approximations may be used for the dynamic viscosity and thermal conductivity of the air:

$$\mu \left({⟨T⟩}^f\right)={\mu }_0{\left(\frac{{⟨T⟩}^f}{{⟨T⟩}_0^f}\right)}^n=1.8\times {10}^{-5}{\left(\frac{{⟨T⟩}^f}{300K}\right)}^{0.7}\left[\mathrm{Pa}\cdot \mathrm{s}\right]$$E17

and

$$k_f\left({⟨T⟩}^f\right)=k_0{\left(\frac{{⟨T⟩}^f}{{⟨T⟩}_0^f}\right)}^n=0.025{\left(\frac{{⟨T⟩}^f}{300K}\right)}^{0.7}\left[\mathrm{W}/\mathrm{mK}\right]$$E18

where the exponent n is 0.7 according to [4]. The specific heat capacity of the air $c_p=1000\left[\mathrm{J}/\mathrm{kgK}\right]$ and the Prandtl number $\mathit{Pr}=1.8\times {10}^{-5}\times 1000/0.025=0.72$ are assumed to be constant.

3. One-dimensional analysis for volumetric receiver

In this section, we perform one-dimensional analysis to obtain analytic solutions for convective-radiative heat transfer in volume receiver. Prior to that, the radiative heat flux $q_R$ needs to be determined in advance. In the literature, there are two models, namely, the Rosseland approximation and the P1 model.

3.1 Analysis based on the Rosseland approximation

In the Rosseland approximation, the radiative heat flux is given by

$$q_R=-\frac{16\sigma }{3\beta }{\left({⟨T⟩}^s\right)}^3\frac{\partial {⟨T⟩}^s}{\partial x_j}$$E19

where $\sigma =5.67\times {10}^{-8}\left[\mathrm{W}/{\mathrm{m}}^2{\mathrm{K}}^4\right]$ is the Stephan-Boltzmann constant while $\beta$ is the mean extinction coefficient.

As schematically shown in Figure 1, the air is flowing through a passage of length L at the rate of the mass flux $G={⟨{\rho }_f⟩}^f⟨u⟩$. Under the low Mach number approximation, namely, ${⟨{\rho }_f⟩}^f\propto 1/{⟨T⟩}^f$, the macroscopic governing equations Eqs. (1), (13) and (14) can be simplified to be a one-dimensional set of equations as follows:

$$-\frac{d{⟨p⟩}^f}{\mathit{dx}}=\frac{{⟨\mu ⟩}^f}{K}\frac{G}{{⟨{\rho }_f⟩}^f}+b\frac{G^2}{{⟨{\rho }_f⟩}^f}=\frac{R}{{⟨p⟩}^f}\left(\frac{{⟨\mu ⟩}^f}{K}G+{\mathit{bG}}^2\right){⟨T⟩}^f$$E20

$$c_{p}G\frac{d{⟨T⟩}^f}{\mathit{dx}}=\frac{d}{\mathit{dx}}\left({\epsilon }^{\ast }{⟨k_f⟩}^f+\epsilon k_{\mathit{disxx}}\right)\frac{d{⟨T⟩}^f}{\mathit{dx}}-h_v\left({⟨T⟩}^f-{⟨T⟩}^s\right)$$E21

$$\frac{d}{\mathit{dx}}\left(\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left({⟨T⟩}^s\right)}^3\right)\frac{d{⟨T⟩}^s}{\mathit{dx}}-h_v\left({⟨T⟩}^s-{⟨T⟩}^f\right)=0$$E22

According to Calmidi and Mahajan [15, 16], Dukhan [17], Kuwahara et al. [18] and Yang et al. [19, 20], the permeability and inertial coefficient of foams are given by

$$K=0.00073{\left(1-\epsilon \right)}^{-0.224}{\left(\frac{1.18}{1-e^{-\left(1-\epsilon \right)/0.04}}\sqrt{\frac{1-\epsilon }{3\pi }}\right)}^{-1.11}{d_m}^2$$E23

and

$$b=\frac{12\left(1-\epsilon \right)}{d_m}$$E24

respectively, where $d_m$ is the pore diameter of foam. The longitudinal dispersion coefficient is roughly about 20 times more than the transverse one. Thus, following Calmidi and Mahajan [16], we may evaluate the longitudinal dispersion coefficient using the following expression:

$$\epsilon k_{{\mathit{dis}}_{\mathit{xx}}}=1.2c_{p}G\sqrt{K}$$E25

With respect to the stagnant thermal conductivity and the volumetric heat transfer coefficient for foams, Calmidi and Mahajan [15, 16] empirically provided the following correlations:

$$k_{\mathit{\text{stag}}}=\epsilon k_f+0.19{\left(1-\epsilon \right)}^{0.763}{k}_s$$E26

$$h_v=8.72{\left(1-\epsilon \right)}^{1/4}{\left(\frac{1-e^{-\left(1-\epsilon \right)/0.04}}{\epsilon }\right)}^{1/2}{\left(\frac{{\mathit{Gd}}_m}{{⟨\mu ⟩}^f}\right)}^{1/2}{\mathit{Pr}}^{0.37}\frac{k_f}{{d_m}^2}$$E27

Kamiuto et al. [21] experimentally affirmed that the Rosseland model is quite effective. Therefore, it can be deduced that the Rosseland model is also applicable for the present case of silicon carbide ceramic foam. Based on the measurements made on cordierite ceramic foams by Kamiuto et al., the mean extinction coefficient $\beta$ is calculated by the following correlation:

$$\beta =8\left(1-\epsilon \right)/d_m$$E28

For a given mass flux $G={⟨{\rho }_f⟩}^f⟨u⟩$, the foregoing three equations along with the equation of state may be solved for the four unknowns, namely, ${⟨T⟩}^f$,${⟨T⟩}^s$, ${⟨p⟩}^f$ and ${⟨{\rho }_f⟩}^f$. The boundary conditions are given as follows:

$x=0$ (inlet):

$${⟨T⟩}^f={⟨T⟩}_0^f=300\left[\mathrm{K}\right]$$E29

$${⟨p⟩}^f={⟨p⟩}_0^f={10}^5\left[\mathrm{P}a\right]$$E30

such that ${⟨{\rho }_f⟩}^f={⟨{\rho }_f⟩}_0^f={⟨p⟩}_0^f/R{⟨T⟩}_0^f={10}^5/\left(287\times 300\right)=1.2\left[\mathrm{kg}/{\mathrm{m}}^3\right]$

$$\begin{array}{l}-\left({\epsilon }^{\ast }{k}_f{}^f+\epsilon k_{disxx}\right){\frac{d{T}^f}{dx}|}_{x=0}-\left(\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left(T^s\right)}^3\right){\frac{d{T}^s}{dx}|}_{x=0}\hfill \\ =I_0\cos \xi -\left(1-\epsilon \right)\left(a\sigma \left({\left(T_0^s\right)}^4-{\left(T_0^f\right)}^4\right)+h_{conv}\left(T_0^s-T_0^f\right)\right)\hfill \end{array}$$E31

where $I_0$ is the intensity of radiation and $\xi$ is the incidence angle. Moreover, $a\cong 0.9$ is the emissivity of the front surface of the receiver, while $h_{\mathit{\text{conv}}}$ is the convective heat transfer coefficient at the frontal surface. The properties of the air depend on the temperature, which makes the integrations of the foregoing governing equations formidable. In order to obtain analytic expressions for the unknown variables, we may approximate these properties by their representative values evaluated at the average air temperature over the receiver as given by

$$\overline{{⟨T⟩}^f}=\frac{1}{L}{\int }_0^L{⟨T⟩}^f\mathit{dx}$$E32

Likewise we shall define the solid phase average temperature as follows:

$$\overline{{⟨T⟩}^s}=\frac{1}{L}{\int }_0^L{⟨T⟩}^s\mathit{dx}$$E33

The two energy equations, that is, Eqs. (21) and (22) may be added together and integrated using the boundary conditions in Eqs. (29) and (31) to give

$$\begin{array}{c}{c}_{p}G\left(T^f-T_0^f\right)=\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{T^f}}{T_0^f}\right)}^n+\epsilon k_{disxx}\right)\frac{d{T}^f}{dx}\\ +\left(\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left(\overline{T^s}\right)}^3\right)\frac{d{T}^s}{dx}+I_0\cos \xi \\ -\left(1-\epsilon \right)\left(a\sigma \left({\left(T_0^s\right)}^4-{\left(T_0^f\right)}^4\right)+h_{\mathrm{cov}}\left(T_0^s-T_0^f\right)\right)\end{array}$$E34

This equation is substituted into Eq. (21) to eliminate ${⟨T⟩}^s$ in favor of ${⟨T⟩}^f$. The resulting ordinary differential equation for ${⟨T⟩}^f$ runs as

$$\begin{array}{ll}{\frac{d^3⟨T⟩}{{\mathit{dx}}^3}}^f& =\frac{{\mathit{Gc}}_p}{{\epsilon }^{\ast }{k}_0{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^n+\epsilon k_{\mathit{disxx}}}{\frac{d^2⟨T⟩}{{\mathit{dx}}^2}}^f\\ & +h_v\frac{k_{\mathit{\text{stag}}}+\epsilon k_{\mathit{disxx}}+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3}{\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^n+\epsilon k_{\mathit{disxx}}\right)\left(\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3\right)}{\frac{d⟨T⟩}{\mathit{dx}}}^f\\ & -h_v\frac{{\mathit{Gc}}_p}{\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^n+\epsilon k_{\mathit{disxx}}\right)\left(\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3\right)}{⟨T⟩}^f\\ & +h_v\frac{I_{0}cos\xi -\left(1-\epsilon \right)\left(\mathit{a\sigma }\left({\left({⟨T⟩}_0^s\right)}^4-{\left({⟨T⟩}_0^f\right)}^4\right)+h_{\mathrm{cov}}\left({⟨T⟩}_0^s-{⟨T⟩}_0^f\right)\right)}{\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^n+\epsilon k_{\mathit{disxx}}\right)\left(\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3\right)}\end{array}$$E35

This ordinary differential equation, with the boundary conditions in Eqs. (29), (30) and (31) and the auxiliary asymptotic condition $d{⟨T⟩}_0^f/\mathit{dx}=d^2{⟨T⟩}_0^f/{\mathit{dx}}^2=0$, yields

$${⟨T⟩}^f=T_{\mathit{eq}}-\left(T_{\mathit{eq}}-{⟨T⟩}_0^f\right)e^{-\mathit{\gamma \lambda x}}$$E36

and

$${⟨T⟩}^s=T_{\mathit{eq}}-\left(T_{\mathit{eq}}-{⟨T⟩}_0^s\right)e^{-\mathit{\gamma \lambda x}}$$E37

where $\gamma$ is the positive real root, which can uniquely be determined from the following cubic equation:

$${\gamma }^3+\frac{{\mathit{Gc}}_p}{\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^n+\epsilon k_{\mathit{disxx}}\right)\lambda }{\gamma }^2-\gamma -\frac{{\mathit{Gc}}_p}{\left(k_{\mathit{\text{stag}}}+\epsilon k_{\mathit{disxx}}+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3\right)\lambda }=0$$E38

where

$$\lambda =\sqrt{\frac{\left(k_{\mathit{\text{stag}}}+\epsilon k_{\mathit{disxx}}+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3\right)h_v}{\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^n+\epsilon k_{\mathit{disxx}}\right)\left(\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3\right)}}$$E39

The solid phase temperature at the inlet ${⟨T⟩}_0^s$ and temperature at the thermal equilibrium, namely, $T_{\mathit{eq}}={⟨T⟩}_{\infty }^f={⟨T⟩}_{\infty }^s$, are given by

$${⟨T⟩}_0^s=T_{\mathit{eq}}+\left(T_{\mathit{eq}}-{⟨T⟩}_0^f\right)\frac{{\mathit{Gc}}_p+\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^n+\epsilon k_{\mathit{disxx}}\right)\mathit{\gamma \lambda }}{\left(\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3\right)\mathit{\gamma \lambda }}$$E40

and

$$T_{\mathit{eq}}={⟨T⟩}_0^f+\frac{I_{0}cos\xi -\left(1-\epsilon \right)\left(\mathit{a\sigma }\left({\left({⟨T⟩}_0^s\right)}^4-{\left({⟨T⟩}_0^f\right)}^4\right)+h_{\mathrm{cov}}\left({⟨T⟩}_0^s-{⟨T⟩}_0^f\right)\right)}{{\mathit{Gc}}_p}$$E41

respectively. Usually, the receiver length L is sufficiently long to reach the local thermal equilibrium. Thus, the average air and solid temperatures are evaluated from

$$\overline{{⟨T⟩}^f}=\frac{1-e^{-\mathit{\gamma \lambda L}}}{\mathit{\gamma \lambda L}}{⟨T⟩}_0^f+\left(1-\frac{1-e^{-\mathit{\gamma \lambda L}}}{\mathit{\gamma \lambda L}}\right)T_{\mathit{eq}}\cong \frac{1}{\mathit{\gamma \lambda L}}{⟨T⟩}_0^f+\left(1-\frac{1}{\mathit{\gamma \lambda L}}\right)T_{\mathit{eq}}$$E42

$$\overline{{⟨T⟩}^s}=\frac{1-e^{-\mathit{\gamma \lambda L}}}{\mathit{\gamma \lambda L}}{⟨T⟩}_0^s+\left(1-\frac{1-e^{-\mathit{\gamma \lambda L}}}{\mathit{\gamma \lambda L}}\right)T_{\mathit{eq}}\cong \frac{1}{\mathit{\gamma \lambda L}}{⟨T⟩}_0^s+\left(1-\frac{1}{\mathit{\gamma \lambda L}}\right)T_{\mathit{eq}}$$E43

As one of the most important performance parameters, the receiver efficiency is defined by

$$\eta =\frac{I_{0}cos\xi -\left(1-\epsilon \right)\left(\mathit{a\sigma }\left({\left({⟨T⟩}_0^s\right)}^4-{\left({⟨T⟩}_0^f\right)}^4\right)+h_{\mathrm{cov}}\left({⟨T⟩}_0^s-{⟨T⟩}_0^f\right)\right)}{I_{0}cos\xi }$$E44

Having established the temperature development, the momentum equation, that is, Eq. (20) along with the equation of state can easily be solved to find out the pressure distribution along the receiver as

$${⟨p⟩}^f=\sqrt{{\left({⟨p⟩}_0^f\right)}^2-2R\left(\frac{{\mu }_0}{K}{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^{n}G+{\mathit{bG}}^2\right)\left(\frac{1-e^{-\mathit{\gamma \lambda x}}}{\mathit{\gamma \lambda }}{⟨T⟩}_0^f+\left(x-\frac{1-e^{-\mathit{\gamma \lambda x}}}{\mathit{\gamma \lambda }}\right)T_{\mathit{eq}}\right)}$$E45

Under the low Mach approximation, the required pumping power per unit frontal area may be evaluated from

$$\begin{array}{c}PP=G\underset{0}{\overset{L}{{\displaystyle \int }}}-\frac{d{p}^f}{{\rho }_f{}^f}=G\underset{0}{\overset{L}{{\displaystyle \int }}}-\frac{1}{{\rho }_f{}^f}\frac{d{p}^f}{dx}dx\\ ={\left(\frac{R}{p_0^f}\right)}^{2}G\left(\frac{{\mu }_0}{K}{\left(\frac{\overline{T^f}}{T_0^f}\right)}^{n}G+b{G}^2\right)\underset{0}{\overset{L}{{\displaystyle \int }}}{\left(T^f\right)}^{2}dx\\ ={\left(\frac{R}{p_0^f}\right)}^{2}G\left(\frac{{\mu }_0}{K}{\left(\frac{\overline{T^f}}{T_0^f}\right)}^{n}G+b{G}^2\right)L(T_{eq}{}^2-\frac{2\left(1-e^{-\gamma \lambda L}\right)}{\gamma \lambda L}{T}_{eq}\left(T_{eq}-T_0^f\right)\\ +\frac{1-e^{-2\gamma \lambda L}}{2\gamma \lambda L}{\left(T_{eq}-T_0^f\right)}^2)\\ \cong \frac{G}{{\left({\rho }_0^f\right)}^2}\left(\frac{{\mu }_0}{K}{\left(\frac{\overline{T^f}}{T_0^f}\right)}^{n}G+b{G}^2\right)\frac{1}{2\gamma \lambda }\left(\left(2\gamma \lambda L-3\right){\left(\frac{T_{eq}}{T_0^f}\right)}^2+2\left(\frac{T_{eq}}{T_0^f}\right)+1\right)\end{array}$$E46

Note that the dynamic pressure change is sufficiently small as compared to the absolute pressure such that ${⟨{\rho }_f⟩}^f\propto 1/{⟨T⟩}^f$.

3.2 Analysis based on the P1 model

Since the Rosseland approximation used in the previous analysis ceases to be valid near boundaries, the validity of applying the Rosseland approximation near the inlet boundary of the receiver should be investigated. In order to examine the validity of the Rosseland approximation, the results based on the Rosseland approximation will be compared with the results obtained from solving the irradiation transport equation based on the P1 model. Since the silicon carbide ceramic foam is optically thick, the radiant energy emitted from other locations in the domain is quickly absorbed such that the radiative heat flux is given by

$$q_R=-\frac{1}{3\beta }\frac{\partial G}{\partial x_j}$$E47

where the diffuse integrated intensity $G_r$ satisfies the irradiation transport equation based on the P1 model as follows:

$$\frac{\partial }{\partial x_j}\left(\frac{1}{3\beta }\frac{\partial G_r}{\partial x_j}\right)+\kappa \left(4\sigma {\left({⟨T⟩}^s\right)}^4-G_r\right)=0$$E48

where $\kappa$ is the absorption coefficient.

Moreover, the effects of turbulence mixing on the heat transfer are also considered. Therefore, the energy equation for the air will be written as

$$\epsilon c_p\frac{\partial }{\partial x_j}{⟨{\rho }_f⟩}^f{⟨u_j⟩}^f{⟨T⟩}^f=\frac{\partial }{\partial x_j}\left({\epsilon }^{\ast }{⟨k_f⟩}^f\frac{\partial {⟨T⟩}^f}{\partial x_j}+\epsilon \left(k_{{\mathit{dis}}_{\mathit{jk}}}+\frac{c_{p_f}{⟨{\mu }_t⟩}^f}{{\sigma }_T}\right)\frac{\partial {⟨T⟩}^f}{\partial x_k}{\delta }_{\mathit{jk}}\right)-h_v\left({⟨T⟩}^f-{⟨T⟩}^s\right)$$E49

where turbulent Prandtl number ${\sigma }_T=0.9$ is assumed to be constant.

Under the low Mach number approximation, namely, we may reduce the macroscopic governing equations namely Eqs. (1), (49), (14) and (48) to a one-dimensional set of the equations as follows:

$$-\frac{d{⟨p⟩}^f}{\mathit{dx}}=\frac{{⟨\mu ⟩}^f}{K}\frac{G}{{⟨{\rho }_f⟩}^f}+b\frac{G^2}{{⟨{\rho }_f⟩}^f}=\frac{R}{{⟨p⟩}^f}\left(\frac{{⟨\mu ⟩}^f}{K}G+{\mathit{bG}}^2\right){⟨T⟩}^f$$E50

$$c_{p}G\frac{d{⟨T⟩}^f}{\mathit{dx}}=\frac{d}{\mathit{dx}}\left({\epsilon }^{\ast }{⟨k_f⟩}^f+\epsilon \left(k_{{\mathit{dis}}_{\mathit{xx}}}+\frac{c_{p_f}{⟨{\mu }_t⟩}^f}{{\sigma }_T}\right)\right)\frac{d{⟨T⟩}^f}{\mathit{dx}}-h_v\left({⟨T⟩}^f-{⟨T⟩}^s\right)$$E51

$$\frac{d}{\mathit{dx}}\left(\left(1-{\epsilon }^{\ast }\right)k_s\frac{d{⟨T⟩}^s}{\mathit{dx}}+\frac{1}{3\beta }\frac{{\mathit{dG}}_r}{\mathit{dx}}\right)-h_v\left({⟨T⟩}^s-{⟨T⟩}^f\right)=0$$E52

$$\frac{d}{\mathit{dx}}\left(\frac{1}{3\beta }\frac{{\mathit{dG}}_r}{\mathit{dx}}\right)+\kappa \left(4\sigma {\left({⟨T⟩}^s\right)}^4-G_r\right)=0$$E53

The turbulence kinetic energy is dropped from the momentum equation since it stays nearly constant within the receiver.

Nakayama and Kuwahara [22] established the macroscopic two-equation turbulence model, which does not require any detailed morphological information for the structure. The model, for given permeability and Forchheimer coefficient, can be used for analyzing most complex turbulent flow situations in homogeneous porous media. For the case of fully developed turbulent flow in an isotropic porous structure, the eddy viscosity is given by

$${⟨{\mu }_t⟩}^f=2\mathit{GbK}$$E54

Note that

$$k_{{\mathit{dis}}_{\mathit{xx}}}/\left(c_{p_f}{⟨{\mu }_t⟩}^f/{\sigma }_T\right)=0.6{\sigma }_T/\mathit{\epsilon b}\sqrt{K}>1$$

such that the dispersion thermal conductivity usually overwhelms the eddy thermal conductivity.

For absorption coefficient $\kappa$, the measurements made on cordierite ceramic foams by Kamiuto et al. [22] give the following correlation:

$$\kappa =4a\left(1-\epsilon \right)/d_m$$E55

The boundary conditions of ${⟨T⟩}^f$ and ${⟨p⟩}^f$ are the same as Eqs. (29) and (30). The other boundary conditions are given as follows:

$$q_{R_x}=-\frac{1}{3\beta }\frac{{\mathit{dG}}_r}{\mathit{dx}}=-\frac{G_r}{2}$$E56

and

$$\begin{array}{c}-\left({\epsilon }^{\ast }{k}_f{}^f+\epsilon \left(k_{disxx}+\frac{c_{p_f}{\mu }_t{}^f}{{\sigma }_T}\right)\right){\frac{d{T}^f}{dx}|}_{x=0}-{\left(1-{\epsilon }^{\ast }\right)k_s\frac{d{T}^s}{dx}|}_{x=0}-\frac{{G_r|}_{x=0}}{2}\\ =\left(1-\left(1-a\right)\left(1-\epsilon \right)\right)I_0\cos \xi \\ -\left(1-\epsilon \right)\left(a\sigma \left({\left(T_0^s\right)}^4-{\left(T_0^f\right)}^4\right)+h_{conv}\left(T_0^s-T_0^f\right)\right)\end{array}$$E57

Furthermore, the streamwise gradients of the dependent variables ${⟨T⟩}^f$, ${⟨T⟩}^s$ and $G_r$ are set to zero sufficiently far downstream at x = L.

The two energy equations, namely, Eqs. (51) and (52) may be added together and integrated using the boundary conditions in Eqs. (29) and (57) to give

$$\begin{array}{l}{c}_{p}G\left({⟨T⟩}^f-{⟨T⟩}_0^f\right)\\ =\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^n+\epsilon \left(k_{{\mathit{dis}}_{\mathit{xx}}}+\frac{c_{p_f}{⟨{\mu }_t⟩}^f}{{\sigma }_T}\right)\right)\frac{d{⟨T⟩}^f}{\mathit{dx}}+\left(1-{\epsilon }^{\ast }\right)k_s\frac{\partial {⟨T⟩}^s}{\partial x_j}+\frac{1}{3\beta }\frac{{\mathit{dG}}_r}{\mathit{dx}}\\ +\left(1-\left(1-a\right)\left(1-\epsilon \right)\right)I_{0}cos\xi -\left(1-\epsilon \right)\left(\mathit{a\sigma }\left({\left({⟨T⟩}_0^s\right)}^4-{\left({⟨T⟩}_0^f\right)}^4\right)+h_{\mathit{\text{conv}}}\left({⟨T⟩}_0^s-{⟨T⟩}_0^f\right)\right)\end{array}$$E58

Eqs. (52) and (53) are combined to give

$$G_r=4\sigma {\left({⟨T⟩}^s\right)}^4+\frac{h_v}{\kappa }\left({⟨T⟩}^s-{⟨T⟩}^f\right)-\frac{\left(1-{\epsilon }^{\ast }\right)k_s}{\kappa }\frac{d^2{⟨T⟩}^s}{{\mathit{dx}}^2}$$E59

This equation, Eq. (59), and Eq. (51) are substituted into Eq. (58) to eliminate $G_r$ and ${⟨T⟩}^s$ in favor of ${⟨T⟩}^f$. The resulting ordinary differential equation for ${⟨T⟩}^f$ runs as

$$\begin{array}{c}{\frac{d^{5}T}{d{x}^5}}^f=\frac{G{c}_p}{{\epsilon }^{\ast }{k}_0{\left(\frac{\overline{T^f}}{T_0^f}\right)}^n+\epsilon \left(k_{di{s}_{xx}}+\frac{c_{p_f}{\mu }_t{}^f}{{\sigma }_T}\right)}{\frac{d^{4}T}{d{x}^4}}^f\\ +\left(3\beta \kappa \frac{\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left(\overline{T^s}\right)}^3+\frac{h_v}{3\beta \kappa }}{\left(1-{\epsilon }^{\ast }\right)k_s}+\frac{h_v}{{\epsilon }^{\ast }{k}_0{\left(\frac{\overline{T^f}}{T_0^f}\right)}^n+\epsilon \left(k_{di{s}_{xx}}+\frac{c_{p_f}{\mu }_t{}^f}{{\sigma }_T}\right)}\right){\frac{d^{3}T}{d{x}^3}}^f\\ -\left(\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left(\overline{T^s}\right)}^3+\frac{h_v}{3\beta \kappa }\right)\frac{3\beta \kappa G{c}_p}{\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{T^f}}{T_0^f}\right)}^n+\epsilon \left(k_{di{s}_{xx}}+\frac{c_{p_f}{\mu }_t{}^f}{{\sigma }_T}\right)\right)\left(1-{\epsilon }^{\ast }\right)k_s}\frac{d^2{T}^f}{d{x}^2}\\ -3\beta \kappa h_v\frac{k_{stag}+\epsilon \left(k_{di{s}_{jk}}+\frac{c_{p_f}{\mu }_t{}^f}{{\sigma }_T}\right)+\frac{16\sigma }{3\beta }{\left(\overline{T^s}\right)}^3}{\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{T^f}}{T_0^f}\right)}^n+\epsilon \left(k_{di{s}_{xx}}+\frac{c_{p_f}{\mu }_t{}^f}{{\sigma }_T}\right)\right)\left(1-{\epsilon }^{\ast }\right)k_s}\frac{d{T}^f}{dx}\\ +\frac{3\beta \kappa h_{v}G{c}_p}{\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{T^f}}{T_0^f}\right)}^n+\epsilon \left(k_{di{s}_{xx}}+\frac{c_{p_f}{\mu }_t{}^f}{{\sigma }_T}\right)\right)\left(1-{\epsilon }^{\ast }\right)k_s}\left(T^f-T_0^f\right)\\ -3\beta \kappa h_v\frac{\left(1-\left(1-a\right)\left(1-\epsilon \right)\right)I_0\cos \xi -\left(1-\epsilon \right)\left(a\sigma \left({\left(T_0^s\right)}^4-{\left(T_0^f\right)}^4\right)+h_{conv}\left(T_0^s-T_0^f\right)\right)}{\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{T^f}}{T_0^f}\right)}^n+\epsilon \left(k_{di{s}_{xx}}+\frac{c_{p_f}{\mu }_t{}^f}{{\sigma }_T}\right)\right)\left(1-{\epsilon }^{\ast }\right)k_s}\end{array}$$E60

This ordinary differential equation, with the boundary conditions in Eqs. (29), (56) and (57) and the zero derivative conditions far downstream ($x\to \infty$: Note L is sufficiently large), yields Eqs. (36) and (37). Note that $\gamma$ is the positive real root which can be determined from the following characteristic equation:

$$\begin{array}{l}{\gamma }^5+\frac{{\mathit{Gc}}_p}{\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^n+\epsilon \left(k_{{\mathit{dis}}_{\mathit{xx}}}+\frac{c_{p_f}{⟨{\mu }_t⟩}^f}{{\sigma }_T}\right)\right)\lambda }{\gamma }^4\\ -\frac{\frac{3\mathit{\beta \kappa }}{h_v}\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^n+\epsilon \left(k_{{\mathit{dis}}_{\mathit{xx}}}+\frac{c_{p_f}{⟨{\mu }_t⟩}^f}{{\sigma }_T}\right)\right)\left(\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3+\frac{h_v}{3\mathit{\beta \kappa }}\right)+\left(1-{\epsilon }^{\ast }\right)k_s}{k_{\mathit{\text{stag}}}+\epsilon \left(k_{{\mathit{dis}}_{\mathit{xx}}}+\frac{c_{p_f}{⟨{\mu }_t⟩}^f}{{\sigma }_T}\right)+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3}{\gamma }^3\\ -\frac{\left(1-{\epsilon }^{\ast }\right)k_s+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3+\frac{h_v}{3\mathit{\beta \kappa }}}{k_{\mathit{\text{stag}}}+\epsilon \left(k_{{\mathit{dis}}_{\mathit{xx}}}+\frac{c_{p_f}{⟨{\mu }_t⟩}^f}{{\sigma }_T}\right)+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3}\frac{3\mathit{\beta \kappa }{\mathit{Gc}}_p}{h_v\lambda }{\gamma }^2+\frac{3\mathit{\beta \kappa }}{{\lambda }^2}\gamma \\ +\frac{3\mathit{\beta \kappa }{\mathit{Gc}}_p}{\left(k_{\mathit{\text{stag}}}+\epsilon \left(k_{{\mathit{dis}}_{\mathit{xx}}}+\frac{c_{p_f}{⟨{\mu }_t⟩}^f}{{\sigma }_T}\right)+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3\right){\lambda }^3}=0\end{array}$$E61

where

$$\lambda =\sqrt{\frac{\left(k_{\mathit{\text{stag}}}+\epsilon \left(k_{{\mathit{dis}}_{\mathit{xx}}}+\frac{c_{p_f}{⟨{\mu }_t⟩}^f}{{\sigma }_T}\right)+\frac{16\sigma }{3\beta }{\left(\overline{{⟨T⟩}^s}\right)}^3\right)h_v}{\left({\epsilon }^{\ast }{k}_0{\left(\frac{\overline{{⟨T⟩}^f}}{{⟨T⟩}_0^f}\right)}^n+\epsilon \left(k_{{\mathit{dis}}_{\mathit{xx}}}+\frac{c_{p_f}{⟨{\mu }_t⟩}^f}{{\sigma }_T}\right)\right)\left(1-{\epsilon }^{\ast }\right)k_s}}$$E62

The solid phase temperature at the inlet ${⟨T⟩}_0^s$ and temperature at the thermal equilibrium, namely,$T_{\mathit{eq}}={⟨T⟩}_{\infty }^f={⟨T⟩}_{\infty }^s$, are determined from the following implicit equations:

$$\begin{array}{ll}{\mathit{Gc}}_p\left(T_{\mathit{eq}}-{⟨T⟩}_0^f\right)=& -\left({\epsilon }^{\ast }{k}_0+\epsilon \left(k_{{\mathit{dis}}_{\mathit{xx}}}+\frac{c_{p_f}{⟨{\mu }_t⟩}^f}{{\sigma }_T}\right)\right)\mathit{\gamma \lambda }\left(T_{\mathit{eq}}-{⟨T⟩}_0^f\right)\\ & -\left(1-{\epsilon }^{\ast }\right)k_s\mathit{\gamma \lambda }\left(T_{\mathit{eq}}-{⟨T⟩}_0^s\right)\\ & -\frac{1}{2}\left(4\sigma {\left({⟨T⟩}_0^s\right)}^4+\frac{h_v}{\kappa }\left({⟨T⟩}_0^s-{⟨T⟩}_0^f\right)+\frac{\left(1-{\epsilon }^{\ast }\right)k_s}{\kappa }{\left(\mathit{\gamma \lambda }\right)}^2\left(T_{\mathit{eq}}-{⟨T⟩}_0^s\right)\right)\end{array}$$E63

$$T_{\mathit{eq}}={⟨T⟩}_0^f+\frac{\left(1-\left(1-a\right)\left(1-\epsilon \right)\right)I_{0}cos\xi -\left(1-\epsilon \right)\left(\mathit{a\sigma }\left({\left({⟨T⟩}_0^s\right)}^4-{\left({⟨T⟩}_0^f\right)}^4\right)+h_{\mathit{\text{conv}}}\left({⟨T⟩}_0^s-{⟨T⟩}_0^f\right)\right)}{{\mathit{Gc}}_p}$$E64

where the boundary condition in Eq. (56) is utilized. Usually, the receiver length L is sufficiently long to reach the local thermal equilibrium. Thus, the average air and solid temperatures are evaluated from Eqs. (42) and (43).

4. Validations of the Rosseland approximation

Smirnova et al. [23] numerically studied the compressible fluid flow and heat transfer within the solar receiver with silicon carbide monolithic honeycombs. In their paper, the following input data were collected to obtain the analytic solutions based on the present local thermal non-equilibrium model:

$$\begin{array}{l}{⟨{\rho }_f⟩}_0^f=1.2\left[\mathrm{kg}/{\mathrm{m}}^3\right],{⟨T⟩}_0^f=300\left[\mathrm{K}\right]({⟨p⟩}_0^f={10}^5\left[\mathrm{Pa}\right]),c_p=1000\left[\mathrm{J}/\mathrm{kgK}\right],\\ G=1.2\left[\mathrm{kg}/{\mathrm{m}}^2\mathrm{s}\right],L=0.05\left[m\right],I_0={10}^6\left[\mathrm{W}/{\mathrm{m}}^2\right],\xi =0,h_{\mathit{\text{conv}}}=0\left[\mathrm{W}/{\mathrm{m}}^2\mathrm{K}\right],\\ k_s=150\left[\mathrm{W}/\mathrm{mK}\right]k_{\mathit{dis}}=0\left[\mathrm{W}/\mathrm{mK}\right],h_v=8.8\times {10}^4\left[\mathrm{W}/{\mathrm{m}}^3\mathrm{K}\right],\epsilon =0.5,\beta =50\left[1/\mathrm{m}\right].\end{array}$$

However, it should be noticed that the porosity of the silicon carbide monolithic honeycombs is not available in Smirnova et al. [23], its value was estimated to be $\epsilon =0.5$ from the figure provided by Agrafiotis et al. [24]. The mean extinction coefficient $\beta$ for silicon carbide monolithic honeycombs is not available in their paper. Finally, the value was estimated to be 50[1/m] by correlating the present results against theirs. It should also be noted that the convective heat transfer coefficient was set to zero since radiation predominates over convection in the receiver front.

As for possible instabilities, the fifth-order characteristic Eq. (61) based on the P1 model should be examined carefully. Figure 2 shows the residual of the fifth-order characteristic equation $f_R\left(\gamma \right)$. The figure clearly shows that the fifth-order characteristic Eq. (61) under a possible range of the silicon carbide parameters yields two positive roots ${\gamma }_h$ and ${\gamma }_l$, which are fairly close to each other. The corresponding temperature variations of both phases however depend strongly on its value, which results in a non-unique value of equilibrium temperature. Since flow instability is inferred by an unexpected nature of the quadratic pressure difference with respect to equilibrium temperature, the existence of two positive roots may be responsible for possible hydrodynamic and thermal instabilities reported previously. A further investigation based on an unsteady procedure is definitely needed to explore possible causes of these instabilities, closely related to the radiative heat transfer mode.

Figure 2.

Residual of the fifth-order characteristic equation.

The third-order characteristic Eq. (38) based on the Rosseland approximation, on the other hand, yields only one positive root ${\gamma }_1$. The corresponding fluid and solid temperature variations generated under the Rosseland approximation are compared with those based on the P1 model with the larger root ${\gamma }_h$. Figure 3 shows that both sets of the temperature developments agree fairly well with each other. Thus, the Rosseland approximation for this case, despite its failure near the inlet boundary, is fairly accurate and may well be used for quick estimations and further analysis.

Figure 3.

Comparison of the temperature developments with the Rosseland approximation and P1 model.

In Figure 4, the present analytic solutions are compared against the large-scale FEM numerical calculations based on COMSOL, reported by Smirnova et al. [23]. It should be mentioned that the direct numerical integrations of Eqs. (20)–(22) were also carried out using the finite volume method code, SAINTS [12]. As the convergence criteria, the residuals of all equations are less than 10⁻⁵. It can be clearly seen that the air temperature increases as receiving heat from the monolithic receiver. Eventually, these two phases reach local thermal equilibrium near the exit. Both sets of solutions agree very well with each other, indicating the validity of the present local thermal non-equilibrium model.

Figure 4.

Axial developments of the fluid and solid phase temperatures: comparison of the present analysis and FEM analysis.

5. Applications to silicon carbide ceramic foam volumetric receiver

In order to overcome the problems associated with thermal spots and flow instabilities, we would like to study fluid flow and heat transfer characteristics in silicon carbide ceramic foams based on the analytical expressions of pressure and temperature fields within a solar volumetric receiver. The performance of the receiver may be assessed in terms of the receiver efficiency $\eta$ under equal pumping power $\mathit{PP}$. Thus, the effects of the pore diameter $d_m$ on the receiver efficiency $\eta$ are presented in Figure 5, since $d_m$ is a crucial geometry parameter affecting hydrodynamic and thermal characteristics of foam shown in Eqs. (23), (24) and (27). The pore diameter $d_m$ is varied whereas the other parameters are fixed as follows:

Figure 5.

Effects of the pore diameter on the receiver efficiency.

$$\begin{array}{l}{⟨{\rho }_f⟩}_0^f=1.2\left[\mathrm{kg}/{\mathrm{m}}^3\right],{⟨T⟩}_0^f=300\left[\mathrm{K}\right]({⟨p⟩}_0^f={10}^5\left[\mathrm{Pa}\right]),c_p=1000\left[\mathrm{J}/\mathrm{kgK}\right],\\ L=0.03\left[m\right],I_0={10}^6\left[\mathrm{W}/{\mathrm{m}}^2\right],\xi =0,h_{\mathit{\text{conv}}}=0\left[\mathrm{W}/{\mathrm{m}}^2\mathrm{K}\right],k_s=150\left[\mathrm{W}/\mathrm{mK}\right],\epsilon =0.9.\end{array}$$

All other parameters are evaluated using Eqs. (17), (18) and from Eq. (23) to (28).

As shown in Figure 5, it is interesting to note that $\eta$ suddenly increases at some critical value of $d_m$ for a given value of PP, which means that the pore diameter $d_m$ must be larger than this critical value to achieve high $\eta$. This finding is useful to design a volumetric receiver, and can be interpreted in what follows.

As indicated in Eq. (46), it can be easily deduced that $G\propto \sqrt{\mathit{PP}}$ for low PP and $G\propto \sqrt[3]{\mathit{PP}}$ for high PP, which results in that the amount of heat carried by the air, $G\left(T_{\mathit{eq}}-{⟨T⟩}_0^f\right)\propto \sqrt{\mathit{PP}}$, increases drastically on increasing the pumping power $\mathit{PP}$ from zero. Nevertheless, its rate of increase diminishes for the higher PP range, in which $G\left(T_{\mathit{eq}}-{⟨T⟩}_0^f\right)\propto \sqrt[3]{\mathit{PP}}$. Moreover, it can also be concluded that the sharp rise in the receiver efficiency occurs around the transition from the Darcy to Forchheimer regime, namely,

$$\frac{{\mu }_0}{K}{\left(\frac{T_{\mathit{eq}}}{{⟨T⟩}_0^f}\right)}^n{G}_{\mathit{tr}}\cong {{\mathit{bG}}_{\mathit{tr}}}^2$$E65

or

$$G_{\mathit{tr}}\cong {\left(\frac{{\mu }_0}{\mathit{bK}}{\left(\frac{I_{0}cos\xi }{c_p{⟨T⟩}_0^f}\right)}^n\right)}^{\frac{1}{1+n}}$$E66

since

$$\frac{T_{\mathit{eq}}}{{⟨T⟩}_0^f}\cong \frac{c_p{G}_{\mathit{tr}}{⟨T⟩}_0^f+I_{0}cos\xi }{c_p{G}_{\mathit{tr}}{⟨T⟩}_0^f}\cong \frac{I_{0}cos\xi }{c_p{G}_{\mathit{tr}}{⟨T⟩}_0^f}$$E67

Thus, Eq. (46) may be written for the case in which the sharp rise in $\eta$ takes place as follows:

$$\mathit{PP}=\frac{G_{\mathit{tr}}}{{\left({⟨\rho ⟩}_0^f\right)}^2}\left(2{{\mathit{bG}}_{\mathit{tr}}}^2\right){\left(\frac{T_{\mathit{eq}}}{{⟨T⟩}_0^f}\right)}^{2}L\cong \frac{2{{\mathit{bLG}}_{\mathit{tr}}}^3}{{\left({⟨\rho ⟩}_0^f\right)}^2}{\left(\frac{I_{0}cos\xi }{c_p{G}_{\mathit{tr}}{⟨T⟩}_0^f}\right)}^2=\frac{2\mathit{bL}}{{\left({⟨\rho ⟩}_0^f\right)}^2}{\left(\frac{I_{0}cos\xi }{c_p{⟨T⟩}_0^f}\right)}^2{\left(\frac{{\mu }_0}{\mathit{bK}}{\left(\frac{I_{0}cos\xi }{c_p{⟨T⟩}_0^f}\right)}^n\right)}^{\frac{1}{1+n}}$$E68

which, for given PP, gives the minimum value of the pore diameter $d_{\mathit{mtr}}$:

$$\frac{d_{\mathit{mtr}}}{L}=f\left(\epsilon \right){\left(\frac{2}{{\left({⟨\rho ⟩}_0^f\right)}^2\mathit{PP}}{\left(\frac{I_{0}cos\xi }{c_p{⟨T⟩}_0^f}\right)}^3\right)}^{\frac{1+n}{2+n}}{\left(\frac{{\mu }_0}{\left(\frac{I_{0}cos\xi }{c_p{⟨T⟩}_0^f}\right)L}\right)}^{\frac{1}{2+n}}$$E69

$$f\left(\epsilon \right)={\left(\frac{{\left({\mathit{bd}}_m\right)}^n}{K/{d_m}^2}\right)}^{\frac{1}{2+n}}={\left(\frac{{\left(12\left(1-\epsilon \right)\right)}^n}{0.00073{\left(1-\epsilon \right)}^{-0.224}{\left(\frac{1.18}{1-e^{-\left(1-\epsilon \right)/0.04}}\sqrt{\frac{1-\epsilon }{3\pi }}\right)}^{-1.11}}\right)}^{\frac{1}{2+n}}$$E70

For PP = 300, 500 and 1000 W/m² studied here, Eq. (69) gives $d_{\mathit{mtr}}$= 0.0022, 0.0016 and 0.0010 m, respectively. It is consistent with what is observed in Figure 5, since an increase in $d_m$ (i.e., decrease in $\beta$) from $d_{\mathit{mtr}}$ makes further penetration of the solar radiation possible. This works to keep the solid temperature at the inlet comparatively low such that heat loss to the ambient by radiation is suppressed. As a result, high receiver efficiency can be achieved. However, the increase in $d_m$ on the other hand results in decreasing the volumetric heat transfer coefficient, as can be seen from Eq. (27). Too large $d_m$ deteriorates interstitial heat transfer from the solid to air. Thus, as can be seen from the figure, the optimal size of $d_m$ exits under the equal pumping power constraint.

In order to achieve local thermal equilibrium for the two phases within the receiver, the length of the receiver is assumed to be sufficiently long in the present study. In view of minimizing the required pumping power, however, it is noticeable that shorter length is better, as clearly seen from Eq. (46). Hence, a minimum length required to approach local thermal equilibrium may be chosen to design a receiver, which would guarantee both maximum receiver efficiency and minimum pumping power. Therefore, we may roughly set the optimal receiver length as

$$L=\frac{3}{\mathit{\gamma \lambda }}$$E71

such that

$$\frac{{\left.{⟨T⟩}^f\right|}_{x=L}-T_{\mathit{eq}}}{{⟨T⟩}_0^f-T_{\mathit{eq}}}=\frac{{\left.{⟨T⟩}^s\right|}_{x=L}-T_{\mathit{eq}}}{{⟨T⟩}_0^s-T_{\mathit{eq}}}=e^{-3}\cong 5\%$$E72

Eq. (71) together with Eq. (69) provides useful information for designing a volumetric solar receiver of silicon carbide ceramic foam.

6. Conclusions

For the first time, the complete set of analytical solutions, which fully considers the combined effects of turbulence, tortuosity, thermal dispersion, compressibility on the convective, conductive and radiative heat transfer within a ceramic foam receiver, is presented based on the two-energy equation model of porous media. Both the Rosseland approximation and the P1 model are applied to account for the radiative heat transfer through the solar receiver, while the low Mach approximation is exploited to investigate the compressible flow through the receiver. Based on the P1 model, two positive roots were found from the characteristic equations of the fifth-order differential equation, indicating possible occurrence of hydrodynamic and thermal instabilities. However, it has been found that the Rosseland approximation for this case, despite its failure near the inlet boundary, is fairly accurate and may well be used for quick estimations and further analysis. Due to their advantages, such as high thermal conductivity and fluid mixing, silicon carbide ceramic foams are considered as a possible candidate for the receiver, which can overcome the problems associated with thermal spots and flow instabilities. The results show that the pore diameter must be larger than its critical value to achieve high receiver efficiency. As a result, there exists an optimal pore diameter for achieving the maximum receiver efficiency under the equal pumping power. The optimal pore diameter yielding the maximum receiver efficiency may be found around the critical value given by Eq. (71). A simple relation is derived for determining the length of the volumetric solar receivers of silicon carbide ceramic foam.

Conflict of interest

The authors declare no conflict of interest.

Nomenclature

A	surface area (m2)
Aint	interfacial surface area between the fluid and solid (m2)
b	inertial coefficient (1/m)
c	specific heat (J/kg K)
cp	specific heat at constant pressure (J/kg K)
dm	pore diameter (m)
G	mass flux (kg/m2 s)
h	specific enthalpy (J/kg)
hv	volumetric heat transfer coefficient (W/m3K)
I0	intensity of radiation (W/m2)
k	thermal conductivity (W/m K)
K	permeability (m2)
L	receiver length (m)
nj	normal unit vector from the fluid side to solid side (−)
PP	pumping power per unit frontal area (W/m2)
Pr	Prandtl number (−)
q	heat flux (W/m2)
R	gas constant (J/kg K)
T	temperature (K)
ui	velocity vector (m/s)
V	representative elementary volume (m3)
xi	Cartesian coordinates (m)
x	axial coordinate (m)
β	mean extinction coefficient (1/m)
γ	dimensionless parameter (−)
ε	porosity (−)
ε∗	effective porosity (−)
ξ	incidence angle (rad)
η	receiver efficiency (−)
λ	characteristic coefficient (1/m)
μ	viscosity (Pa s)
ν	kinematic viscosity (m2/s)
ρ	density (kg/m2)
σ	Stephan-Boltzmann constant (W/m2K4)
κ	absorption coefficient (1/m)
τij	stress tensor (Pa)
Special symbols
φ˜	deviation from intrinsic average
ϕ	Darcian average
ϕf,s	intrinsic average
Subscripts and superscripts
dis	dispersion
eq	equilibrium
f	fluid
s	solid
stag	stagnation
0	reference