Transient radiative heat transfer within a suspension of coal particles undergoing steam gasification

Abstract

Transient radiative heat transfer in chemical reacting media is examined for a non-isothermal, non-gray, absorbing, emitting, and Mie-scattering suspension of coal particles, whose radiative properties vary with time as the particles undergo shrinking by endothermic gasification. A numerical model that incorporates parallel filtered collision-based Monte Carlo ray tracing, finite volume method, and explicit Euler time integration scheme is formulated for solving the unsteady energy equation that couples the radiative heat flux with the chemical kinetics. Variation of radiative properties, attenuation characteristics, temperature profiles, and extent of the chemical reaction are reported as a function of time. It is found that radiation in the visible and near IR spectrum incident on a cloud of coal particles greater than 2.5μm is more likely to be forward scattered than absorbed, but the opposite is true as the particles shrink below 1.3μm. The medium becomes optically thinner as the particles shrink and this effect is more pronounced for smaller initial coal particles because these offer higher volume fraction to particle diameter ratio and, consequently, attain higher temperatures, reaction rates, and shrinking rates.

Appendix

1.1 Accuracy of the combined Monte Carlo method with time integration schemes

The accuracy of MC combined with the different time integration schemes is determined by solving for the temperature distribution within a plane layer of a gray-isotropic participating medium, and comparing the results with those obtained using a semi-analytical method. A 1D plane layer of a non-isothermal absorbing, emitting, and isotropically scattering gray medium is considered. The medium is contained within black plane boundaries at constant temperature T_b. Its initial temperature is T₀. The equations developed are valid for either T₀ > T_b or T₀ < T_b, i.e. the medium undergoing either cooling or heating, respectively. Medium properties are listed in Table 3 and are assumed arbitrary. Neglecting convection, conduction, pressure work, and internal heat generation, the general energy equation for one-component medium is given by:

$$\rho c_{\text{v}} \frac{{\partial T}}{{\partial t}} = - \nabla \cdot \vec q_{\text{r}} ^{\prime \prime } $$

(24)

where $\vec q_{\text{r}}^{\prime\prime}$ is the radiative heat flux across the medium. Two approaches for finding the divergence of $\vec q_{\text{r}}^{\prime\prime}$ are presented: (1) the semi-analytical method, and (2) the MC method. Finite volume technique is employed in both approaches for discretization of Eq. 24. In particular, for a 1D geometry divided into a large number of sub-layers, each of same thickness and at uniform temperature:

$$T_j^{n + 1} = T_j^n + \int\limits_{\Delta t} {\frac{1}{{\rho _j c_{{\text{v}}j} \Delta x_j }}\int\limits_{\Delta x} {\left( { - \frac{{\partial q_{\text{r}} ^{\prime \prime } }}{{\partial x}}} \right)} } _j {\text{d}}x{\text{d}}t$$

(25)

where j denotes a sub-layer and Δt= t ^{n +1}−tⁿ is the time step interval.

Table 3 Baseline parameters used for the gray-isotropic medium

1.1.1 The semi-analytical method

The radiative flux divergence can be analytically derived from the equation of radiative transfer [25]:

$$\frac{{{\text{d}}q_{\text{r}} ^{\prime \prime } }}{{{\text{d}}x}} = 4\frac{k}{\omega }\left[ {\sigma T^4 (x,t) - \pi S(x,t)} \right]$$

(26)

The source function can be expressed in terms of the exponential integral functions,

$$S(x,t) = (1 - \omega )\frac{{\sigma T^4 (x,t)}}{\pi } + \frac{\omega }{2}\left\{ {\frac{{\sigma T_{\text{b}}^{\text{4}} }}{\pi }\left[ {E_2 (x) + E_2 (L - x)} \right] + \int\limits_0^L {S(x^ {*} } ,t)E_1 \left( {|x^ {*} - x|} \right){\text{d}}x^ {*} } \right\}$$

(27)

Inserting (26) into (25) and further assuming that the flux divergence and field variables are constant over each sub-layer, an expression for the temperature in each sub-layer can be derived:

$$T^{n + 1} = T^n + \frac{4}{{\rho c_{\text{v}} }}\frac{k}{\omega }\int\limits_{t_n }^{t_{n + 1} } {\left[ {\pi S\left( {x,t} \right) - \sigma T^4 \left( {x,t} \right)} \right]} {\text{d}}t$$

(28)

where the shorthand notation Tⁿ=T(t) ⁿ) is introduced. Note that for simplicity, spatial discretization in sub-layers is omitted from the notation of Eq. 28 and succeeding equations. The temperature distribution is obtained using Eq. 28, where the integral is computed using the different time integration schemes presented below.

1.1.2 The MC method

The integral of the divergence of the radiative flux over a sub-layer, shown in the right-hand side of Eq. 25, represents the net radiative power absorbed by that sub-layer, i.e. the difference between absorbed and emitted power:

$$q_{\text{r}} = q_{\text{a}} - q_{\text{e}} = \int\limits_{\Delta x} {\left( { - \frac{{\partial q_{\text{r}} ^{\prime \prime } }}{{\partial x}}} \right){\text{d}}x} $$

(29)

Inserting Eq. 29 into Eq. 25 yields,

$$T^{n + 1} = T^n + \frac{1}{{\rho c_{\text{v}} \Delta x}}\int\limits_{t^n }^{t^{n + 1} } {\left( {q_{\text{a}} - q_{\text{e}} } \right){\text{d}}t} $$

(30)

where

$$q_{\text{e}} = 4kV\sigma T^4 $$

(31)

and q _a is found by MC as the number of rays absorbed in a sub-layer times the power carried by a single ray. At the boundaries, the emitted heat fluxes are:

$$q_{{\text{e,b}}_{\text{1}} } ^{\prime \prime } = - q_{{\text{e,b}}_{\text{2}} } ^{\prime \prime } = \varepsilon _{\text{b}} \sigma T_{\text{b}}^{\text{4}} $$

(32)

Temperatures and, consequently, the total power emitted vary with time. In contrast, the total number of rays for each MC iteration is set constant with time. Thus, the power carried by each ray is calculated at each MC iteration as the ratio of the total power emitted to the total number of rays. The temperature distribution is obtained using Eq. 30, where the integral is computed using the different time integration schemes presented below.

1.1.3 Time integration schemes

Five schemes for time integration are considered [7, 5]:

$${\text{Explicit}}\;{\text{Euler}}:\;T^{n + 1} = T^n + f\left( {t^n ,T^n } \right)\Delta t$$

(33)

$${\text{Implicit}}\;{\text{Euler}}:\;T^{n + 1} = T^n + f\left( {t^{n + 1} ,T^{n + 1} } \right)\Delta t$$

(34)

$${\text{Crank}}\;{\text{Nicolson:}}\;T^{n + 1} = T^n + \left[ {f\left( {t^n ,T^n } \right) + f\left( {t^{n + 1} ,T^{n + 1} } \right)} \right]\frac{{\Delta t}}{2}$$

(35)

$${\text{2nd - order}}\;{\text{Runge - Kutta:}}\;T^{n + \frac{1}{2}} = T^n + f\left( {t^n ,T^n } \right)\frac{{\Delta t}}{2}$$

(36)

$$T^{n + 1} = T^n + f\left( {t^{n + \frac{1}{2}} ,T^{n + \frac{1}{2}} } \right)\Delta t$$

(37)

$$ {\text{4th - order}}\;{\text{Runge - Kutta:}}\;T^{n + \frac{1} {2},a} = T^n + f\left( {t^n ,T^n } \right)\frac{{\Delta t}} {2} $$

(38)

$$T^{n + \frac{1}{2},b} = T^n + f\left( {t^{n + \frac{1}{2}} ,T^{n + \frac{1}{2},a} } \right)\frac{{\Delta t}}{2}$$

(39)

$$T^{n + 1,a} = T^n + f\left( {t^{n + \frac{1}{2}} ,T^{n + \frac{1}{2},b} } \right)\Delta t$$

(40)

$$T^{n + 1} = T^n + \left[ {f\left( {t^n ,T^n } \right) + 2f\left( {t^{n + \frac{1}{2}} ,T^{n + \frac{1}{2},a} } \right) + 2f\left( {t^{n + \frac{1}{2}} ,T^{n + \frac{1}{2},b} } \right) + f\left( {t^{n + 1} ,T^{n + 1,a} } \right)} \right]\frac{{\Delta t}}{6}$$

(41)

where f denotes a time-dependent function. The explicit Euler and both Runge-Kutta schemes are self-starting. They only require an initial temperature distribution at t = 0. The implicit Euler and Crank-Nicolson schemes require, in addition, the temperature distribution at time level n+1, which is found iteratively. The MC method requires ray-tracing iterative runs at each time level, i.e. each time the integral on the right-hand side of (30) is evaluated. As long as such a procedure is followed, the generalization of the MC method to problems involving time-dependent radiation properties is straightforward and only requires the additional computation of these properties at each time level.

1.1.3.1 Parallelization and smoothing filters for MC

The described MC procedure requires significant computational resources when applied in conjunction with time integration schemes. Parallelization and smoothing filter concepts help alleviate this problem. Parallel methods are especially suitable for MC because they take advantage of the fact that tracing a generic ray constitutes a Markov chain’s event and, consequently, can be treated fully independent of tracing other generic rays. Different concepts of parallelization of MC are outlined in [4]. In this work we use parallelization by spatial region for the energy balance and by ray for ray tracing. Parallel MC computing concept has been previously applied for steady-state radiation problems involving iterative ray tracing [13], and it is here applied for the first time to transient radiation heat transfer. The proposed mechanism encompasses a master process and several slave processes, each responsible for tracing assigned rays within the whole domain. Material properties and grid configuration are obtained at the master process and forwarded to all slave processes, where the temperatures of the assigned sub-layers are computed by summing up q _a obtained from all processes. Subsequently, the time integration schemes are used to move forward in time. Critical to the implementation of this mechanism is the correct generation of pseudo random numbers and this study uses the algorithm developed by Press et al. [22]. The spatial noise level in the MC solution was reduced by employing a 33-point quadratic lower-pass Savitzky–Golay smoothing filter [22].

Fig. 10

Temperature profiles calculated by the “reference” semi-analytical method and by the filtered and unfiltered MC method using n_rays=10⁵ and the 4th-order Runge-Kutta time integration scheme

1.1.4 Results

The MC was performed using all of the aforementioned time schemes, for samples containing n_rays = 10⁴, 10⁵, and 10⁶ rays, for time steps Δt = 2, 4, and 8× 10⁻³ s with and without smoothing filter. In addition, the semi-analytical method was applied using the 4th-order Runge-Kutta time integration scheme for Δt = 10⁻³ s, yielding the most accurate solution, referred to as the “reference” solution. Common baseline parameters are listed in Table 3. Since for this specific example material properties are assumed constant with time, results are presented in a non-dimensional form:

$$ \tilde t = \frac{{4\sigma T^{3}_{0} }} {{f{\text{v}}\rho c_{{\text{v}}} L}}t $$

(42)

$$\tilde T = \frac{T}{{T_0 }}$$

(43)

Figure 10 shows the temperature profiles for the reference solution and for the filtered/unfiltered MC results using n_rays=10⁵ with the 4th-order Runge-Kutta integration scheme. The unfiltered results are indicated by the dots scattered around the reference solution. Scattering increases as steady state is approached. The temperature profile for the filtered results coincides well with the reference curve. The accuracy of the MC method is determined by calculating the error in the temperature profile based on the lumped squared relative differences between the actual and the reference solution, according to:

$$\delta ^2 (t) = \frac{1}{L}\int\limits_L {\left[ {1 - \frac{{T\left( {x,t} \right)}}{{T_{{\text{ref}}} \left( {x,t} \right)}}} \right]^2 {\text{d}}x} $$

(44)

where T_ref denotes the temperature distribution obtained for the reference (semi-analytical) solution and T the one obtained for the MC solution. This error is caused by the approximation in the time and space discretization due to a finite time and space intervals, and by the statistical approximation due to a finite sample of rays. Computations with refined space grid were carried out until the error due to space discretization could be neglected for a space grid resolution L/Δx = 250. Figure 11 shows the error in the filtered MC solution, calculated by Eq. 44, as a function of time for the cases listed in Table 4.

Fig. 11

Error of the MC method (cases are listed in Table 4)

Table 4 Parameters for the MC method

Obviously, the accuracy is improved with shorter time step intervals, larger samples of rays, and with the use of a filter. The explicit Euler scheme exhibited the best stability with longer time steps and minimum computational time. As expected, the 4th-order Runge-Kutta scheme gives the highest accuracy, with an error of less than 1% for n_rays ≥ 10⁵ and $\Delta\tilde t \leqslant 1.12.$ Increasing the sample of rays by an order of magnitude decreases the error roughly by half an order of magnitude. The introduction of filters further improves the accuracy such that the error for n_rays = 10⁵ and without filter is comparable to the one for n_rays = 10⁴ and with filter. For the implicit Euler or the Crank-Nicolson schemes, the MC method does not converge. For the 2nd-order or the 4th-order Runge-Kutta schemes with $\Delta \tilde t \geqslant 1.12,$ n_rays≤10⁴, and without low-pass filter, the MC method does not converge either.