**Advances in Design, Modelling, and Applications of Heat Transfer Equipment**

Special Issue Editors

**Zden ˇek Jegla Petr Stehl´ık**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Special Issue Editors* Zdenek Jegla ˇ Brno University of Technology Czech Republic

Petr Stehl´ık Brno University of Technology Czech Republic

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Energies* (ISSN 1996-1073) (available at: https://www.mdpi.com/journal/energies/special issues/advances in heat transfer equipment).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03936-481-7 (Hbk) ISBN 978-3-03936-482-4 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**

## **About the Special Issue Editors**

**Zden ˇek Jegla** (Assoc. Prof. Dr.) is currently an Associate Professor and Senior Researcher in the field of Process Engineering at the Institute of Process Engineering of the Faculty of Mechanical Engineering (FME) at the Brno University of Technology (BUT) and a Head of the Heat Transfer Research Group at the NETME Research Centre FME BUT. He has many years of research and industrial experience, and he is author or co-author of numerous papers and contributions in international journals and conferences. His research and development as well as application activities are aimed at applied and enhanced heat transfer; waste heat recovery systems and equipment; process fired heaters, boilers and heat exchangers; simulation, optimization and CFD applications in process and power industry; process and equipment design and integration for energy savings and emissions reduction.

**Petr Stehl´ık** (Prof. Dr.) is a Professor and Director of the Institute of Process Engineering of the Faculty of Mechanical Engineering at the Brno University of Technology. He is the Vice President of the Czech Society of Chemical Engineers and member of AIChE and ASME. He is an executive editor or guest editor of several international journals; member or chairmen of scientific committees of international conferences; reviewer of PhD theses at reputable foreign universities; coordinator or contractor of international and national research projects; and the author and co-author of numerous publications and several patents. His research, development, and application activities involve heat transfer and its applications (heat exchangers, combustion systems, waste heat utilization); thermal treatment and energy utilization of waste (waste-to-energy); and energy savings and environmental protection.

## **Preface to "Advances in Design, Modelling, and Applications of Heat Transfer Equipment"**

Heat-transfer equipment, typically represented by, for example, heat exchangers, process furnaces, and steam boilers, is among the essential equipment used for production processes in a number of industries (e.g., chemical and petrochemical, food, pharmaceutical, power, aviation and space) as well as for processes and applications in the communal sphere (e.g., waste incineration plants, heating plants, laundries, hospitals, server rooms, agriculture applications). Increasing demands for economical and efficient heat energy management can only be met when not only the layout of the whole system but also the individual heat-transfer equipment and its details are designed according to state-of-the-art knowledge. The purpose of this Special Issue is to present the latest advances in designing, modeling, testing, and operating heat-transfer equipment, including unconventional and innovative designs of heat-transfer equipment and their applications.

> **Zden ˇek Jegla, Petr Stehl´ık** *Special Issue Editors*

## *Review* **Review of Design and Modeling of Regenerative Heat Exchangers**

#### **Bohuslav Kilkovský**

Institute of Process and Environmental Engineering, Faculty of Mechanical Engineering, Brno University of Technology, Technická 2, 616 69 Brno, Czech Republic<EMAIL_ADDRESS>
Received: 28 December 2019; Accepted: 5 February 2020; Published: 9 February 2020

**Abstract:** Heat regenerators are simple devices for heat transfer, but their proper design is rather difficult. Their design is based on differential equations that need to be solved. This is one of the reasons why these devices are not widely used. There are several methods for solving them that were developed. However, due to the time demands of calculation, these models did not spread too much. With the development of computer technology, the situation changed, and these methods are now relatively easy to apply, as the calculation does not take a lot of time. Another problem arises when selecting a suitable method for calculating the heat transfer coefficient and pressure drop. Their choice depends on the type of packed bed material, and not all available computational equations also provide adequate accuracy. This paper describes the so-called open Willmott methods and provides a basic overview of equations for calculating the regenerative heat exchanger with a fixed bed. Based on the mentioned computational equations, it is possible to create a tailor-made calculation procedure of regenerative heat exchangers. Since no software was found on the market to design regenerative heat exchangers, it had to be created. An example of software implementation is described at the end of the article. The impulse to create this article was also to broaden the awareness of regenerative heat exchangers, to provide designers with an overview of suitable calculation methods and, thus, to extend the interest and use of this type of heat exchanger.

**Keywords:** regenerative heat exchanger; packed bed; heat transfer; pressure drop

#### **1. Introduction**

Regenerative heat exchangers are devices used for indirect heat exchange between hot and cold media. In these devices, heat is first transferred from a hot medium to a storage material and then is transferred to a cold medium. Thus, the hot and cold media are alternately in contact with the solid material forming the packed bed. In the hot cycle, heat is transferred from the medium to the packed bed, and, in the cold cycle, the cold medium absorbs the heat stored in the solid material. This cycle is the reason why regenerative heat exchangers must operate in pairs (they must have two beds) to work continuously. Regenerative heat exchangers are used mainly in the metallurgical industry, in air treatment, air preheating, or recovery of waste heat, and in turbine applications. However, the complexity of the calculations resulted in their limited expansion.

Although there are several approaches to calculating these devices, more accurate methods require the solution of differential equations. Their solution is quite time-demanding and they need the use of a computer. This paper describes the solution of the calculation of regenerative heat exchangers using the Willmott open method, which seems to be the best for computer use. This method shows great stability (convergence) and allows the inclusion of calculation of the equations describing heat transfer and pressure drop. The solution is performed by an iterative calculation, and the result is the distribution of gas temperatures and pressure over time and along the bed.

Another difficulty when calculating these devices is finding suitable computational equations to describe the heat transfer between the gas and the packed bed material. Although these equations can be found in the literature, they are usually given only for spherical shapes. In addition, different equations give different results. It is, therefore, necessary to choose a suitable calculation equation in order to get the best results. A better situation is in the calculation of pressure drops, for which many computational equations can be found in the literature. However, some of them give different results from those measured. This paper, therefore, provides an overview of suitable computational equations so that a complete computational algorithm of regenerative heat exchangers can be compiled.

#### **2. Description of the Regenerator**

#### *2.1. Classification of Regenerative Heat Exchangers*

Regenerators can be divided into two categories: fixed-bed and rotary regenerators. In the fixed-bed regenerator, a single fluid stream has cyclical or reversible flow. Valves are employed to switch the flow the hot and cold gas streams. In the rotary regenerator, the storage material rotates continuously through two counterflowing streams of fluids. Only one stream flows through a section of the storage material at a given time. However, both streams eventually flow through all sections of the storage material during one rotation.

Fixed-bed regenerators are commonly run in pairs. It means that two or more regenerators are used in parallel because of the requirement for a continuous stream of the gas. During one part of a cycle, the hot gas flows through one of the regenerators and heats up the storage material, while the cold gas flows through and cools down the storage material in the second regenerator. Both gases directly contact storage material in the regenerators, although not both at the same time, since each is in a different regenerator at any given time. After a sufficient amount of time, the cycle is switched such that the cooler storage material in the second regenerator is preheated with the hot gas, while the hot storage material in the first regenerator exchanges its heat into the cold gas. This cycle is permanently repeated.

The advantage of regenerators over recuperators is that they have a much higher surface area for a given volume. Hence, the regenerator usually has a smaller volume and weight than an equivalent recuperator. This means that regenerators are more economical in terms of materials and manufacturing. The storage material of regenerators also has a degree of self-cleaning characteristics, reducing fluid-side fouling and corrosion. Disadvantages include mixing the media as a result of alternating the passage of hot and cold media through the packed bed. Regenerators are, thus, ideal for gas–gas heat exchange.

Various materials and shapes can be used as storage materials. Because solids have a very large heat capacity compared to gases, they are used as intermediary storage of the heat. Their selection depends on given conditions and requirements, especially temperature. For very high temperature, ceramic storage material should be used. For low or moderate temperatures, the heat storage material can be made of metal, e.g., steel or aluminum. There exist several types of storage shapes (see Figure 1). For large regenerators, bricks can be used. For smaller regenerators, honeycombs, spherical particles, monoliths, saddles, rings, or Raschig rings can be used.

**Figure 1.** Some types of geometry of storage materials [1].

Regenerative heat exchangers can be used in various processes. The most common applications include the following:

The benefit of the exchanger is seen in the potential of its current multiple function. The exchanger could be used, for example, for simultaneous gas purification. This means that the storage material also serves as a catalyst on which the chemical reaction takes place.

#### *2.2. Basic Geometric Characteristics of Packed Bed*

Packed bed calculations use various geometric characteristics describing the packed bed. The main ones are specified in this section.

#### **Voidage, -**

An important parameter in the calculation of the flow in a packed bed is voidage. The voidage ε is defined as the ratio of the free volume of the packed bed to the total volume.

$$
\varepsilon = \frac{V\_b - V\_p}{V\_b} \times 100 = \frac{V\_m}{V\_b} \times 100,\tag{1}
$$

where *Vb* is the total volume of the packed bed (m3), *Vp* is the volume of the packed bed material (m3), and *Vm* is the free volume of the packed bed (m3).

#### **Particle diameter, m**

The particle diameter can be defined as the diameter of a sphere that has the same volume as the particle,

*dV* = - 6 π*Vp* 1/3 , (2)

or as the equivalent particle diameter *dp* (m), according to the specific surface given by Ergun [2], which has the same ratio of the surface to the volume as the given particle and is given by

$$d\_p = \frac{\theta \sum V\_p}{A\_p} \,\tag{3}$$

where *Ap* is the particle surface area (m2).

#### **Sphericity, -**

Sphericity is defined in Reference [3] as the ratio of the surface area of the sphere to the surface area of the particle. The sphericity is 1 for a sphere and is less than 1 for any particle that is not a sphere.

$$
\psi = \frac{A\_s}{A\_p} = \left[\frac{36\pi V\_p^2}{A\_p^3}\right]^{1/3} = \frac{\pi d\_p^2}{A\_p},\tag{4}
$$

where *As* is the surface area of a sphere that has the same volume as the particle (m2), *Ap* is the particle surface area (m2), and *Vp* is the volume of the particle (m3).

**Hydraulic diameter of packed bed, m**

$$d\_h = 4r\_h = 4\frac{\varepsilon}{a} = \frac{4\varepsilon}{a\_r(1-\varepsilon)},\tag{5}$$

where *rh* is the hydraulic radius (m), *a* is the absolute specific surface (m2), and *ar* is the relative specific surface (m−1).

#### **Absolute specific surface, (m**−**1)**

Absolute specific surface is the ratio between the particle surface area and the volume of the packed bed.

$$a = \frac{A\_p}{V\_b} = \frac{A\_p(1-\varepsilon)}{V\_p}.\tag{6}$$

#### **Relative specific surface, (m**−**1)**

Relative specific surface is the ratio between the particle surface area and the volume of the particle in the packed bed.

$$a\_l = \frac{A\_p}{V\_p}.\tag{7}$$

The relationship between relative and absolute specific surface is

$$a = a\_l (1 - \varepsilon). \tag{8}$$

#### **3. Mathematical Model of the Regenerative Heat Exchanger**

#### *3.1. Energy Balance of the Regenerator*

If the gas flows through the packed bed of the regenerator and the total heat transfer area between the bed and the gas is *Ap*, then the mean temperature change of packed bed *Tb* at time *t* can be expressed in the form

$$M\_b \mathcal{C}\_{p,b} \frac{\partial T\_b}{\partial t} = h\_t A \left( T\_\mathcal{g} - T\_b \right) \tag{9}$$

where *Mb* is the mass of packed bed (kg), *Tb* is the mean temperature of packed bed (◦C), *Tg* is the temperature of the gas flowing through the bed (◦C), *ht* is the total heat transfer coefficient (W·m−2·K−1) between the flowing gas and the bed material, *Cp*,*b* is the heat capacity of packed bed (J·kg−1· K−1), and *A* is the total heat transfer area (m2). In a cooling period, where the gas temperature, *Tg,* is lower than the bed temperature, *Tb*, the gas temperature increases over time while the bed temperature decreases *dTb dt* < 0. During the heating period (*Tg* > *Tb*), the gas outlet temperature decreases with time, while the bed temperature increases *dTb dt* > 0.

Heat is recovered or absorbed by the flowing gas through the packed bed of the regenerator. Since the gas flowing through the regenerator changes its temperature over time, we consider a change in the *y*-axis (along the height of the regenerator).

$$m\_{\mathcal{S}} \mathbf{C}\_{p,\mathcal{S}} L \frac{\partial T\_{\mathcal{S}}}{\partial y} + M\_{\mathcal{S}} \mathbf{C}\_{p,\mathcal{S}} \frac{\partial T\_{\mathcal{S}}}{\partial t} = h\_t A \{T\_b - T\_{\mathcal{S}}\}\_\prime \tag{10}$$

where *mg* is the mass flow rate of gas (kg·s−1), *Cp*,*g* is the heat capacity of gas (J·kg−1· K−1), *Mg* is mass of gas resident in the regenerator (kg), and *L* is the height of regenerator (m).

The most important assumption in this model is that the thermal conductivity of the packing material is infinite in a direction perpendicular to gas flow (and zero in a direction parallel to the gas flow). This implies that, at any level in the regenerator, the solid material is isothermal in a direction perpendicular to gas flow, and this may be true or approximately true where the packing is thin or is made of materials of high conductivity. In this case, the coefficient *ht* is the surface heat transfer coefficient, usually a convective coefficient to which may be added a radiative component.

However, if the packing of the regenerator is constructed of material of low thermal conductivity and/or the thickness of packing around the channels through which the gases flow is comparatively large, then it is necessary to incorporate the resistances to heat transfer at the solid surface and within the solid into a lumped or total heat transfer coefficient. Hausen [4] developed an equation to calculate this heat transfer coefficient in the following form:

$$\frac{1}{h\_{\rm th}} = \frac{1}{h\_{\rm lum}} + \frac{1}{h\_{\rm r}} = \frac{1}{h\_{\rm c}} + \frac{d}{2(n+2)\lambda\_{\rm b}}\phi\_{\rm H} + \frac{1}{h\_{\rm r}}.\tag{11}$$

where *n* = 1 for slabs (plane walls) of thickness *d* in (m), *n* = 2 for solid cylinders of diameter *d* in (m), and *n* = 3 for spheres of diameter *d* in (m), λ*b* is the thermal conductivity of packing material of regenerator (W·m−1·K−1), *hlum* is the lumped heat transfer coefficient (W·m−2·K−1), *hc* is the convective heat transfer coefficient (W·m−2·K−1), and *hr* is the radiative heat transfer coefficient (W·m−2·K−1). The lumped heat transfer coefficient incorporates the surface convective heat transfer coefficient, *hc*, and the resistance to heat transfer within the regenerator packing, as represented by the *d* 2(*n*+2)*k*φ*H* therm. The total heat transfer coefficient can be used in the conventional model of the thermal performance of the regenerator, set out in the differential equations.

The function φ*H*, called Hausen factor, attempts to reproduce the effect of the very rapid temperature changes within the packing, immediately after a reversal, at the start of a hot or cold period. According to Reference [5], this factor can be calculated in the case that *d*2 4α 1 *P* + 1 P  ≤ 5(*n* + 1)/2 using equation

$$\phi\_H = 1 - \frac{d^2}{4\alpha \left(n+3\right)^2 - 1} \left\{ \frac{1}{P'} + \frac{1}{P''} \right\} \tag{12}$$

and, for other values,

$$
\phi\_H = \pi (n+2) / \sqrt{\left(\varepsilon + \frac{d^2}{4\alpha} 18 \{\frac{1}{\mathbf{P'}} + \frac{1}{\mathbf{P''}}\}\right)} \tag{13}
$$

where α is the thermal diffusivity (m2·s−1), *P* is the length of period for heating and cooling process (s), ε = 2.7 for plates, ε = 9.9 for cylinders, and ε = 27.0 for spheres, Ω is the reduced time for hot period, and Ω" is the reduced time for cold period.

This problem is described in more detail in Reference [5].

#### *3.2. Di*ff*erential Equations*

Equations (9) and (10) were rearranged by Hausen [6] to the form

$$\frac{\partial T\_{\mathcal{S}}}{\partial \underline{\mathcal{S}}} = T\_b - T\_{\mathcal{S}'} \tag{14}$$

and

$$\frac{\partial T\_b}{\partial \eta} = T\_\S - T\_{b\prime} \tag{15}$$

where ξ is the dimensionless length, and η is the dimensionless time.

$$\eta = \frac{h\_t A}{M\_b \mathbb{C}\_{p,b}} \Big( t - \frac{M\_\%}{m\_\% L} y \Big),\tag{16}$$

and

$$\xi = \frac{h\_t A}{m\_\wp \mathbb{C}\_{p,\xi} L} \mathcal{Y}.\tag{17}$$

When *t* = *P* and *y* = *L*, each period of regenerator operation is defined in terms of two dimensionless parameters given by Hausen [6] "reduced period", Π, and "reduced length", Λ.

$$
\Pi = \frac{h\_t A}{M\_b \mathbb{C}\_{p,b}} \Big| P - \frac{M\_\mathcal{S}}{m\_\mathcal{S}} \Big) \tag{18}
$$

and

$$
\Lambda = \frac{h\_{\text{fl}} A}{m\_{\text{g}} \mathbb{C}\_{p, \text{g}}}.\tag{19}
$$

The effectiveness of regenerator behavior may be measured in terms of the thermal ratios, η*REG*. The thermal ratio for the heating period is η *REG* and that for the cooling period is η *REG*. Ideally, the exit gas temperature in the hot or cold period should be equal to the inlet gas temperature in the opposite period. The thermal ratios, which are defined below, measure the degree to which this ideal is achieved.

For the heating period,

$$\eta'\_{\rm REC} = \frac{T'\_{\rm g,o, \, m} - T'\_{\rm g,i}}{T'\_{\rm g,i} - T'\_{\rm g,i}},\tag{20}$$

and, for the cooling period,

$$\eta''\_{\rm REC} = \frac{T'\_{\rm g.i} - T'\_{\rm g.o.m}}{T'\_{\rm g.i} - T'\_{\rm g.i}}.\tag{21}$$

*Energies* **2020**, *13*, 759

The term thermal ratio can be misleading in the sense that it is not always a measure of efficiency, and perhaps the term temperature ratio might be more appropriate.

In fixed-bed regenerators, the exit gas temperatures vary with time. The chronological average exit temperatures are, therefore, computed (*T g*,*o*, *m* and *T g*,*o*,*m*).

For the symmetric regenerator (Λ = Λ = Λ", Π = Π = Π", and η*REG* = η'*REG* = η"*REG*), an estimate of the thermal ratio is given by

$$
\eta\_{\rm REC} = \frac{\Lambda}{\Lambda + 2}.\tag{22}
$$

It can be seen that a larger reduced length leads to a greater thermal ratio.

Let us suppose that we have a heat balance, once cyclic equilibrium is attained,

$$m\_{\mathcal{g}}' \mathbf{C}\_{p,\mathcal{g}}' P' \left(T\_{\mathcal{g},i}' - T\_{\mathcal{g},o}'\right) = m\_{\mathcal{g}}' \mathbf{C}\_{p,\mathcal{g}}' P'' \left(T\_{\mathcal{g},o}' - T\_{\mathcal{g},i}'\right). \tag{23}$$

Dividing this equation by - *T g*,*i* − *T g*,*i* , we get

$$m\_{\mathcal{g}}' \mathbf{C}\_{p,\mathcal{g}}' \mathbf{P}' \boldsymbol{\eta}\_{REG}' = m\_{\mathcal{g}}' \mathbf{C}\_{p,\mathcal{g}}' \mathbf{P}'' \boldsymbol{\eta}\_{REG}'.\tag{24}$$

If *m gC p*,*gP* = *m g C p*,*gP* , the regenerator is said to be balanced and both thermal ratios are equal. Equation (24) can be converted to the following form:

$$\frac{\Pi'}{\Lambda'} \eta'\_{\rmREG} = \frac{\Pi''}{\Lambda''} \eta''\_{\rmREG}.\tag{25}$$

We can say that, in general, a regenerator is balanced if

$$\frac{\Pi'}{\Pi''} = \frac{\Lambda'}{\Lambda''} = k.\tag{26}$$

If *k* = 1, the regenerator is said to be symmetric. Equation (25) can be modified to the following form:

$$\frac{\Pi'}{\Lambda'} \times \frac{\Lambda''}{\Pi''} = \mathcal{y}.\tag{27}$$

If γ = 1, the regenerator is balanced.

When *m gC p*,*gP m g C p*,*gP* , this corresponds to the most general case where γ - 1, η *REG* η *REG*, and the regenerator is said to be unbalanced. A summary of these classifications is presented in Table 1.

**Table 1.** Possible types of regenerators.

The previous equations apply if the regenerative heat exchanger is symmetric. Hausen [4] proposed that the performance of a balanced nonsymmetric regenerator can be accurately estimated using the symmetric regenerator model employing the harmonic reduced length Λ*H* and the harmonic reduced period Π*H* in both hot and cold periods.

$$\frac{2}{\Lambda\_H} = \frac{1}{\Pi\_H} \left( \frac{\Pi'}{\Lambda'} + \frac{\Pi'}{\Lambda''} \right) \tag{28}$$

and

$$\frac{2}{\Pi\_H} = \frac{1}{\Pi'} + \frac{1}{\Pi''}.\tag{29}$$

This proposal was verified as acceptable by Iliffe [7], and its applicability was extended to

unbalanced regenerators by Razelos [8]. The thermal ratio than can be calculated using Equation (22). A factor *K*/*K*0 describing the effect of nonlinear variations of temperature using a factor is given in Reference [4].

$$\frac{K}{\mathbb{K}\_0} = \frac{\eta\_{\rm REC}}{1 - \eta\_{\rm REC}} \frac{2}{\Lambda\_H}.\tag{30}$$

A smaller value of *K*/*K*0 results in a greater effect of both the nonlinear variations of temperature and the corresponding truncation error.

#### *3.3. Calculation Methods*

Several methods were developed to determine the solution of regenerative heat exchangers. Some of them are shown in Figure 3. These methods can be divided into two groups: rapid and precise methods. Rapid methods are intended only for fast and preliminary calculation and do not provide sufficient accuracy. In addition, only the media outlet temperatures are the result. These methods are not suitable for the application of real calculations and designs. Rapid methods were created mainly before the expansion of computers. An example of a rapid method is given in Reference [4].

**Figure 3.** Diagram of calculation methods used for the solution of regenerative heat exchangers [4,7,9–14].

The most suitable methods for the design of these heat exchangers are precise methods which involve the solution of differential equations, and the result is the variation of the media temperatures at the outlet over time. These methods can be classified into two groups: open and closed. Each of these two methods comprises two subgroups, namely, linear and nonlinear methods, according to the way in which thermophysical properties of the gas and the storage material are calculated.

In the open methods, the gas and solid temperatures are evaluated by solving differential equations over successive cycles of regenerator operation. The temperature profiles in each cycle and time are the result. We know the number of cycles to equilibrium at the end of the calculation. The closed methods are those in which the steady-state performance is calculated directly without the consideration of any previous cycles. In general, the closed methods are faster than open methods; however, when using modern computers, the difference is slight. The closed methods appear to be suitable for solving a linear problems, but are becoming extremely complicated (and sometimes unstable) in solving nonlinear problems, i.e., the thermophysical properties of both fluid and solids, including heat transfer coefficients, can vary spatially and temporally, depending on temperature and/or when mass flow rates of fluids in one or both periods of regeneration operation may vary with time. Open methods show great stability even when solving nonlinear problems. In these cases, open methods are preferable. Another advantage of open methods is their easy modification by including equations for calculating heat transfer coefficient and pressure drops.

For computer solutions, the open methods proposed by Willmott, given in References [9,10], seem to be the most suitable. These methods were chosen because they are relatively easy to program, and they involve the calculation of important parts (i.e., calculation of heat transfer coefficient, calculation of radiation influence, calculation of pressure losses, geometric characteristics). In addition, these methods enable including the effect of changes in the thermophysical properties of media flowing through the packed bed in time (nonlinear model) and are more stable.

#### *3.4. Selected Mathematical Model*

The selected Willmott open method can be further divided into a linear, quasi-linear, and nonlinear model. The difference between these models is in the calculation of the thermophysical properties of gas and bed.

The linear model calculates the thermophysical properties of the gas and bed based on the reference temperature, i.e., the properties are constant at each packed bed point and at every moment. The calculation is fast, but it may not be accurate enough in some cases, especially when the temperature changes significantly. The fluid properties are calculated at a reference temperature.

$$T\_{\rm g,ref} = \frac{T\_{\rm g,i}^{\prime} + T\_{\rm g,i}^{\prime\prime}}{2}. \tag{31}$$

The quasi-linear model was described in Reference [13]. This model uses a different reference temperature for hot and cold media.

$$T'\_{\rm g.ref} = \frac{T'\_{\rm g.i} + T'\_{\rm g.o}}{2}; \quad T''\_{\rm g.ref} = \frac{T''\_{\rm g.i} + T'\_{\rm g.o}}{2}. \tag{32}$$

The outlet temperatures are replaced by the newly calculated ones, and the reference temperatures are recalculated after each period or cycle. Based on these temperatures, all thermophysical properties of fluids, packed bed material, and related values are recalculated.

The nonlinear model considers the change in gas and bed properties at the place and time as a function of temperature. The heat transfer coefficient is then calculated from these properties at every moment and place. This is important for the most realistic simulation of high-temperature regenerators. The calculation is much more accurate but time-consuming. However, this is not a considerable problem with up-to-date computers.

Reference [9] proposed a method for solving Equations (14) and (15) so that it could be used for the solution on computers. This method uses a trapezoidal method for the numerical solution of differential equations. This method of calculating the regenerator is known as the Willmott open method. The simplifications introduced to the derivation and calculation of these differential equations are as follows, according to Reference [9]:

#### **Boundary conditions**

There are determined two boundary conditions.

$$T\_b'(y,0) = T\_b'^\prime(L - y\_\prime P^{\prime\prime}),$$

$$T\_b'^\prime(y,0) = T\_b'(L - y\_\prime P^\prime).$$

For the first cycle, the temperature along the bed is set arbitrarily, e.g., to ambient temperature. In the following cycles, the mentioned condition is already applied.

Using the trapezoidal numerical method, which provides excellent properties of numerical stability, the differential equations move to the following forms:

$$T\_{\mathcal{g}\,\,r+1,\mathcal{S}} = T\_{\mathcal{g}\,\,r\,\mathcal{S}} + \frac{\Delta\xi}{2} \left\{ \left(\frac{\partial T\_{\mathcal{S}}}{\partial\xi}\right)\_{r+1,\mathcal{S}} + \left(\frac{\partial T\_{\mathcal{S}}}{\partial\xi}\right)\_{r,\mathcal{S}} \right\}\_{\mathcal{r},\mathcal{S}} \tag{33}$$

and

$$T\_{b\ r,S+1} = T\_{b\ r,s} + \frac{\Delta \eta}{2} \left\{ \left(\frac{\partial T\_b}{\partial \eta}\right)\_{r+1,S} + \left(\frac{\partial T\_b}{\partial \eta}\right)\_{r,S} \right\}\_{r,S} \tag{34}$$

where *r* refers to distance (*r* = 0..*m*) and *S* refers to the time (*S* = 0..*P*). Thus, the distance step is *L*/*m* and the time step is mostly 1 s.

#### 3.4.1. Linear Model

In the case of the linear model, these equations can be further adjusted to form

$$T\_{\mathbb{g}\,\,r+1,\mathbb{S}} = A\_1 T\_{\mathbb{g}\,\,r,\mathbb{S}} + A\_2 \Big(T\_{\mathbb{b}\,\,r+1,\mathbb{S}} + T\_{\mathbb{b}\,\,r,\mathbb{S}}\Big),\tag{35}$$

and

$$T\_{b\ r,S+1} = B\_1 T\_{b\ r,S} + B\_2 (T\_{\ll r,S+1} + T\_{\ll r,S}) \, \tag{36}$$

where

$$A\_1 = \frac{1-\alpha}{1+\alpha},\ A\_2 = \frac{\alpha}{1+\alpha},\ B\_1 = \frac{1-\beta}{1+\beta},\ B\_2 = \frac{\beta}{1+\beta},$$

$$\alpha = \frac{1}{2}\Delta\xi = \frac{\Lambda}{2m},\ \beta = \frac{1}{2}\Delta\eta = \frac{\Pi}{2P}.$$

In these equations, Λ is the reduced length, Π is the reduced time, *m* is the number of sections, and *P* is the cycle time.

These equations are the main equations needed to calculate the regenerative heat exchanger. For the proper solution of these equations, the bed needs to be divided into a suitable number of sections, in which the temperature calculations are gradually performed. The calculation is then a combination of the above and other supplementary equations. The calculation is solved by an iterative procedure for the given time over the entire length of the bed. This procedure is then repeated for all times and all periods (cooling and heating) until equilibrium is reached, i.e., for consecutive cooling cycles, abs (Φ(*n*) − Φ(*n* − 1)) must be less than the value given by the user. Then, the cycle (*n* + 1) is the equilibrium cycle.

$$\Phi(n) = \frac{T'\_{\mathbb{g},o,\,m} - T'\_{\mathbb{g},i}}{T'\_{\mathbb{g},i} - T'\_{\mathbb{g},i}}.\tag{37}$$

This is not an exhaustive description of the solution of differential Equations (33) and (34) using a trapezoidal rule for a linear model. See Reference [9] for more information how to apply this method to the computer. This basic computational model was developed in GNU Octave software and then applied to several examples differing in reduced period value and reduced length. Subsequently, the influence of these reduced quantities on the calculation was evaluated.

**Case 1—taken from Reference [9]**

The input data of the calculation are given in Table 2. The task was to calculate the outlet temperature of the hot and cold medium for the different number of sections *m*. If we divide the regeneration bed into a small number of sections, the outlet temperature does not correspond to reality (see the results).

**Table 2.** Input data for calculation.

The results of the calculation are shown in graphical form in Figure 4.

**Figure 4.** Influence of number sections on outlet temperature obtained by the basic calculation model: (**a**) for the heating period; (**b**) for the cooling period.

These calculations showed that the accuracy of the outlet temperature is dependent on the number of sections on which the regenerator (packed bed) is divided. In this case, the limit number of sections is 10. However, this is not true in all cases. From further calculations, it was possible to prove that, with a higher value of reduced length, the required number of sections on which the packed bed must be divided decreases; however, at the same time, the number of cycles needed to reach the equilibrium increases. On the other hand, as the reduced period increases, the number of sections (steps) increases, but the number of cycles necessary to establish equilibrium decreases. I suggest having a minimum of 50 sections. This is no problem for computers.

#### 3.4.2. Nonlinear Model

The linear model considers only the constant properties of fluids. This means faster calculation; however, in some cases, the results may differ significantly from reality. Such a procedure is only suitable for preliminary calculations. More accurate results are obtained when changing the properties of fluid with a change in temperature. This model is called the nonlinear model.

In the case of the nonlinear model, Equations (33) and (34) can be further adjusted to form

$$T\_{\mathbb{S}^r}r + 1,\mathbb{S} = A\_{1,r,\mathbb{S}}T\_{\mathbb{S}^r,\mathbb{S}} + A\_{2,r+1,\mathbb{S}}T\_{\mathbb{b}^r,r+1,\mathbb{S}} + A\_{3,r,\mathbb{S}}T\_{\mathbb{b}^r,\mathbb{S}^r} \tag{38}$$

and

$$T\_{b\ r\_{r}S+1} = B\_{1,r,S}T\_{b\ r\_{r}S} + B\_{2,r,S+1}T\_{\text{g. }r\_{r}S+1} + B\_{3,r,S}T\_{\text{g. }r\_{r}S} \tag{39}$$

where

$$A\_{1,r,S} = \frac{1 - \alpha\_{r,S}}{1 + \alpha\_{r+1,S}}, \ A\_2 = \frac{\alpha\_{r+1,S}}{1 + \alpha\_{r+1,S}}, \ A\_{3,r,S} = \frac{\alpha\_{r,S}}{1 + \alpha\_{r+1,S}},$$

$$B\_{1,r,S} = \frac{1 - \beta\_{r,S}}{1 + \beta\_{r,S+1}}, \ B\_{2,r,S+1} = \frac{\beta\_{r,S+1}}{1 + \beta\_{r,S+1}}, \ B\_{3,r,S} = \frac{\beta\_{r,S}}{1 + \beta\_{r,S+1}},$$

$$\alpha = \frac{1}{2} \Delta \xi = \frac{\Lambda}{2m'}, \ \beta = \frac{1}{2} \Delta \eta = \frac{\Pi}{2P}.$$

These constants are calculated in every point and time with changing temperature. This is not an exhaustive description of the solution of differential Equations (33) and (34) using a trapezoidal rule for a nonlinear model. See Reference [10] for more information how to apply this method to the computer.

#### **4. Pressure Drops**

A very important value in designing a regenerative heat exchanger or, generally, a device using a packed or structured material layer is the amount of pressure drop of the flowing gas. The general equation for its calculation is mostly defined as

$$
\Delta p = \lambda\_k \frac{L}{d\_p} \overline{c}^2 \times \frac{(1-\varepsilon)}{\varepsilon^3} \times \rho\_\prime \tag{40}
$$

where *c* is the velocity of gas based on the empty cross-section of the bed (m·s−1), and λ*k* is the friction factor, most often given in the form

$$
\lambda\_k = \frac{k\_1}{\mathcal{R}c\_m} + \frac{k\_2}{\mathcal{R}e\_m^b}'\tag{41}
$$

where *k*1 and *k*2 are constants, and *b* is the exponent. *Rem* is the modified Reynolds number given by Ergun as follows:

$$\mathrm{Re}\_{\mathrm{ll}} = \frac{d\_{\mathbb{P}} \rho \tilde{\mathbf{c}}}{\eta (1 - \varepsilon)} = \frac{\mathrm{Re}}{(1 - \varepsilon)}. \tag{42}$$

Pressure drop in the packed bed is commonly calculated using the Ergun equation [2].

$$
\lambda\_k = \frac{-\Delta P \, d\_p}{L \rho \overline{\varepsilon}^2} \frac{\varepsilon^3}{(1 - \varepsilon)} = \frac{150}{R \varepsilon\_m} + 1.75. \tag{43}
$$

According to Reference [15], this equation very often over-predicts the respective value, and it is better to use special equation tailored to specific storage materials which are based on experimental data.

Ergun also defined the friction coefficient in the shape

$$
\lambda\_v = \frac{-\Delta P \, d\_p^2}{\mu L \tilde{\varepsilon}} \frac{\varepsilon^3}{\left(1 - \varepsilon\right)^2} = \lambda\_k \frac{\text{Re}}{1 - \varepsilon} = \lambda\_k \text{Re}\_m. \tag{44}
$$

Sometimes, friction coefficient proposed in Reference [16] is used.

$$
\lambda\_p = \frac{-\Delta P \, d\_p}{L \rho \overline{\varepsilon}^2} = \, \lambda\_k \frac{(1-\varepsilon)}{\varepsilon^3} = \, \lambda\_v \frac{(1-\varepsilon)^2}{\varepsilon^3 \, Re}. \tag{45}
$$

The Reynolds number is often also given in a modified form, which is defined by the equation

$$Re\_l = \frac{d\_l \rho \tilde{c}}{6\eta (1 - \varepsilon)} = \frac{Re}{6(1 - \varepsilon)}.\tag{46}$$

There are many correlations for calculating the friction coefficient. The best known is the Ergun equation mentioned above, but it overestimates pressure drop for a bed of randomly arranged smooth spheres for *Rem* > 700.

Based on an extensive comparison of the computational equations for determining the pressure drop through the packed bed performed in Reference [17], the computational equations listed in Table 3 were selected.

**Table 3.** Suitable equations for calculating of friction coefficient for the packed bed of spheres.

There are several other computational equations. These are mentioned, for example, in Reference [15], which deals with the effect of particle shape, size distribution, packing arrangement, and roughness of particles on pressure drop.

The Eisfeld and Schnitlein equation (see Table 3) using the constants *K*1 = 190, *k*1 = 2.00, and *k*2 = 0.77 can be used for the cylindrical shapes of the particles. For other particles, the constants *K*1 = 155, *k*1 = 1.42, and *k*2 = 0.83 can be used.

For non-spherical particles, Reference [25] proposed altering the constants in the Ergun equation by means of the particle sphericity ψ. For cylindrical particles, the friction factor is altered to

$$
\lambda\_k = \frac{150}{\psi^{3/2} Re\_m} + \frac{1.75}{\psi^{4/3}}.\tag{47}
$$

This equation is based on experimental data for *Rem* < 400.

Singh et al. [26] presented a correlation for pressure drop through beds of differently shaped particles in the form

$$
\lambda\_k = \frac{\varepsilon^3}{(1-\varepsilon)} 4.466 Re^{-0.2} \psi^{0.696} \varepsilon^{-2.945} \mathcal{e}^{11.85 \left(\log \psi\right)^2}.\tag{48}
$$

This is not an exhaustive list of all available computational equations to determine the friction coefficient as the medium flows through the packed bed. Other equations can be obtained in References [15,17].

For calculating the pressure drop across a layer consisting of a set of straight channels, e.g., honeycomb, the formulas for calculating pressure drop in a straight channel can be used.

$$
\Delta p = \lambda \frac{\overline{\mathcal{E}}}{2\varepsilon} \frac{h}{d\_h} \rho\_\prime \tag{49}
$$

where the friction factor can be calculated using the Swamee–Jain equation [27].

$$\lambda = \frac{0.25}{\left[\log\left(\frac{\frac{\epsilon}{d\_h}}{3.7} + \frac{5.74}{\mathrm{Re}^{0.9}}\right)\right]^2} \tag{50}$$

where *e* is the effective roughness height (m), and *dh* is the hydraulic diameter (m).

An important value is the pressure drop of the sieve on which the bed is placed. This pressure drop must also be determined and included in the total pressure drop.

#### **5. Voidage Calculation**

The pressure drop is strongly influenced by the mean voidage in the packed bed. In industrial practice, a large number of different non-spherical shaped particles are used, which form the bed and have different voidage. The most cited and experimentally verified effect on voidage is the ratio of regenerator diameter to equivalent particle diameter. This effect is most significant for a ratio of diameters less than 10 (see Figure 5). If the bed consists of non-spherical particles, their diameter is expressed by the equivalent diameter *dp*. In the randomly packed bed, other than spherical particles, the influence of particle orientation affects the local distribution of the voidage, which affects the mean voidage. This means that different medium voidage can be achieved with each filling of the packing.

**Figure 5.** Dependence of mean voidage on the ratio *D*/*dp*.

Among the most used empirical formulas for calculating the voidage are those from References [28–32]. These equations tend to be awkward and have too many parameters. Due to the complexity of using the mentioned equations, simpler but sufficiently accurate equations were developed to calculate the voidage. Some of them are listed below.

The equation designed by Zou and Yu [33] and recommended by Di Felice and Gibilaro [34] is in the form

$$\varepsilon = \varepsilon\_b + 0.01 \left[ \exp\left(\frac{10.686}{D/d\_p}\right) - 1\right] \tag{51}$$

for *dp D* ≤ 0.256.

By analyzing experimental data, Zou and Yu suggested a coefficient of ε*b* = 0.4, while Reference [35] suggested that ε*b* = 0.373, a value derived from their own data (subscript *b* means a bulk region).

Benyahia et al. [36] developed the following equations based on the measured data:

For spherical particles (1.5 ≤ *D*/*dp* ≤ 50),

$$
\varepsilon = 0.390 + \frac{1.740}{\left(D/d\_p + 1.140\right)^2}.\tag{52}
$$

For solid cylinders (1.7 ≤ *D*/*dp* ≤ 26.3),

$$
\varepsilon = 0.373 + \frac{1.703}{\left(D/d\_p + 0.611\right)^2}.\tag{53}
$$

For hollow cylinders (1.9 ≤ *D*/*dp* ≤ 14.5),

$$
\varepsilon = 0.465 + \frac{2.030}{\left(D/d\_p + 1.033\right)^2}.\tag{54}
$$

Pro general particles (1.5 ≤ *D*/*dp* ≤ 50; 0.42 < ψ < 1.0),

$$
\varepsilon = \left(0.1504 + \frac{0.2024}{\psi}\right) + \frac{1.0814}{\left(D/d\_p + 0.1226\right)^2}.\tag{55}
$$

The course of the voidage versus the *D*/*dp* ratio is shown in Figure 5. It can be seen that, up to value 10, the voidage strongly depends on this ratio. From this value, the voidage is approximately constant.

#### **6. Heat Transfer Calculation**

Determination of the heat transfer coefficient of the gas flowing through the packed bed is the most sensitive point of the calculation. Even if we set the correct number of sections and consider the variable properties of the fluids and the packed bed material, if we determine the wrong heat transfer coefficient, we get poor results.

The heat transfer coefficient mostly includes the effect of both convection and radiation for both streams. Heat transfer coefficient depends predominantly on the type of storage material. The respective equations for convective heat transfer were published, e.g., in Reference [37], while those for radiative heat transfer were discussed, e.g., in Reference [38].

The heat flux density *q* between the gas and the packed bed material can be expressed as follows:

$$q = q\_c + q\_r = \overline{h\_t} (T\_\mathcal{g} - T\_b)\_\prime \tag{56}$$

where *qc* is the convective heat flux density (W·m−2), *qr* is the radiative heat flux density (W·m−2), and the total heat transfer coefficient is defined as

$$h\_t = h\_\mathcal{c} + h\_\mathcal{r}.\tag{57}$$

The influence of the conductive transfer is negligible compared to the others. As mentioned previously, the convective component must be, in some cases, replaced by the effective heat transfer coefficient *hlum*.

$$h\_t = h\_{lmu} + h\_r.\tag{58}$$

#### *6.1. Convective Heat Transfer Coe*ffi*cient*

The choice of the heat transfer coefficient is a crucial factor. Many computational equations can be found to calculate the heat transfer coefficient in the literature. However, they can give different results. It is necessary to choose a computational equation that was verified for the particles that form the packed bed.

The equation for a randomly packed bed (Equation (59)), which was verified on experimental data, was recommended in Reference [39].

$$Nu = \frac{h\_t d\_p}{\lambda\_\text{g}} = 2 + 1.8 Re^{1/2} Pr^{1/3} \,\text{s} \tag{59}$$

for *Pr* = 0.7–0.8, *Re* > 100. Here, λ*g* is the thermal conductivity (W·m−1·K−1), *Re* is the Reynolds number (-), and *Pr* is the Prandtl number (-).

The heat transfer coefficient in the packed bed of spheres can be calculated according to Reference [40] using the following equation:

$$Nu = \frac{h\_{\text{c}}d\_p}{\lambda\_{\text{g}}} = 0.584 Re^{0.7} Pr^{1/3} \text{\AA} \tag{60}$$

for *Re* = 500–50,000. This equation was also recommended in Reference [41]. A significant wall effect was observed at ratios of bed diameter to particle diameter less than 20. Above that ratio, this equation is valid.

An equation for a packed bed formed with spheres was recommended in Reference [42].

$$Nu = \frac{h\_t d\_p}{\lambda\_\mathcal{g}} = 1.09 Re^{0.68} Pr^{1/3} \,\text{s} \tag{61}$$

for *Re* = 200–10,400.

Another suitable equation was mentioned in Reference [40].

$$Nu = \frac{h\_{\varepsilon}d\_{p}}{\lambda\_{\mathcal{g}}} = \varepsilon Re Pr^{1/3} \left( 0.0108 + \frac{0.929}{Re^{0.58} - 0.483} \right) \tag{62}$$

This equation takes into account the voidage, ε, and it is valid for a Reynolds number of 20–10,000. The equation for flow in a packed bed can also be used [43].

$$Nu = \frac{h\_{\varepsilon}d\_{p}}{\lambda\_{\text{g}}} = \left(0.5 Re^{0.5} + 0.2 Re^{2/3}\right) Pr^{1/3},\tag{63}$$

for 20 < *Re* < 100,000.

As mentioned in Reference [44], the heat transfer performances in a packed bed are usually formulated by traditional correlation as follows:

$$Nu = \frac{h\_{\text{eff}} d\_p}{\lambda\_{\text{g}}} = a\_1 + a\_2 Pr^{1/3} Re^n \left(\frac{d\_p}{d\_h} \varepsilon\right)^n \tag{64}$$

where *a*1, *a*2, and *n* are model constants. The values of these constants given by Reference [44] are *a*1 = 2.0, *a*2 = 1.1, and *n* = 0.6. Values of these constant for different packed cells were mentioned in Reference [45].

An experimental study was carried out to investigate the heat transfer of packed bed solar energy storage system having large-sized elements of storage material of different shapes [26].

$$Nu = \frac{h\_{\text{c}}d\_p}{\lambda\_{\text{g}}} = 0.437 \text{Re}^{0.75} \psi^{3.35} \varepsilon^{-1.62} e^{29.03 \left(\log \psi\right)^2} . \tag{65}$$

The heat transfer in the structured packed bed formed by straight channels can be determined using the following basic equation:

$$Nu = \frac{h\_{\text{c}}d\_p}{\lambda\_{\text{g}}} = 0.023 Re^{4/5} Pr^{1/3}.\tag{66}$$

#### *6.2. Radiative Heat Transfer Coe*ffi*cient*

In certain high-temperature applications, the heat exchanger is heated up by a flue gas containing significant proportions of CO2 and H2O vapor. In these circumstances, the radiation heat transfer must be considered.

The heat flux density due to gas–solid radiation is given by

$$q\_{\mathcal{I}} = \frac{\varepsilon\_{\mathcal{b}} + 1}{2} \sigma \Big[ \varepsilon\_{\mathcal{S}} T\_{\mathcal{S}}^{4} - \alpha\_{\mathcal{S}} T\_{\mathcal{b}}^{4} \Big],\tag{67}$$

where ε*g* and α*g* are the emissivity and mean absorptivity of gases (-), respectively, ε*b* is the emissivity of bed material (-), and σ is Stefan Boltzmann's constant, 5.67 × 10−8 (W·m−2·K−1). The fraction of ε*b*+1 2 is sometimes called the emissivity correction factor.

The values of emissivity and absorptivity are proportional to the percentage of carbon dioxide and water vapor in the gas and to the beam length. Both are functions of the gas temperature. These values can be obtained from Hottel's charts [46] or using the methods given in the literature.

$$q\_r = h\_r \Big(T\_\% - T\_b\rceil. \tag{68}$$

According to Reference [47], it is possible to calculate the radiative heat transfer *hr* using the equation mentioned in Reference [48].

$$
\hbar h\_r \approx 4 \sigma \varepsilon\_b T\_b^3. \tag{69}
$$

More specifically, however, the emissivity and gas absorption can be determined based on the mean beam length (m), *Lb*, in a given space and on the partial pressures of the respective radiant gas components. Several methods can be found in the literature.

The mean beam length *Lb* is determined from a known equation for the case when *Lb* > 1 m,

$$L\_b = \frac{3.6V}{P},\tag{70}$$

and, for the case when *Lb* < 1 m,

$$L\_b = \frac{3.4V}{P},\tag{71}$$

where *V* is the mean channel volume (m3), and *P* is the mean inner surface of the channel in the packed bed (m2).

The emissivity of gas can be calculated using the equation from Reference [49].

$$
\varepsilon\_{\mathcal{S}} = \left[ \varepsilon\_{CO\_2} + \left( \varepsilon\_{H\_2O} \mathbb{C}\_{H\_2O} \right) \right] (1 - \mathbb{C}\_{SO}),
\tag{72}
$$

where ε*CO*2 is the emissivity of carbon dioxide (-), ε*H*2*O* is the uncorrected emissivity of water vapor (-), *CH2O* is Beer's law correction factor for water vapor (-), and *CSO* is the spectral overlap correction factor (-).

The emissivity of carbon dioxide (ε*CO*2) is a function of temperature and the product (*pCO*2*L*), where *pCO*2 is the partial pressure of carbon dioxide in the gases (Pa) and *Lb* is the mean beam length (m). Thus,

$$
\varepsilon\_{\rm CO\_2} = f(p\_{\rm CO\_2} L\_{b\prime} T\_{\rm \%}).\tag{73}
$$

The emissivity of water vapor (ε*H*2*O*) is a function of temperature and (*pH*2*OLb*), where *pH*2*O* is the partial pressure of water vapor (Pa) and *Lb* is the mean beam length (m). Thus,

$$
\varepsilon\_{H\_{2\bullet}} = f(p\_{H\_2\bullet} L\_{b\star} T\_{\S}).\tag{74}
$$

Mean absorptivity of gas (α*g*) is estimated using the equation from Reference [49].

$$\alpha\_{\mathcal{S}} = \left(a\_{CO2} + a\_{H\_2O}\right) \left(1 - \mathcal{C}\_{\mathcal{SO}}\right) = \left(\left[\varepsilon\_{CO\_2} \left(\frac{T\_{\mathcal{S}}}{T\_b}\right)^{0.65}\right] + \left[\varepsilon\_{H\_2O} \mathcal{C}\_{H\_2O} \left(\frac{T\_{\mathcal{S}}}{T\_b}\right)^n\right]\right) \tag{75}$$

where *n* = 0.5 for *Tb* < 500 ◦C, *n* = 0.4 for *Tb* > 900 ◦C, and *n* = 0.45 for 500 ◦C < *Tb* < 900 ◦C.

The emissivity of carbon dioxide (ε*CO*2) and water vapor (ε*H*2*O*), and the absorptivity of gases can be computed using methods available from the literature [50–53].

#### *6.3. Heat Losses*

If the wall of the regenerator is not sufficiently insulated, the heat losses through the wall of the equipment should also be included in the calculation.

These are defined as follows:

$$Q\_w = h\_o A\_w (T\_w - T\_{inf})\_\prime \tag{76}$$

where *Aw* is the wall area (m2), *Tw* is the wall temperature (◦C), *Tinf* is the ambient temperature (◦C), and *ho* is the outside convective heat transfer coefficient (W·m−2·K−1).

For vertical cylinders, the convective heat transfer coefficient *ho* is calculated from a correlation of Churchill and Chu [54] for plane surfaces as

$$Nu = \frac{h\_0 H}{\lambda\_o} = \left\{ 0.825 + \frac{0.387 Ra^{\frac{1}{6}}}{\left[ 1 + \left( 0.492 / Pr \right)^{9/16} \right]^{8/27}} \right\}^2,\tag{77}$$

where *H* is the height of the regenerator (m), and *Ra* is the Rayleigh number given as follows:

$$Ra = \frac{\lg \beta (T\_w - T\_\infty) H^3}{a\nu}.\tag{78}$$

where *H* is the height of the regenerator (m), *g* is the gravity acceleration (m·s−2), β = 1 *Tf* = 2 *Tw*+*Tin f* is the thermal expansion coefficient (K−1), α is the thermal diffusivity (m2·s−1), and ν is the kinematic viscosity of gas (m2·s−1) as defined in the temperature of *Tf*. Validity of the correlation is in the range 10−1 < *Ra* < 1012. This correlation can be used provided the curvature effect is not too significant. This represents the limit where boundary layer thickness is small relative to cylinder diameter *D*. The correlations for vertical plane walls can be used when *D*/*L* ≥ *35*/*Gr0.25* where *Gr* is the Grashof number.

#### **7. Software Implementation of the Model**

Since no software was found on the market to design regenerative heat exchangers, it had to be created. The abovementioned Willmott methods for the linear and nonlinear model, together with mentioned equations for determination of heat transfer and pressure drop, were used in creating the computational software. To make the software user-friendly, the JAVA environment was used to build the software. The created calculation software (see Figure 6) enables effectively solving these types of heat exchangers, and it provides results in the form of text output and graphical dependencies of temperatures and pressure drops along the packed bed and over time. The software enables a user to print input data, text results, and various graphical dependencies. It is also possible to save these data and charts in various graphic file formats and copy dependencies in text form to Excel.

**Figure 6.** Screenshot of created software.

In order to perform the calculation, the user must specify the geometric characteristics of the packed bed, the media properties, the time of hot and cold period, the number of sections into which the packed bed is divided, and so on. Furthermore, the method of calculation and suitable equations for calculating heat transfer and pressure drop are selected.

The code of the created software has thousands of lines and it consists of many classes (for various methods, heat transfer, pressure drop, display results, saving, etc.); therefore, it is not possible to show all the code here. As mentioned in the article, the trapezoidal method is used to solve differential equations. A small part of the code for solving differential Equations (33) and (34) for the hot period is shown in the gray rectangle below.

```
/* -------------------------------------------------------------------------
* start of hot period—packed bed is heated up—linear model
* --------------------------------------------------------------------------
*/
//reverse of temperatures
for (s=0; s < getSection()+1; s++)
{
Tb_H[0 ][s] = Tb_C[getP_C()][getSection()-s];
}
//start—The constants a1, a2, b1, b2, k1, k2 must be calculated before the run this code.
for (s = 0; s < getSection(); s++)
{
Tg_H[0][s+1] = a1h*Tg_H[0][s] + a2h*(Tb_H[0][s+1] + Tb_H[0][s]);
}
for (t = 0; t < getP_H(); t++)
{
Tb_H[t+1][0] = b1h*Tb_H[t][0] + b2h*(Tg_H[t+1][0] + Tg_H[t][0]);
for (s = 0; s < getSection(); s++)
{
Tb_H[t+1][s+1] = k1h*Tb_H[t][s+1] + k2h*Tg_H[t][s+1] + k3h*Tb_H[t+1][s] + k4h*Tg_H[t+1][s];
Tg_H[g+1][s+1] = a1h*Tg_H[t+1][s] + a2h*(Tb_H[t+1][s+1] + Tb_H[t+1][s]);
}
}
for (t=0; t<getP_H(); t++) {
Tgh=Tgh+getTg_H()[t][getSection()]; // Calculate total outlet temp. of gas to get average temp.
Tbh=Tbh+getTb_H()[t][getSection()];
}
Tgh=Tgh/getP_H(); // Calculate average outlet temp. in heating process
Tbh=Tbh/getP_H();
/* -------------------------------------------------------------------------
* end of hot period
* --------------------------------------------------------------------------
*/
```
One possibility of how to improve the created software is to implement an open-source database of fluid properties. Currently, it is necessary to specify fluid properties using interpolation curves, which is not quite effective and delays the calculation. The user must firstly create interpolation curves

and then insert their constants into the software. The next step would be to add the possibility of using the regenerative heat exchanger for flue gas cleaning. This means the software would be able to divide the packed bed into a part where heat accumulation or heat output occurs and a part where the chemical reactions (catalytic bed) take place.

The developed software was used to solve the regenerative heat exchanger. It should be noted the software is still being developed and improved.

#### *Case Study*

The task was to utilize the heat from the gas stream at an inlet temperature of 727 ◦C and to heat up the gas at an inlet temperature of 27 ◦C with the accumulated heat. Flow rates, hot and cold periods, and media properties are given in Table 4. A regenerative heat exchanger with a packed bed diameter of 0.2 m and a height of 1 m should be used. The geometry of the regenerator and properties of the packed bed are given in Table 5. The Willmott linear method was used for calculating the heat exchanger.

**Table 4.** Properties of hot and cold fluid.

It should be noted that the cooling period means that the bed is cooling down and the heating period means that the bed is heating up.

Convergence is attained when the pseudo-thermal ratios given by Equation (37) in two subsequent cycles are numerically equal, with a difference of less than 1 × 10−6. According to the calculation results, the equilibrium was reached after 18 cycles, and the calculation time was 48 milliseconds. It is obvious that the calculation is fast. The main results are shown in Table 6, and graphical dependencies of temperatures and pressures are shown in the figures below. These are only basic graphical outputs provided by the software. Similarly, the output text protocol contains much more data than shown in Table 6.

**Table 6.** Basic results of the calculation.

The basic graphical output is the course of the gas temperatures and the packed bed at the outlet as a function of time (see Figure 7). It can be seen the temperature of the cold medium at the outlet is the highest at the beginning of the period. Gradually, the stored heat decreases and the outlet temperature of the cold gas falls. The opposite situation occurs with hot gas. The hot gas transfers most heat at the start of the period and, gradually, the ability of the bed to absorb energy decreases and the hot gas outlet temperature increases.

**Figure 7.** Dependences of temperatures on the time.

Figure 8 shows the temperature variations of the packed bed material along with its height for the hot and cold periods at steady state. It can be seen that the packed bed has the highest temperature in the first sections. With an increase in the height of the packed bed, the temperature decreases significantly. This is the reason why both media flow gradually in opposite directions.

**Figure 8.** Dependences of temperatures on the location.

Figure 9a,b illustrates that the temperature profiles of the gaseous phase along the bed gradually change with time. In this calculation, they were obtained after a simulation of 18 cycles from the initial condition. Moreover, these figures predict the temperature variations at the same location of the bed but different time of the period.

**Figure 9.** The graphical dependences of gas temperatures along the bed: (**a**) for the cooling period; (**b**) for the heating period.

The pressure drop variations for the hot and cold periods over time is shown in Figure 10. When calculating the pressure drop, the variable physical properties of the gases as a function of temperature were considered. The pressure drop increases with the temperature, while the gas volume increases with temperature and, hence, medium velocity increases. The pressure drop increases during the heating process and decreases during the cooling process with time. Thus, the pressure drop variations correspond to the temperatures shown in Figures 7 and 9.

**Figure 10.** Dependences of pressure drop.

The text and graphical results shown in the paper are not complete results provided with the software. They represent only a selection of the main results.

#### **8. Summary**

Our task was to create software for calculating regenerative heat exchangers with a fixed bed at the request of an engineering office in the Czech Republic. Research revealed there are several computational methods that differ in accuracy and computational stability. It was shown that the open computational methods proposed by Willmott seem to be the most suitable for use on a computer. Although these methods are several years old, they were not suitable because of their complexity at the time of their creation. In subsequent years, the methods were simplified in order to facilitate and accelerate the calculation. This is because, for the exact solution of these types of exchangers, it is necessary to solve the system of differential equations. With the development of computers, however, the situation changed, and these methods are again up to date. The advantage of these methods is their stability and possibility to adjust them by adding suitable calculation formulas for the determination of heat transfer and pressure drop.

The impulse to create this paper was also to broaden the awareness of regenerative heat exchangers, to provide designers with an overview of suitable calculation methods and, thus, to extend the interest and use of these types of heat exchangers. This is the reason why not only calculation methods, but also equations for determining the heat transfer coefficient of convection and radiation and equations for predicting pressure drops are mentioned in the paper. Potential candidates for the calculations of these devices are given as an overview of suitable calculation methods and equations, and one does not need to extensively search for them in the available literature.

**Funding:** This research has been supported by the project No. CZ.02.1.01/0.0/0.0/16\_026/0008413 "Strategic partnership for environmental technologies and energy production", which has been co-funded by the Czech Ministry of Education, Youth and Sports within the EU Operational Programme Research, Development and Education.

**Acknowledgments:** The authors gratefully acknowledge the financial support provided by the EU project Strategic Partnership for Environmental Technologies and Energy Production, funded as project No. CZ.02.1.01/0.0/0.0/16\_026/0008413 by the Czech Republic Operational Program Research, Development, and Education, Priority Axis 1: Strengthening capacity for high-quality research.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Nomenclature**

#### **Substricpts**

#### **References**

© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Online Ash Fouling Prediction for Boiler Heating Surfaces Based on Wavelet Analysis and Support Vector Regression**

**Shuiguang Tong 1,2, Xiang Zhang 1,2, Zheming Tong 1,2,\*, Yanling Wu 3, Ning Tang 1,2 and Wei Zhong 3**

Received: 1 December 2019; Accepted: 16 December 2019; Published: 20 December 2019

**Abstract:** Depending on its operating conditions, traditional soot blowing is activated for a fixed time. However, low-frequency soot blowing can cause heat transfer efficiency to decrease. High-frequency soot blowing not only wastes high-pressure steam, but also abrades surface pipes, reducing the working life of a heat exchange device. Therefore, it is necessary to design an online ash fouling monitoring system to perform soot blowing that is dependent on the status of ash accumulation. This study presents an online monitoring model of ash-layer thermal resistance that reflects the degree of ash fouling. A wavelet threshold denoising algorithm was applied to smooth the thermal resistance of the ash layer calculated by the heat balance mechanism model. Thus, the variation in thermal resistance becomes more visible, which is more conducive to optimizing the operation of soot blowing. The designed Support Vector Regression (SVR) model could achieve the online prediction of thermal resistance denoising for low-temperature superheaters. Experimental analysis indicates that the prediction accuracy was 98.5% during the testing phase. By using the method proposed in this study, online monitoring of heating surfaces during the ash fouling process can be realized without adding complicated and expensive equipment.

**Keywords:** utility boiler; thermal resistance; wavelet threshold denoising algorithm; support vector regression (SVR)

#### **1. Introduction**

In the Chinese economy, the coal-fired boiler continues to play a significant role in converting heat energy into steam or hot water [1]. The energy conservation and environmental protection policies of coal-fired boilers in China have been continuously strengthened in recent years [2]. Many cities are carrying out energy-saving and emission-reduction renovations of boilers. In 2015, China pledged to reduce its CO2 emissions per unit of GDP by 60–65% by 2030 compared with the levels in 2005 [3,4]. Therefore, China's coal-fired boilers have rapidly developed in the direction of high parameters and large capacities in recent years [5]. Supercritical and ultra-supercritical units with rated power generations greater than 600 MW have gradually become the dominant direction of coal-fired units.

Since the end of July 2017, China has put into operation 101 units of 1000 MW ultra-supercritical units. However, with the increase of boiler capacity and parameters, especially the improvement of steam parameters, the problem of ash and slag in boiler furnaces and on convection heating surfaces has become more serious. In the furnace combustion process, ash fouling and slagging on heating surfaces has always been one of the critical factors affecting heat transfer efficiency. After the pulverized coal particles are burned, part of the fly ash passes through the platen superheater, high-temperature superheater, low-temperature superheater, reheater, economizer, air preheater and so on. Ash and slag are deposited on the heating surfaces of all levels. The slope and superheater areas are sometimes covered by massive ash deposits, which in severe cases can cause undesirable unit shutdowns [6]. Due to the small thermal conductivity of the ash, ash accumulation increases the thermal resistance of the heating surface and deteriorates the heat exchange, which in turn causes the exhaust temperature to increase and the boiler efficiency to decrease. When the ash accumulation is serious, ash blockage will occur in the flue, resulting in increased ventilation resistance and reduced boiler output. Sometimes, a power plant is even forced to shut down its furnaces. In addition, high-temperature corrosion and wear of the heating-surface metal tube (caused by ash and slagging) are main factors that can result in boiler tube explosion. Dong et al. [7] introduced a comprehensive experiment in which fly ash particles were impacted under controlled conditions against a flat steel surface to understand the ash deposition process in a pulverized coal boiler system. A systematic study was conducted at a coal-fired power plant, including detailed fuel analysis, boiler and soot-blower characteristics and optimization of slagging in the furnace [8].

The most common solution is to clean the heating surface using various types of soot blowers, as soot blowing can ensure the safety and economy of boiler operation. Shi et al. [9] firstly proposed a cleanliness factor model to monitor the ash deposition status of air preheaters. The analysis of fouling kinetics and optimization of the soot-blowing strategy are then performed to optimize steam consumption and heat transfer efficiency. Shi et al. [10] also established an optimization model of soot blowing on a boiler's economizers. The measurement data and basic thermodynamic calculation data of the distributed control system (DCS) of a thermal power plant are used to calculate the fouling rate of the heated surface in real time. The traditional soot blowing method is purged for a fixed time, depending on operating conditions [11]. The low-frequency soot blowing process can cause heat transfer efficiency to decrease, diminishing a boiler's performance. High-frequency soot blowing, on the other hand, can not only result in the waste of high-pressure steam, but also abrade surface pipes, causing enormous thermal stress and reducing the working life of the heat exchange device while increasing the maintenance cost of the soot blower. Therefore, regular soot blowing will not meet the economic and efficiency requirements of soot blowing.

Different methods for automatic soot blowing that are dependent on the severity of ash in the heated area have become a hot research topic. However, the fouling degree of the heating surface includes many factors, including the characteristics of nonlinearity and strong coupling. This is difficult to simulate with a precise physical model [12], so it is particularly important to establish an online monitoring model that uses thermal resistance to characterize the fouling degree of the heating surface. Through accurate online monitoring of the fouling level on a boiler's heating surface, soot blowing optimization can be implemented to improve the economy and efficiency of boiler operation.

There are three primary methods for predicting the fouling level of a heating surface, including ash fouling monitoring based on instrumental measurements (direct method), ash fouling monitoring based on a mechanism model (indirect method) and ash fouling monitoring based on a data-driven model (indirect method).

In boiler furnaces, instrumentation-based ash fouling monitoring systems include an acoustic pyrometry system. Multi-path measurement of the temperature field in the boiler furnace is carried out by an acoustic generator and an acoustic receiver. The geometric pixel segmentation and regularization algorithm are used to reconstruct the two-dimensional temperature field [13]. Acoustic pyrometry technology is used to monitor the temperature change of flue gas on the local heating surface in the boiler furnace, and a new cleaning factor is defined to characterize the shift in the ash fouling degree [14]. However, acoustic pyrometry technology has its limitations: it can only measure the temperature change of the local heating surface; most of the heating surface is unmeasurable; and additional acoustic sensors are required.

An ash fouling monitoring system based on a mechanism model mainly uses the heat transfer principle to establish the heat transfer model of heating surfaces, and it selects the cleaning factor to reflect the fouling degree of the heating surface [15]. The boiler-side heat absorption models are established to calculate the heat absorption rate of the working fluid in the boiler, including heat exchanger, metal wall and heat loss [16]. The superheater model is used to analyze the effect of the ash layer on the outer tube surfaces and the scale-layer deposit on the inner tube surfaces [17]. Digital signal processing tools are also widely used in the data processing and analysis of boiler monitoring to optimize the ash fouling monitoring system. A dual-extended Kalman filter has been used to estimate the influencing factors of ash deposition, and a cleaning factor indicator has been applied to reflect the fouling degree of heating surfaces [18]. Sobota introduced a method for determining the heat-flow parameters of a steam boiler [19] that can perform online calculations of heating flow rates absorbed by boiler furnaces and superheaters. These parameters can determine the degree of slagging in the furnace. Monitoring the performance of a unit and determining the degree of slagging in the furnace is of great importance. Pronobis studied the influence of biomass co-combustion on the fouling of boiler convection surfaces [20], determining that the co-combustion of straw can cause severe chloride corrosion in superheater tubes. Meanwhile, adopting properly organized soot-blowing technology can eliminate the more serious problem of ash accumulation during the biomass fuel combustion process. Feng et al. [21] built a model of multi-pressure heat recovery steam generator (HRSG) based on the laws of thermodynamics and incorporated energy balance equations for heat exchangers while analyzing flue gas and water/steam parameters such as temperature, pressure, steam mass flow rate and the heat efficiency of different heat exchangers. Tong et al. [22] developed a model to solve the heat transfer calculation of multi-sectional regenerative air heaters (which usually have several layers), and they designed a computational procedure through which to carry out the model. Feng et al. [23] analyzed the steam mass flow rate and outlet temperature of HRSG based on several parameters. This analysis was based on the laws of thermodynamics, and incorporated into the energy balance equations for heat exchangers.

An ash fouling monitoring system based on a data-driven model uses machine learning algorithms to identify and map the relationships between data, avoiding the establishment of complex nonlinear coupled physical models. Enrique et al. (2005) introduced an artificial neural network (ANN) into the prediction of a fouling degree on the boiler heating surface [24]. A probabilistic prediction model for soot blowing based on an artificial neural network and an adaptive neuro-fuzzy inference system was proposed that could avoid the establishment of the nonlinear complex theoretical models involved in fouling dynamics [25]. Temperature differences between the two sides of the tube and shell and the efficiency of the heat exchanger have been predicted based on a local linear wavelet neural network model [26]. The particle swarm optimization (PSO) algorithm can be applied to the hyper-parameter optimization of the support vector machine (SVM) to establish a fouling prediction model for the heat exchanger, which provides another method for the fouling prediction of the heat exchanger [27]. However, the robustness of the pure data-driven model is weak, and easily deviates from the actual value under small disturbances.

To enhance the stability and prediction accuracy of ash fouling monitoring systems, this study proposes a method of combining the mechanism model and the data-driven model for the fouling prediction of boiler heating surfaces.

#### **2. Online Ash Fouling Monitoring Model of a Low-Temperature Superheater**

Figure 1 shows the overall framework of the training process for the ash accumulation monitoring model. Firstly, according to the principle of heat balance, a heat transfer mechanism model is established to calculate the thermal resistance of the ash layer, after which the wavelet threshold denoising algorithm is used to filter out the noise fluctuation in the thermal resistance. Finally, this paper establishes an online fouling prediction model for heating surfaces based on support vector regression (SVR) [28]. This method can filter jagged fluctuations in the thermal resistance data and

predict the ash fouling severity of low-temperature superheaters online, providing more-favorable conditions for the research of further soot blowing optimization strategies.

#### *2.1. The Thermal Resistance Model of a Low-Temperature Superheater*

As shown in Figure 2, heat balance and the heat transfer calculation of heating surfaces are based on the basic theory of boiler thermal balance and performed one by one along opposite directions of the flue gas flow from the outlet of the air preheater. The inlet working-fluid temperature of the low-temperature superheater is calculated according to the heat balance principle—Equations (1) and (2) for the working fluid side and the flue gas side based on the known inlet and outlet gas temperatures of each heating surface, the working fluid outlet parameters and the pipe arrangement structure of each heating surface. Heat transfer Equation (7) can then be applied to calculate the heat transfer coefficient of the boiler's low-temperature superheater.

**Figure 2.** The flue gas flow process on the heating surfaces of the boiler.

Heat absorption on working fluid side:

$$Q\_{\oplus^x} = \frac{D(i' - i' + \Delta i\_{jw})}{B\_j} \tag{1}$$

Heat release on flue gas side:

$$Q\_{yq} = \varphi \left( l' - l'' + \Delta \alpha \cdot l\_{lk}^0 \right) \tag{2}$$

where *Qgz* is the convection heat absorption on the working fluid side; *Qyq* is the convection heat release on the flue gas side; *D* is the quantity of the working fluid flow; *Bj* denotes the calculating fuel quantity; *i* and *i* represent the enthalpy of the inlet working fluid and outlet working fluid, respectively; *I* and *I* are the enthalpy of the inlet flue gas and outlet flue gas, respectively; Δ*ijw* is the reduced value of steam enthalpy in the desuperheater; ϕ is the heat retention coefficient; Δα is the air leakage coefficient; and *I* 0 *lk* is the theoretical cold air enthalpy.

Physical parameters such as *D*, *Bj*, *I* , *i* , *I* , Δ*ijw* and so on are collected in the distributed control system (DCS), and ϕ, *I* 0 *lk*, Δα and so on can be obtained through manual boiler design and thermal calculation. According to the heat balance equation *Qgz* = *Qyq*, the enthalpy of the inlet working fluid at the low-temperature superheater is as follows:

$$\dot{\mathbf{r}}' = \dot{\mathbf{r}}'' + \Delta \dot{\mathbf{r}}\_{jw} - \frac{qB\_j(I' - I'' + \Delta \alpha I\_{lk}^0)}{D} \tag{3}$$

According to the enthalpy temperature table of superheated steam, as shown in Figure 3, the corresponding inlet working fluid temperature *T* is expressed as follows:

$$T = f(I, P) \tag{4}$$

where *I* is superheated steam enthalpy, and *P* is superheated steam pressure.

**Figure 3.** The enthalpy diagram of superheated steam.

Next, the logarithmic mean temperature difference Δ*t* is obtained from the inlet and outlet working fluid temperatures and the inlet and outlet flue gas temperatures, which can be described as follows:

$$
\Delta t = \frac{\Delta t\_d - \Delta t\_x}{\ln \frac{\Delta t\_d}{\Delta t\_x}} \tag{5}
$$

where Δ*td* is the difference between the inlet temperature of the flue gas side and the inlet temperature of the working fluid side of the low-temperature superheater, and Δ*tx* is the difference between the outlet temperature of the flue gas side and the outlet temperature of the working fluid side of the low-temperature superheater.

When the maximum temperature difference Δ*td* and the minimum temperature difference Δ*tx* are satisfied using Δ*td* Δ*tx* ≤ 1.7, the logarithmic mean temperature difference Δ*t* can be simplified as follows:

$$
\Delta t = \frac{1}{2} (\Delta t\_d + \Delta t\_\mathbf{x}) \tag{6}
$$

Finally, the actual heat transfer coefficient of the heating surface should be

$$K = \frac{Q\_{\mathcal{S}^2} B\_j}{\Delta t \cdot H} \tag{7}$$

where *H* is the heat transfer area of the low-temperature superheater.

Taking the low-temperature superheater as an example, Figure 4 depicts the heat transfer in a single superheater tube. The tubular convective heating surface is regarded as a heat transfer model of multi-layer cylindrical wall, and the heat transfer coefficient *K* is calculated as follows:

$$K = \frac{1}{\frac{1}{a\_1} + \frac{1}{2\pi L \lambda\_\xi} \ln \frac{r\_2}{r\_1} + \frac{1}{2\pi L \lambda\_m} \ln \frac{r\_2}{r\_2} + \frac{1}{2\pi L \lambda\_h} \ln \frac{r\_4}{r\_3} + \frac{1}{a\_2}} = \frac{1}{\frac{1}{a\_1} + R\_\xi + R\_m + R\_h + \frac{1}{a\_2}}\tag{8}$$

where α1 is the heat release coefficient of the flue gas side; α2 is the heat absorption coefficient for the working fluid side; λ*h* is the ash-layer thermal conductivity; λ*m* is the metal tube thermal conductivity; λ*g* is the scale-layer thermal conductivity; *r*1 is the inner radius of the scale layer; *r*2 is the inner radius of the metal tube; *r*3 is the outer radius of the metal tube; *r*4 is the outer radius of the ash layer; *L* is the length of the metal tube; *Rh* is the thermal resistance of the ash layer; *Rm* is the thermal resistance of the metal tube; and *Rg* is the thermal resistance of the scale layer.

**Figure 4.** Heat transfer of the low-temperature superheater tube wall.

Since the influence of *Rm* is small, the effect of *Rm* is ignored. Before the raw water is replenished into the boiler, the power plant must treat the boiler feed water to remove the salts, impurities and gases, so that the quality of the supply water meets certain requirements. Therefore, the thermal resistance of the scale layer *Rg* is small and can be ignored to simplify the calculation of *R*m = *Rg* = 0. Thus, the thermal resistance of the ash layer *Rh* is calculated as follows, according to the Equation (8):

$$R\_h = \frac{1}{K} - \frac{1}{\alpha\_1} - \frac{1}{\alpha\_2} \tag{9}$$

where the heat release coefficient of the flue gas side α1 is the convective heat transfer coefficient α*d* and the radiation heat release coefficient α*f* , and is calculated as:

$$
\alpha\_1 = \alpha\_d + \alpha\_f \tag{10}
$$

For the low-temperature convection heating surface, the radiation heat release coefficient α*f* = 0. When tube bundles on the convection heating surface are arranged in parallel, the heat release coefficient on the flue gas side becomes

$$
\alpha\_1 = \alpha\_d = 0.2 \text{C}\_s \text{C}\_z \frac{\lambda\_y}{d} \left(\frac{\omega\_y d}{\nu\_y}\right)^{0.65} \left(\frac{\mu\_y \text{C}\_{p,y}}{\lambda\_y}\right)^{0.33} \tag{11}
$$

When the tube bundles on the convection heating surface are arranged in staggered rows, the convective heat transfer coefficient on the flue gas side is calculated as

$$\alpha\_d = \mathbb{C}\_s \mathbb{C}\_z \frac{\lambda\_y}{d} (\frac{\omega\_y d}{\nu\_y})^{0.6} (\frac{\mu\_y \mathbb{C}\_{p,y}}{\lambda\_y})^{0.33} \tag{12}$$

The heat absorption coefficient α2 on the working fluid side is given by

$$\alpha\_2 = 0.023 \text{C}\_t \text{C}\_l \frac{\lambda\_f}{d\_{dl}} (\frac{\alpha\_f d\_{dl}}{\nu\_f})^{0.8} (\frac{\mu\_f \text{C}\_{p,f}}{\lambda\_f})^{0.4} \tag{13}$$

where *Cs*,*Cz*,*Ct*,*Cl* are the correction factors determined by the structural dimensions of tube bundles, and represent the correction factor related to pipe pitch, the correction factor associated with the vertical tube row number, the correction factor related to airflow and wall temperature and the relative length correction factor, respectively; *d* is the outer diameter of the low-temperature superheater tube; *ddl* is the inner diameter of the low-temperature superheater tube; ω*f* and ω*y* are the flow rate of the working fluid and flue gas at the average temperature, respectively; λ*f* and λ*y* are the thermal conductivity of the working fluid and flue gas at the average temperature, respectively; ν*f* and ν*y* are the kinematic viscosity of the working fluid and flue gas at the average temperature, respectively; μ*f* and μ*y* are the dynamic viscosity of the working fluid and flue gas at the average temperature, respectively; and *Cp*, *f* and *Cp*,*y* are the constant-pressure specific heat of the working fluid and flue gas at the average temperature.

The thermal resistance of the ash layer *Rh* is used to characterize the fouling degree of the convection heating surface. Generally, the larger the value is, the more serious the fouling of the heating surface is.

#### *2.2. Wavelet Decomposition Model*

In 1988, Mallat proposed a concept of multi-resolution with a fast algorithm for wavelet decomposition and reconstruction [29,30]. The original signal *f*(*t*) ∈ *L*2(*R*) Ψa,*b*(*t*) is a continuous wavelet. Then, the inner product of *f*(*t*) and Ψa,*b*(*t*), which is called the continuous wavelet transform, is described as

$$\mathcal{W}\_f(a,b) = \left\langle f(t), \Psi\_{\mathbf{a},b}(t) \right\rangle = \frac{1}{\sqrt{a}} \int\_{\mathbb{R}} f(t) \* \Psi(\frac{t-b}{a}) dt \tag{14}$$

However, in actual engineering calculations, the signals are discrete. As a result, it is necessary to discretize the wavelet transform. The binary discrete wavelet transform (DWT) can be expressed as

$$\mathcal{W}\_f(2^j, b) = \left\langle f(t), \Psi\_{j,b}(t) \right\rangle = 2^{-j/2} \int\_{\mathbb{R}} f(t) \ast \Psi(\frac{t-b}{2^j}) dt \tag{15}$$

The basic principle is as follows. A signal *f*(*t*) in space *Vj* can be represented by basic functions in two orthogonal subspaces, *Vj*+1 and *Wj*+1, which can be determined using Equation (16). According to Equation (17), the first level decomposes *A*0 into a low-frequency part *A*1 and a high-frequency part *D*1, the second level decomposes *A*1 into a low-frequency part *A*2 and a high-frequency part *D*2 and so on, until the multi-resolution decomposition of signals can be realized.

$$V\_j = V\_{j+1} \oplus \mathcal{W}\_{j+1} \tag{16}$$

$$\begin{array}{l} f(t) = \sum c A\_0(k) \boldsymbol{\uprho}\_{j,k}(t) \\ = \sum c A\_1(k) \boldsymbol{\uprho}\_{j+1,k}(t) + \sum c D\_1(k) \boldsymbol{\uprho}\_{j+1,k}(t) \\ = \sum c A\_2(k) \boldsymbol{\uprho}\_{j+2,k}(t) + \sum c D\_2(k) \boldsymbol{\uprho}\_{j+2,k}(t) + \sum c D\_1(k) \boldsymbol{\uprho}\_{j+1,k}(t) \\ = \dots \end{array} \tag{17}$$

In the scale metric space *Vj*, the coefficient *A*0(*k*) is decomposed to two wavelet coefficients *A*1(*k*) and *D*1(*k*) in the scale metric spaces *Vj*+1 and *Wj*+1. Similarly, the two wavelet coefficients *A*1(*k*) and *D*1(*k*) can be used for reconstruction to get *A*0(*k*). The reconstruction algorithm and the decomposition algorithm are corresponding and mutually inverse.

#### *2.3. VisuShrink Soft Threshold Denoising*

Generally speaking, high-frequency signals contain noise details. The wavelet threshold denoising algorithm is an effective filtering method. After wavelet decomposition, the threshold wavelet method is used to weight the decomposed wavelet coefficients. For high-frequency wavelet coefficients, the VisuShrink soft threshold denoising method [31] is adopted for signal processing. All low-frequency wavelet coefficients representing the global trend of the original signal are retained. Then, all small signals are reconstructed to obtain denoising signals, as shown in Figure 5.

**Figure 5.** Wavelet threshold denoising process.

*Energies* **2020**, *13*, 59

The VisuShrink method selects the high-frequency wavelet coefficients to estimate the standard deviation of the noise:

$$\sigma = \frac{median\left(\left|\omega\_{\bar{j}}(k)\right|\right)}{0.6745} \tag{18}$$

The global threshold λ is determined as

$$
\lambda = \sigma \sqrt{2 \log n} \tag{19}
$$

Finally, the soft threshold denoising equation is written as follows

$$\omega\_{\lambda} = \begin{cases} \begin{bmatrix} \operatorname{sgn}(\omega)(|\omega| - \lambda) \end{bmatrix} & |\omega| \ge \lambda \\ 0 & |\omega| < \lambda \end{cases} \tag{20}$$

#### *2.4. SVR Theory*

An SVR [32–34] method is usually realized through regression and prediction. The principle of SVR is to learn a function *f*(*x*), so that the function value is as close as possible to the real value. Given the training samples *D* = ) (*x*1, *y*1),(*x*2, *y*2),(*xm*, *ym*) \* , with *xi* ∈ *Rn* for the input, *yi* ∈ *R* for the target output, *m* as the number of training samples. The regression model equation is calculated as

$$f(\mathbf{x}) = w\phi(\mathbf{x}) + b \tag{21}$$

where *w*, *x* ∈ *Rn*; b ∈ *R*, φ(*x*) is kernel function. SVR transforms low-dimensional nonlinear problems into high-dimensional linear problems by introducing kernel functions. Solving the dual Lagrangian problem of SVR:

$$\max\_{\alpha, \hat{\alpha}} \sum\_{i=1}^{m} y\_i (\hat{\alpha}\_i - \alpha\_i) - \varepsilon (\hat{\alpha}\_i + \alpha\_i) - \frac{1}{2} \sum\_{i=1}^{m} \sum\_{j=1}^{m} (\hat{\alpha}\_i - \alpha\_i)(\hat{\alpha}\_j - \alpha\_j) \mathcal{K}(\mathbf{x}\_i, \mathbf{x}\_j) \tag{22}$$

This is subject to

$$\sum\_{i=1}^{m} \left(\vec{\alpha}\_{i} - \alpha\_{i}\right) = 0\tag{23}$$

$$0 \le \alpha\_{i\prime} \land\_i \le \mathbb{C}, i = 1, 2, \dots, m \tag{24}$$

The above process must satisfy Karush–Kuhn–Tucker (KKT) conditions:

$$\begin{cases} \alpha\_i (f(\mathbf{x}\_i) - y\_i - \varepsilon - \underline{\xi}\_i) = 0 \\ \quad \pounds\_i (y\_i - f(\mathbf{x}\_i) - \varepsilon - \underline{\hat{\xi}}\_i) = 0 \\ \quad \alpha\_i \pounds\_i = 0, \underline{\xi}\_i \underline{\xi}\_i = 0 \\ (\mathcal{C} - \alpha\_i) \underline{\xi}\_i = 0, (\mathcal{C} - \alpha\_i) \underline{\xi}\_i = 0 \end{cases} \tag{25}$$

Eventually, the solution to SVR is

$$f(\mathbf{x}) = \sum\_{i=1}^{m} (\pounds\_i - \alpha\_i) \mathsf{K}(\mathbf{x}\_i, \mathbf{x}) + b \tag{26}$$

$$b = y\_i + \varepsilon - \sum\_{i=1}^{m} (\alpha\_i - \alpha\_i) \mathcal{K}(\mathbf{x}\_i, \mathbf{x}) \tag{27}$$

where *K*(*xi*, *x*) = φ(*xi*) *T*φ(*x*) is the kernel function.

The kernel function can map the nonlinear problem of low-dimensional space to high-dimensional space, which will then become a linear problem. However, constructing the kernel function *K*(*xi*, *x*) is a significant problem. The most crucial step is to determine the mapping of input space to feature space, which can only be achieved if we know the distribution of data within the input space. However, in most cases, we do not know the specific distribution of data being processed. Therefore, it is generally challenging to construct a kernel function that conforms precisely to the input space, and as a result this paper uses the radial basis function (RBF) instead of rebuilding a new kernel function.

$$K(\mathbf{x}, \mathbf{x}\_{\mathrm{i}}) = \exp(-\frac{\left\|\mathbf{x} - \mathbf{x}\_{\mathrm{i}}\right\|^2}{\delta^2})\tag{28}$$

The Gaussian radial basis function is one of the most widely used kernel functions. It offers a better performance regardless of whether a sample is large or small, and it has fewer parameters than the polynomial kernel function. Therefore, the Gaussian kernel function is preferred in most cases.

Vapnike's [35] research demonstrates that the kernel parameter δ and penalty factor *C* are the key factors affecting the performance of SVM. A larger *C* will reduce the training error, but will also result in over-fitting at the same time, which will increase the test error. When the kernel parameter δ is small, the regression prediction has better accuracy. However, if the kernel parameter δ is much lower, the accuracy of the model will drop significantly.

#### **3. Case Study and Data Collection**

This paper takes a 170 t/h tangentially fired pulverized coal boiler from the Huilian Power Plant in Wuxi, Jiangsu Province, China, as the research object. The horizontal section size of the furnace chamber is 7090 × 7090 mm (width × depth). Table 1 gives the design parameters of the boiler.

**Table 1.** Design parameters of the boiler (rated load).

This type of case study boiler has a typical heating surface structure similar to that of many power station boilers in China. Figure 6 shows the schematic diagram of the case study boiler system. It can be seen that after the pulverized coal is burned by the furnace to produce flue gas, the flue gas flows along the flue gas passage all the way through the platen superheater, the high-temperature superheater, high-temperature reheater, low-temperature superheater, low-temperature reheater, economizer, air preheater and many other heating surfaces to the desulfurization tower and bag filter before finally being discharged by the chimney. The heat transfer mode of the heating surfaces is convection, except for the platen superheater which is half radiation and half convection. In the actual operation of the boiler, the flue gas contains a large number of fly ash particles. When fly ash is deposited on the heated surface, the heat transfer performance is degraded.

The data used in this study was collected from the DCS of power plants, as shown in Table 2. By analyzing the formation mechanism and influence factors of ash fouling, 20 characteristic parameters affecting the ash fouling were selected, including the eight powder feeder speeds, the main steam flow rate, the main steam pressure, the main steam temperature, the outlet steam temperature, the

inlet flue gas temperature, the outlet flue gas temperature, the steam flow rate, the feedwater flow rate, the feedwater temperature, the boiler oxygen amount, the air supply rate of blower A and the air supply rate of blower B. The function of the powder feeder is to continuously and uniformly send pulverized coal from the pulverized coal bin to the furnace for combustion according to the coal burning amount required by the boiler load. The speed can control the amount of pulverized coal entering the furnace. Blower A is a primary fan of the boiler that is used to dry the fuel and send it into the furnace. Blower B is a secondary fan of the boiler that is used to overcome resistance from the air preheater, air duct and burner, and to maintain full combustion of fuel. The air blowers are located at the air inlet of the air preheater (Figure 6, number 11). In summary, this paper uses 20 characteristic parameters that formed the input set of the SVR model. The input data were collected every three minutes from 25 May 2019 at 00:00 to 29 May 2019 at 23:57. After the normalization of these 2400 sets of data, all samples were divided into an 80% training set and a 20% test set.

**Figure 6.** Schematic diagram of the boiler system case study. 1—water wall, 2—steam-water separator, 3—steam turbine, 4—platen superheater, 5—high-temperature superheater, 6—high-temperature reheater, 7—low-temperature superheater, 8—low-temperature reheater, 9—economizer, 10—air preheater, 11—air blower, 12—desulfurization tower, 13—bag filter.

In this study, the thermal resistance of the ash layer was chosen to be an indicator by which to detect the degree of fouling on the heated surface. After the original thermal resistance data were de-noised by wavelet analysis, the focus was to set up correlations between the 20 characteristic parameters and the thermal resistance of the ash layer using the SVR model.

 Results and Discussion**

#### *4.1. Analysis of the Ash-Layer Thermal Resistance*

Taking advantage of the thermal resistance mechanism model for the low-temperature superheater established in Section 2.1, the thermal resistance of the ash layer *Rh* combined with the operation data from the Huilian power plant was calculated. Figure 7 is the calculation result of the ash-layer thermal resistance for the low-temperature superheater. The smaller the thermal resistance of ash-layer value, the better the heat transfer performance of the low-temperature superheater was.

**Figure 7.** The changing trend of the thermal resistance for the low-temperature superheater.

It can be seen that there was a sudden drop every day from 06:00–08:00 as the result of a scheduled ash blowing operation by the power plant, and that the thermal resistance of the ash layer of the heating surface decreased significantly after soot blowing.

On the morning of 27 May 2019, the degree of ash fouling on the heating surface was not serious enough, but the operator still performed the blowing operation, which exposed the shortcomings of the scheduled soot blowing action. Scheduled soot blowing does not consider whether ash accumulation is severe, and wastes steam that causes severe pipe wall wear, resulting in a decrease in boiler efficiency. The improper operation of the soot blower will not only reduce the economic benefit of unit operation, but will also lead to pipe wall wear or even an explosion of the pipe, affecting the safety of the boiler unit.

This paper calculates the thermal resistance of the ash layer in order to observe and visualize degrees of ash fouling on the heating surface. Based on quantification, the original soot blowing method can be changed from scheduled soot blowing to on-demand soot blowing.

On the other hand, it can be observed that there are a large number of sawtooth fluctuations in Figure 7. These were caused by noise data and unstable working conditions during operations. Changes in working conditions can affect the heat transfer performance of the heating surface. These fluctuations are very unfavorable for both the operator's observations and the development of an ash blowing strategy.

#### *4.2. Wavelet Threshold Denoising Results*

The wavelet threshold denoising algorithm was performed to smooth and filter the thermal resistance data, reducing the influence of sawtooth fluctuations on the degree of ash fouling. This method not only preserved the changing trend of thermal resistance, but also lessened the sawtooth jitter. It also provided cleaner raw data for the next training of the SVR model.

In this experiment, the sym5 wavelet was selected with two to five decomposition levels. Figure 8 shows the waveform diagram for each high-frequency and low-frequency signal after wavelet decomposition. The comparison results for the two to five wavelet decomposition levels and the denoising are shown in Figure 9a–d.

The signal-to-noise ratio (SNR) and root mean square error (RMSE) are shown after denoising in Table 3, which demonstrates the accurate quantization of the denoising effect. The higher the SNR of the signal, the smaller the RMSE, and the closer the denoising signal was to the original signal. It can be seen that the denoised curve not only retained the changing trend of the original signal, but also eliminated the influence of partial noises.

$$SNR = 10 \log \left[ \frac{\sum\_{n} \mathbf{x}^2(n)}{\sum\_{n} \left[ \mathbf{x}(n) - \mathbf{x}'(n) \right]^2} \right] \tag{29}$$

$$RMSE = \sqrt{\frac{1}{n} \sum\_{n} \left[ \mathbf{x}(n) - \mathbf{x}'(n) \right]^2} \tag{30}$$

It can be seen in Figure 9 and Table 3 that after calculation with the four-level wavelet threshold denoising algorithm, the results of SNR and RMSE were 31.8 and 0.000865, respectively. This indicates that thermal resistance data after a four-level wavelet threshold denoising filtering not only maintained the overall trend, but also reduced sawtooth fluctuations to a minimum. Thus, the wavelet threshold denoising algorithm was rendered valid by filtering out some noise information.

**Figure 8.** Results of four-level wavelet decomposition using the sym5 wavelet.

**Figure 9.** Wavelet threshold denoising results.

**Table 3.** Root mean square error (RMSE) and signal-to-noise ratio (SNR) of wavelet denoising signals with different levels.

However, the five-level wavelet threshold denoising lost some of the details of the original signal, thus the four-level wavelet threshold denoising results were applied as the output data sets of the SVR model described in Section 2.4. After wavelet decomposition and denoising, the changing trend of thermal resistance was more visible, which is more conducive to the design of blowing ash optimization strategy.

#### *4.3. Fouling Prediction Based on SVR*

The SVR model described in Section 2.4 was utilized for training the mapping relationship between characteristic parameters and the thermal resistance of the ash layer. The grid search method was adopted to optimize the two hyper-parameters of the SVR, including penalty factor *C* and kernel parameter δ.

Predictive simulation experiments were performed in the Python 3.6 environment. The first 80% of samples were divided into a training set and the remaining 20% into a testing set. The RBF kernel function was then presented. The performance evaluation criterion selected for the SVR model was the coefficient of determination *R*2 (optimum was the maximum value). The grid search method was used for hyper-parametric optimization to improve the prediction accuracy of the SVR model. Figure 10 shows the model prediction accuracy under different hyper-parameters (*C*, δ). Training and testing accuracy of the SVR model using four-level wavelet decomposition and denoising data are provided in Table 4. The optimum hyper-parameters of the SVR model were (29, 0.13), and the accuracy of the training set and testing set were 0.994 and 0.985, respectively.

**Figure 10.** Prediction accuracy with different hyper-parameters (*C*, δ) based on the grid search algorithm.

**Table 4.** Training and testing accuracy of the SVR model by using the four-level wavelet decomposition denoising signal.

The output results for the thermal resistance of the ash layer were predicted by simulating the SVR model and comparing it with the actual results of the mechanism model. The original data and forecasted data of the thermal resistance for the training and testing phases are illustrated in Figure 11. Figure 12 is a graph of the absolute error residual. It can be observed that the predicted value of the SVR model was within a range of ±5% of the actual value of the mechanism model, as shown in Figure 13.

**Figure 11.** Comparison between the original value and the predicted value of thermal resistance.

**Figure 13.** Comparison of predicted and actual thermal resistance.

The SVR model can predict the fouling process of the heating surface with precision by using 20 characteristic parameters and avoiding the establishment of complex equations such as heat transfer and radiation. The ash fouling prediction model of boiler heating surfaces based on wavelet analysis and SVR has good generalization and can meet the actual soot blowing operation requirements.

#### **5. Conclusions**

In this study, a thermal resistance network of the ash layer is developed to characterize the degree of boiler ash fouling for the first time, and a physical model was established based on the heat transfer balance principle of the low-temperature superheater. To reduce the adverse effects of sawtooth fluctuations, a method based on wavelet analysis and SVR was proposed to improve the online prediction of thermal resistance. The case study showed that the SNR and RMSE obtained using a four-level wavelet threshold denoising algorithm to filter the thermal resistance data were 31.8 and 0.000865, respectively. This indicates that the filtered thermal resistance data not only maintained the overall trend, but also reduced the sawtooth fluctuation to a minimum. With the filtered thermal resistance data, an SVR model was trained to map the relationship between characteristic parameters and the thermal resistance of the ash layer. The maximum prediction accuracy was 99.4% during the training phase and 98.5% during the testing phase. The accuracy of this method can meet real engineering needs and provide an innovative way for the analysis and prediction of ash accumulation for a low-temperature superheater.

The wavelet analysis and SVR model make it possible to predict ash fouling online without significant sawtooth fluctuation. Therefore, considering the degree of agreement between the predicted results and the actual calculated results, the online prediction method of thermal resistance based on wavelet analysis and SVR can be used as an appropriate prediction tool for estimating and predicting the ash fouling behavior of a low-temperature superheater without adding complicated and expensive equipment.

**Author Contributions:** Conceptualization, S.T. and X.Z.; data curation, Y.W.; investigation, N.T.; methodology, X.Z and W.Z.; software, X.Z; supervision, S.T. and Z.T.; validation, X.Z., Z.T. and N.T.; visualization, Y.W.; writing—original draft, S.T., X.Z., Z.T. and W.Z.; writing—review and editing, S.T., X.Z., Z.T. and W.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Key R & D Program of China (Grant No. 2017YFA0700305 and 2018YFB0606105), Zhejiang Provincial Natural Science Foundation (Grant No. LR19E050002, LY17F030007), and Zhejiang University & Huaguang Smart Energy System Joint Research Center.

**Acknowledgments:** We are very grateful for the contributions of Jianyun Zhao from Hangzhou Boiler Group Co., Ltd., Hangzhou, China, Junhua Mao from Wuxi Huaguang Boiler Co., Ltd., Wuxi, China, Liangwei Xia from Harbin Boiler Co., Ltd., Harbin, China, and Wen Yang from Jianglian Heavy Industry Group Co., Ltd., Nanchang, China. These industrial collaborators mentioned above made significant contributions to this paper although they were not co-authored.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **E**ff**ect of Saturated Steam Carried Downward on the Flow Properties in the Downcomer of Steam Generator**

#### **Yan Liu 1,2, Hui-lie Shi 2, Chun Gui 2, Xian-yuan Wang 1,2,\* and Rui-feng Tian 2**

Received: 2 September 2019; Accepted: 20 September 2019; Published: 24 September 2019

**Abstract:** The saturated water separated by the steam separator in a natural circulation steam generator may carry a small amount of saturated steam into the downcomer. The steam contacts subcooled water and condenses directly in the downcomer causing the variations in the pressure and steam quality and likely affecting the stability of the water cycle in the secondary loop. It is not conducive to core heat extraction and thus affects nuclear safety. The mathematical model of the downcomer was established in this study based on the internal structures of a natural circulation steam generator. The volume-of-fluid (VOF) and large-eddy-simulation (LES) models were used for analysis on FLUENT software (ANSYS, Pittsburgh, PA. USA) platform. The influence of direct contact condensation of top-down flowing steam on the flow properties in the downcomer of the steam generator under high pressure was studied. The trend of the temperature, pressure, and the void fraction were obtained by combining these models with the condensation model. Further, a one-dimensional calculation program based on the differential drop was also developed to assess the flow field in the downcomer. The calculation results are in good agreement with the experimental results which indicated that, when affected by the saturated steam carried downward, the flow temperature close to the exit of the downcomer rises slightly due to the absorption of the heat released by the steam condensation. Furthermore, the density corrected by the pressure-drop is more reliable than that corrected by the temperature. After the velocity in the downcomer has increased to a certain value, the sensitivity of steam quality to the subcooling degree in the downcomer begins to decline. The results in this paper can be used to perform stability analyses and to design steam generators. The results of research are helpful to the stability analysis and the design of a steam generator, and to improve the accuracy of the measurement of the steam generator operating parameters, thus enhancing the safety of Pressurize Water Reactor operating system.

**Keywords:** direct contact and condensation of steam; natural circulation steam generator; downcomer; condensation model

#### **1. Introduction**

A natural circulation steam generator is a critical part of the nuclear power plant equipment [1]. It is mainly used to absorb the heat of reactor coolant as well as to generate the saturated steam that drives a steam turbine. The basic circulation loop of a natural circulation steam generator is composed of an ascending channel, a steam separator, and a downcomer [2]. As shown in Figure 1, the feed water flows from the feed pipe into the steam generator and from the top to bottom along the downcomer. The feed water flows over the support plate into the ascending channel and absorbs heat in the ascending channel to generate saturated steam. This saturated steam is then filtered out by the

steam separator, while the water filtered by the steam separator, known as recycled water, carries part of the steam, mixes it with the feed water, and flows along the downcomer. Since the density of the water in the downcomer is greater than the steam mixture in the ascending channel, the water flows down the downcomer, whereas the mixture flows up the ascending channel.

**Figure 1.** Schematic diagram of the steam generator simulator.

The following relationship can be assumed under the condition that the water circulation is stable [2]:

$$\rm{P} + \rho\_{\rm{m}} \rm{g} \rm{H} - \Delta \rm{P}\_{\rm{dc}} = \rm{P} + \rho\_{\rm{x}} \rm{g} \rm{H}\_{\rm{x}} + \rho\_{\rm{x} \rm{s}} \rm{g} \rm{H}\_{\rm{xs}} - \Delta \rm{P}\_{\rm{r}} - \Delta \rm{P}\_{\rm{sp}} \tag{1}$$

where P is the steam pressure, ρm is the density of the downcomer, H is the height of the downcomer, ρx is the density of the pre-heating section in the ascending channel, Hx is the height of the pre-heating section, ρxs is the density of the vapor section in the ascending channel, Hxs is the height of the vapor section, ΔPr is the ascending channel resistance, and ΔPsp is the steam separator resistance.

The left side of Equation (1) represents the driving force for the circulation of the steam generator's secondary loop, while the right side is its resistance. The saturated steam directly contacts the downcomer walls and condenses. This changes the density and resistance in the downcomer and might lead to the decrease of the density in the downcomer and the decrease of the driving force, thus affecting the stability of the water cycle in the secondary loop. Therefore, the study of this topic has certain research significance.

Woods et al. conducted six groups of different experiments based on the advanced plant experimental test facility (APEX) for AP600 in the Oregon State University [3]. These experiments were aimed to evaluate the steam condensation rate of an established scaled model [3,4] of steam generators under different pressures and inlet flow rates of primary and secondary loops. Brucker and Sparrow researched the direct contact and condensation of steam bubbles under high-pressure conditions [5]. The process was simulated from below through the experiments at 1-6 MPa. Their study illustrated that the time for steam bubbles to collapse increases with increasing pressure but decreases when the subcooling degree rises.

Numerical simulations have been carried out as well. Yang et al. established a three-dimensional physical model of a steam generator unit tube based on similarity theory using the steam generator of Daya Bay Nuclear Power Plant as a prototype [6]. The numerical simulations of the two-phase flow and boiling heat transfer properties at the secondary loop of the steam generator were performed via the CFX software(ANSYS, Pittsburgh, PA, USA) [7,8] using a particle model and a phase-transformation model [9]. However, the trend of the void fraction on the cross-section of a steam generator ascending channel was also investigated. Indeed, Jiang et al. established a computational model of a three-dimensional flow at the secondary loop of the steam generator via FLUENT [10] software using a porous-medium model under the condition that it is single-phase flow at the secondary loop of the steam generator [11]. The three-dimensional flow field at the secondary loop of a steam generator during the steady-state operation of the nuclear power plant was computed, and the distribution of pressure and velocity for the entire flow field was obtained. Li et al. [12] performed a numerical simulation on the fluid flow in the feed-water pipe and downcomer via the CFX software using the re-normalization group (RNG) k-εturbulence model [13]. Furthermore, the mixing of feed and recirculation water in the steam generator downcomer as well as the distribution of the flow rate and temperature were mainly investigated under a single-phase condition. Zhou et al. [14] investigated the shape variation of the steam column under different steam pressures using the Euler multiphase flow model [15], and the realizable k-ε turbulence [13] and thermal equilibrium change [13] models. Gulawani et al. carried out Computational Fluid Dynamics simulations of the direct contact and condensation of steam using a double-resistance model [16]. In their model, there was no resistance on the vapor-phase side, and the Ranz-Marshall model [17] was adopted in the liquid phase. The computational results were in good agreement with the experimental data. Li et al. conducted numerical simulations and mechanism analyses of the steam direct contact condensation in the t-shaped circular tubes by using the gas-liquid two-phase flow and large eddy simulation turbulence model [18] combined with the double-resistance model on the FLUENT platform. None of the above studies investigated the influence of the steam direct-contact condensation on the internal flow field of the steam generator under the top-down flow state in the downcomer of the steam generator.

This work aims to investigate the condensation state of saturated steam carried to the narrow downcomer of the steam generator in the direction of the top to bottom and its influence on the flow characteristics. The mathematical model of the downcomer was established based on the internal structures of natural circulation steam generator, and the volume-of-fluid (VOF) and large-eddy-simulation (LES) models combined with the condensation model were used for the analysis. A small steam generator simulator and an experiment loop were designed in accordance with a typical power plant system. To this end, the experimental and computational analyses were carried out.

#### **2. Computational Analysis**

According to the internal structure of the steam generator, the equation of mass, momentum, and energy conservation was established for the direct contact condensation process of steam and subcooled water in the downcomer of the steam generator. The phase change models of condensation and mass transfer at the vapor-liquid condensation interface were therefore built. An appropriate condensation heat and mass transfer model was also built and a user define function (UDF) was written and embeded in the Computational Fluid Dynamics software. The changes in the velocity field, temperature field, and pressure field during direct contact condensation were obtained.

#### *2.1. Physical Model*

As shown in Figure 2, the actual physical model was reconstructed in equal proportions. The inlet section of the vapor mixture, the width of the downcomer, and the size of the subcooled water inlets are consistent with the real object. The inlet diameters of the three subcooled water inlets are 3.74 mm; the inlet inner diameter of the soda mixture is 966 mm, while the inlet outer diameter is 1266 mm; finally, the width of the downcomer is 23 mm.

**Figure 2.** Schematic diagram of the physical model.

The inlet velocity of subcooled water was set as 11 m/s, while the flow mixed with saturated water and vapor was 0.12 m/s. The temperature of the mixture was the saturation temperature under the operating pressure, and the vapor rate of the mixture was 0.5%.

The commercial software ANSYS/ICEM [19] was used to draw three kinds of hexahedral structured grids with different quantities and the mesh quantities were 723,030, 1,584,664, and 2,441,390, respectively. The grids of the subcooled water inlet were encrypted, the number of boundary layer grids was 10, and the y+ values corresponding to three different numbers of grids in the numerical calculation process were 70-90, 15-20, and 10, respectively. Based on the above three different amounts of grids, the temperature at the center of the section at the top of the downcomer was monitored over time. As shown in Figure 3, when the mesh quantity was 1,584,664 and 2,141,390, the temperature changed with time with the same trend with most of the errors being within 1%. In order to save time and cost, the second grid number was adopted for further simulations.

**Figure 3.** Temperature curve of the top section at different times.

Where Vs is the inlet velocity of the vapor mixture; Vm is the inlet velocity of the subcooled water; Ts is temperature of the vapor mixture; Tw is temperature of the subcooled water; and P is the pressure at the coordinate point.

#### *2.2. Mathematical Model*

A VOF multiphase model tracked the location of the vapor-liquid interface in real time, while a LES model [18] was employed to capture the pulsation characteristics of turbulence.

#### (1) Control equation

The VOF multiphase flow model control equation includes a continuity equation, a momentum equation, and an energy equation, as follows.

The continuity equation is [18]:

$$\frac{1}{\rho} \left[ \frac{\partial}{\partial t} (\alpha\_{\text{q}} \rho\_{\text{q}}) + \nabla \cdot (\alpha\_{\text{q}} \rho\_{\text{q}} \mathbf{u}\_{\text{q}}) \right] = \mathbf{S}\_{\alpha\_{\text{q}}} + \sum\_{\mathbf{P}=1}^{n} (\dot{\mathbf{m}}\_{\text{P}\mathbf{q}} - \dot{\mathbf{m}}\_{\text{P}\mathbf{P}}) \tag{2}$$

$$\text{With } \sum\_{\mathbf{q}=1}^{n} \alpha\_{\mathbf{q}} = 1 \tag{3}$$

where uq is the velocity of phase q, m/s; αq is the volume fraction of phase q; . mpq is the mass transfer from phase p to phase q, kg/s; . mqp is the mass transfer from phase q to phase p, kg/s; and Sαq is the quality source term.

The momentum equation is [18]:

$$
\left[\frac{\partial}{\partial \mathbf{t}} (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \mathbf{u})\right] = -\nabla \mathbf{P} + \nabla \cdot \left[\mu \left(\nabla \mathbf{u} + \nabla \mathbf{u}^{\mathsf{T}}\right)\right] + \rho \mathbf{g} + \mathbf{F} \tag{4}
$$

where ρ is the density, kg/m3; μ is the dynamic viscosity, kg/(m·s); u is the velocity, m/s; F is the volume force. Only the surface tension of two phases was considered in this numerical simulation.

The continuous surface force (CSF) model was adopted to determine the tension of the liquid phase:

$$\mathbf{F} = \mathbf{F}\_{\text{CSF}} = \sigma \kappa \nabla \mathbf{a} \tag{5}$$

where κ is the curvature and σ is the surface tension coefficient.

The density and dynamic viscosity [18] used in Equation (4) are defined in Equations (6) and (7):

$$
\mathfrak{p} = \mathfrak{a}\_2 \mathfrak{p}\_2 + (1 - \mathfrak{a}\_2)\mathfrak{p}\_1 \tag{6}
$$

$$
\mu = \alpha\_2 \mu\_2 + (1 - \alpha\_2)\mu\_1 \tag{7}
$$

The energy equation can be written as [18]:

$$\left[\frac{\partial}{\partial t}(\rho \mathbf{E}) + \nabla \cdot (\mathbf{u}(\rho \mathbf{E} + \mathbf{P}))\right] = \nabla \cdot [(\mathbf{k\_{eff}} \nabla \mathbf{T})] + \mathbf{S}\_{\mathbf{h}} \tag{8}$$

where E is the energy, p is the pressure, Pa; T is the temperature, K; keff is the effective thermal conductivity of the fluid, w/(m·K); and Sh is the energy source term.

(2) LES turbulence model

In the LES model, the n-s equation of incompressible fluid is filtered and its corresponding equation is obtained as follows [18].

$$\frac{\partial \rho}{\partial t} + \frac{\partial \rho \overline{\mathbf{u}}\_{\mathrm{i}}}{\partial \mathbf{x}\_{\mathrm{i}}} = 0 \tag{9}$$

$$\frac{\partial \rho \overline{\mathbf{u}\_{\mathrm{i}}}}{\partial \mathbf{t}} + \frac{\partial \rho \overline{\mathbf{u}\_{\mathrm{i}}} \overline{\mathbf{u}\_{\mathrm{j}}}}{\partial \mathbf{x}\_{\mathrm{j}}} = -\frac{\partial \overline{\mathbf{P}}}{\partial \mathbf{x}\_{\mathrm{i}}} - \rho\_{0} \mathbb{B} \left(\mathbf{T} - \mathbf{T}\_{0}\right) \mathbf{g} + \frac{\partial}{\partial \mathbf{x}\_{\mathrm{j}}} \left[ (\mu + \mu\_{\mathrm{t}}) \left(\frac{\partial \overline{\mathbf{u}\_{\mathrm{i}}}}{\partial \mathbf{x}\_{\mathrm{j}}} + \frac{\partial \overline{\mathbf{u}\_{\mathrm{i}}}}{\partial \mathbf{i}}\right) \right] \tag{10}$$

where μ is the dynamical viscosity, kg/(m·s); μt is the turbulent viscosity; and β is the coefficient of thermal expansion.

The Smagorinsky–Lilly [18] model was used for small scale vortices:

$$
\mu\_t = \rho \mathcal{L}\_s^2 \left| \overline{\mathbf{S}} \right|^2 \tag{11}
$$

where & & &S & & & = , 2SijSji. Ls is the sublattice-scale mixing length described as [18]:

$$\mathcal{L}\_{\mathsf{s}} = \min\{\mathsf{xd}, \mathsf{C}\_{\mathsf{s}} \mathsf{V}^{1/3}\} \tag{12}$$

where κ is the Von Karman constant, set as 0.42; d is the distance between the cell and the nearest wall; V is the volume of the cell; and Cs is the Smagorinsky constant, set as 0.1.

#### (3) Condensation model

The condensation process in the downcomer was determined by using the condensation double-resistance model, which is based on the heat balance principle and can be used to calculate the condensation phenomenon caused by heat transfer between the two phases of pure material flow. The condensation double-resistance model assumes that one phase (vapor phase) is dispersed as a small spherical bubble with a variable diameter at the intersection interface. The heat and mass transfer occur on the bubble surface, while thermal resistance exists simultaneously at the intersection interface. The energy transfer of steam direct contact condensation depends on the phase-change interfacial area, the interfacial heat transfer coefficient, and the interfacial mass transfer rate.

In the liquid phase, the vapor-liquid interface area is [18]:

$$\mathbf{A}\_{\rm ifg} = \frac{6\alpha\_{\rm ig}}{\mathbf{d}\_{\rm ig}} \tag{13}$$

where αig is the steam volume fraction and dig is the average bubble diameter, which is given by Anglart and Nylundl [20]:

$$\mathbf{d}\_{\rm ig} = \begin{cases} 1.5 \times 10^{-4} & \theta\_{\rm i} > 13.5\\ \frac{\Phi\_{\rm i}(\theta - \theta\_{\rm 0}) + \Phi\_{\rm 0}(\theta\_{\rm 1} - \theta)}{\Phi\_{\rm i} - \Phi\_{\rm 0}} & 0 < \theta\_{\rm i} < 13.5\\ 1.5 \times 10^{-3} & \theta\_{\rm i} \le 0 \end{cases} \tag{14}$$

The mean bubble diameter is used to calculate the liquid Reynolds number and Nusselt number in the condensing double-resistance model.

The heat transfer coefficient of the liquid phase can be determined as [18]:

$$\mathbf{h}\_{\rm if} = \frac{\mathbf{k}\_{\rm if} \mathbf{N} \mathbf{u}\_{\rm fi}}{\mathbf{d}\_{\rm ig}} \tag{15}$$

Nusselt number of the liquid phase is determined according to the correlation given by Hughmark [21]:

$$\mathrm{Nu\_{\bar{\mathrm{fl}}}} = \begin{cases} 2.0 + 0.6 \mathrm{Re}^{0.5} \mathrm{Pr}^{0.33} & 0 \le \mathrm{Re} < 766.06 \\\\ 2.0 + 0.27 \mathrm{Re}^{0.62} \mathrm{Pr}^{0.33} & \mathrm{Re} \ge 766.06 \end{cases} \tag{16}$$

The heat transfer quantity on the liquid phase side of per unit volume [15] can be calculated as:

$$\mathbf{H}\_{\rm if} = \mathbf{h}\_{\rm if} \mathbf{A}\_{\rm ifg} \tag{17}$$

The changes in fluid properties can cause an abrupt increase in the condensation rate. Koncar and Mavlo [22] proposed the so-called umbrella restriction, through the constraint of which the sudden change in the condensation rate can be limited in some special cases, approaching that of actual situations:

$$\mathrm{pH\_{if}} = \min\left[\mathrm{H\_{if}}, 1753, 9 \max\left[4.724, 472.4 \alpha\_{\mathrm{ig}} \left(1 - \alpha\_{\mathrm{ig}}\right)\right]\right] \times \max\left[0, \min\left(1, \frac{\alpha\_{\mathrm{ig}} - 1.0 \times 10^{-10}}{0.1 - 1.0 \times 10^{-10}}\right)\right] \tag{18}$$

*Energies* **2019**, *12*, 3650

The zero-resistance [18,20] model is applied from the interface at the vapor phase. The temperature on the vapor phase side of the interface approximates the saturation temperature at the operating pressure.

The heat flux of convective heat transfer from the phase-change interface to the liquid phase side is [18]:

$$\mathbf{q}\_{\rm if} = \mathbf{h}\_{\rm if}(\mathbf{t}\_{\rm is} - \mathbf{t}\_{\rm if}) \tag{19}$$

The total heat flux from the phase change interface to the liquid-phase side is [18]:

$$\mathbf{Q}\_{\rm il} = \mathbf{q}\_{\rm if} + \mathbf{m}\_{\rm isl} \mathbf{H}\_{\rm il} \tag{20}$$

where misl is the condensation rate, Hil is the saturated water enthalpy, whereas the total heat flux from the phase-change interface to the vapor phase side corresponds to [18]:

$$\mathbf{Q\_{is}} = \mathbf{m\_{is}} \mathbf{H\_{is}} \tag{21}$$

His is the saturated steam enthalpy. Therefore, based on the total heat flux equilibrium at the phase change interface, it can be assumed that:

$$\mathbf{Q}\_{\rm li} = \mathbf{Q}\_{\rm is} \tag{22}$$

The steam condensation rate can be obtained from the above equations:

$$\mathbf{m}\_{\rm isl} = \frac{\mathbf{q}\_{\rm if}}{\mathbf{H}\_{\rm is} - \mathbf{H}\_{\rm il}} \tag{23}$$

#### *2.3. Boundary Conditions*

In the calculation of the boundary conditions, the void fraction of the inlet distributed evenly was set as 0.5%, while the inlet velocity was determined according to the flow and flow area. The effect of the wall thickness was neglected in the numerical model, and the adiabatic boundary condition of the solid wall was adopted [23].

The double-precision implicit unsteady solver, which was developed based on pressure, was selected for the numerical calculations. The SIMPLE [13] algorithm was adopted for the coupling of pressure and velocity, while the PRESTO [13] discrete format was adopted for the pressure term. The momentum equation, energy equation, and volume fraction were discretized by a bounded central difference, second-order upwind, and geometric reconstruction, respectively. The VOF two-phase flow model used an explicit discrete format. Since surface tension was involved in the calculation, the implicit body force had been chosen in the body force formulation to improve the convergence of solutions. The Smagorinsky-Lilly [18] subgrid model was selected for the LES turbulence model. The residual convergence standard of the energy equation was set as 1 × 10−6, while the residual convergence standard of other variables was set as 1 × 10−3.