Heat Transfer Analysis Using Thermofluid Network Models for Industrial Biomass and Utility Scale Coal-Fired Boilers

Abstract

Integrated whole-boiler process models are useful in the design of biomass and coal-fired boilers, and they can also be used to analyse different scenarios such as low load operation and alternate fuel firing. Whereas CFD models are typically applied to analyse the detail heat transfer phenomena in furnaces, analysis of the integrated whole-boiler performance requires one-dimensional thermofluid network models. These incorporate zero-dimensional furnace models combined with the solution of the fundamental mass, energy, and momentum balance equations for the different heat exchangers and fluid streams. This approach is not new, and there is a large amount of information available in textbooks and technical papers. However, the information is fragmented and incomplete and therefore difficult to follow and apply. The aim of this review paper is therefore to: (i) provide a review of recent literature to show how the different approaches to boiler modelling have been applied; (ii) to provide a review and clear description of the thermofluid network modelling methodology, including the simplifying assumptions and its implications; and (iii) to demonstrate the methodology by applying it to two case study boilers with different geometries, firing systems and fuels at various loads, and comparing the results to site measurements, which highlight important aspects of the methodology. The model results compare well with values obtained from site measurements and detail CFD models for full load and part load operation. The results show the importance of utilising the high particle load model for the effective emissivity and absorptivity of the flue gas and particle suspension rather than the standard model, especially in the case of a high ash fuel. It also shows that the projected method provides better results than the direct method for the furnace water wall heat transfer.

1. Introduction

Integrated whole-boiler process models are useful in the analysis of biomass and coal-fired boilers. It can support the design process, serve as surrogate models for on-line monitoring and control, and allow investigation of different operational scenarios, such as the impact of fuel type and quality, boiler cleanliness, or off-design performance such as low load operation. Although computational fluid dynamics (CFD) can resolve the complex heat transfer and fluid flow characteristics in detail, it is not feasible for whole-boiler performance analysis. This requires one-dimensional thermofluid network models. These incorporate zero-dimensional furnace models combined with the solution of the fundamental mass, energy, and momentum balance equations for the different heat exchangers and fluid streams such as the combustion air, flue gas and steam.

This approach is not new, and there is a large amount of information already available in textbooks and technical papers. However, it is often fragmented and incomplete and therefore difficult to follow and apply. It also contains errors that create uncertainty unless there is a good understanding of the underlying fundamentals. The aim of this paper is therefore to: (i) provide a review of recent literature to show how the different approaches to boiler modelling have been applied; (ii) to provide a review and clear description of the thermofluid network modelling methodology, including the simplifying assumptions and its implications; and (iii) to demonstrate the methodology by applying it to two case study boilers with different geometries, firing systems and fuels at various loads and comparing the results to site measurements, which highlight important aspects of the methodology. There are three major modelling approaches used in boiler analysis, namely CFD modelling, one dimensional thermofluid network modelling, and coupled simulations where the detailed three-dimensional CFD models are coupled to one-dimensional network models and solved simultaneously. The rest of Section 1 below provides a review of recent literature to show how these different approaches have typically been applied, of which some papers are by the current authors. Following this, Section 2 provides a detailed description of the network modelling methodology that highlights the underlying assumptions. Section 3 reports new results of two case studies by the current authors where the network methodology has been applied and the results compared to site measurements, which highlight the impact of the various assumptions. Finally, Section 4 provides a summary and conclusions.

1.1. CFD Modelling

CFD has been used extensively for the analysis of solid fuel combustion systems. It can adequately resolve the complex fluid flow and heat transfer phenomena, as well as the combustion processes [1]. This has allowed various researchers to model full-scale boilers operating under steady-state conditions, provided that sufficient boundary conditions and computational resources are available. Laubscher and Rousseau [2] utilised a CFD modelling methodology, developed in ANSYS Fluent^® 2019, to evaluate the thermal performance of the evaporator and various superheaters of a 620 MW_e subcritical boiler at full and partial boiler loads. (Some of the results of this study are used in Section 3.1 to validate the application of the thermofluid network modelling methodology). Laubscher and van der Merwe [3] successfully used a CFD model to resolve the combustion phenomena and heat transfer to the various heat exchanger surfaces of a semi-suspension fired bagasse boiler. The model was also developed in ANSYS Fluent^® 2019 and included custom sub-models developed to capture the fuel particle grate interactions and non-spherical fuel particle drag effects in the boiler furnace. (Some of the results of this study are also used in Section 3.2 to demonstrate the applicability of the thermofluid network modelling methodology).

Constenla et al. [4] used the CFD modelling software ANSYS Fluent^® to resolve the flow characteristics and real operating conditions of a 350 MW_e tangentially fired pulverised coal furnace, with the results providing valuable insight into to combustion stability. In addition, CFD provides insights into the fireside operational aspects, such as pollution control and its mitigation. Modlinski et al. [5] utilised a CFD modelling approach to validate its capabilities in predicting the particle and gas distribution temperatures in the combustion chamber of a large utility scale boiler, since knowledge of the temperature field and the elimination of high-temperature zones during operation is essential in reducing a power plants NOx emissions. Rago et al. [6] investigated the effects of using swirl stabilised burners on NOx formation distributions. This was achieved with CFD modelling that incorporated the detailed chemical and kinetic characteristics found in the combustion processes.

CFD plays an important role in investigating the retrofitting of boilers using new technological advances in the industry. Li et al. [7] made use of CFD to investigate the combustion characteristics and NOx formation of a burner retrofitted in a 330 MW_e low volatile coal fired boiler. The model was validated using site data, with the simulated results showing a 15% reduction in the NOx formation for multiple boiler loads when using the retrofitted burners. Using an open source CFD code, Madjeski [8] modelled the pulverised coal combustion in a 225 MW_e power plant to illustrate the optimal positioning of the over firing air nozzles. The CFD model was validated against site data under normal operating conditions and the simulated results confirmed the correct position of the over-fired air nozzles in the boiler furnace. Li et al. [9] presented a novel burner arrangement that was proposed to improve the performance of a 660 MW_e tangentially fired boiler. A CFD model was used to investigate the effects and characteristics of the burner arrangement under various operating conditions, with the results indicating the optimal configurations. Furthermore, He et al. [10] utilised the CFD package, ANSYS Fluent^®, to analyse the overheating issues experienced in the reheater and superheater pendants of a pulverised coal power plant. Various retrofitting cases were performed, with the best case being observed to modify the flue gas flow field entering the reheater and superheater zones, thus minimising the surface overheating effects.

The Eulerian–Lagrangian approach is acknowledged as the conventional CFD approach for modelling solid fuel combustion systems [11]. The approach uses a Eulerian description of the gaseous phase, where the finite volume approach is used to derive the governing equations in a stationary frame of reference. The solid fuel particles are modelled using a Lagrangian frame of reference whereby tracking of the particle trajectories is conducted for each particle injection. The particle source terms are calculated by averaging over a multitude of particle trajectories [12], typically leading to variability of the continuous phase source terms. This requires significant under-relaxation of the solid-phase source terms in the continuous phase equations, resulting in extended computation times [11].

Many researchers have used CFD as a tool to provide insight into the fire side and particulate effects of pulverised fuel combustion systems. Utilising the conventional CFD approach, Laubscher and Rousseau [13] conducted a study to investigate the effects of secondary air swirl direction on the heat absorption of the furnace, platen and high-temperature superheaters in a 620 MW_e utility scale boiler. Six swirl arrangements were investigated. It was found that swirling the burners in the same direction increases the furnace heat absorption compared to the current operational arrangement. Modlinski [14] utilised a CFD model to conduct a numerical investigation of a tangentially fired pulverised boiler to investigate the retrofitting of the boiler with rapid ignition jet swirl burners. The study focused on the effects of the new firing systems on the combustion performance, flow, and heat transfer in the furnace. The results highlighted the effect of the new burners on reducing the likelihood of corrosion and slagging occurring in the furnace. In addition, a stable flame operation was achieved even without the use secondary air swirl.

Liu et al. [15] investigated the co-firing of low-volatile and high-volatile coals for a 330 MW_e tangentially fired pulverised coal boiler. The CFD model was validated against experimental results. The primary objective was to simulate multiple cases that consider the co-firing ratio, particle sizes, injection position and the unburnt carbon content. The findings show that a substantial increase in the unburnt content and NOx emission was observed when the co-firing ratio increases from 0 to 100%. It also showed that unburnt carbon content could be reduced by 83% with only marginal costs involved by addressing the particle size associated with the grinding of the low-volatile coal in the mill.

Low-load operation studies have recently been performed using CFD as the primary simulation tool. The studies have focused on ensuring combustion stability, minimising harmful emissions, and investigating the gas and solid flow interactions at low loads. Belosevic et al. [16] conducted a numerical study into the combustion behaviour of pulverised coal in a full scale 350 MW_e tangentially fired utility boiler. Experiments highlighted the complex fluid and interphase exchange of the dispersed particles. Findings show that the particle combustion, depending on particle size, greatly influenced the vertical position of the flame ball. In addition, the coal quality affects areas such the furnace exit gas temperature (FEGT), products of combustion, and thermal loading on the evaporator walls. Similarly, Hernik et al. [17] investigated the effects of using different mill system configurations at a minimum boiler load of 40%. The most favourable mill system configuration was selected based on the case that exhibited suitable combustion stability, which was defined for the study as the effective mixing of the fuel and oxidiser as well as maintaining sufficient temperature distributions.

CFD modelling approaches have been used successfully by researchers to provide engineering insight into boiler operational conditions. Using CFD does not interfere with the regular operation of the plant, and it can resolve the complex heat transfer and fluid flow characteristics with sufficient accuracy and detail for the fire side components. However, the initial design and analysis of boiler systems as well as addressing what-if questions in an operational environment require fast simulation tools that can resolve both the fire side and working fluid components. Using CFD to resolve the integrated whole-boiler performance is not feasible at present due to the long computational times involved. Therefore, integrated one-dimensional thermofluid network models that can resolve the fundamental mass, energy and momentum balance equations together with the component characteristics for the various heat exchangers are a valuable alternative.

1.2. One-Dimensional Process Modelling

One-dimensional process modelling techniques have been developed by various researchers to capture the steady-state and dynamic response of coal fired power plants, since the modelling approach provides efficient use of computational resources [18]. Kuronen et al. [19] utilised the one-dimensional process modelling software APROS to model a 750 MW_e once-through ultra-supercritical boiler. The results were validated for both steady state and transient applications using the measured site data for a load change from 87% to 100% maximum continuous rating. The model was empirically tuned using a lumped variable approach. The primary focus of the study was to develop a model for testing and design purposes for future flexible operation.

The furnaces of coal fired power plants boilers are comprised of a multitude of complex and interacting phenomena, such as the combustion dynamics, the gas-solid interactions, fluid dynamics, and radiation heat transfer. As a result, the heat transfer in the furnaces incorporate many non-uniformities, thus, many researchers utilise a lumped parameter analysis when dealing with the furnace section, such as the Gurvich/Blokh method. This method was first proposed by Blokh in 1984 [20] and has subsequently been used to create boiler design manuals such as those provided in references [21,22]. A one-dimensional process model of the start-up system of an existing 600 MW_e supercritical once-through boiler was developed by Deng et al. [23]. The model incorporates the two-phase homogeneous flow model to resolve the water/steam networks. The simulation software APROS was used to build the entire model, including the flue gas and fire side interactions. The model was able to provide sufficient resolution of the steam characteristics during start-up with a 2.29% relative error between the design and simulated steam exit temperatures and flow rate.

Zima et al. [24] analysed a subcritical boiler using a process modelling approach, with the furnace heat transfer being captured using an energy and mass balance combined with the empirical relations described in the Gurvich/Blokh method. Downstream components were modelled in a similar manner using mass and energy control volumes. The model heat fluxes and exit flue gas temperature were verified using a CFD model for a steady state boiler loads of 40% and 100%. The study aided in development of an in-house program to be used in support of the boiler design process.

Starkloff et al. [25] utilised a one-dimensional process modelling approach to capture the full water and flue gas interactions of a coal fired power plant using the APROS software package. The objective of the study was to develop a fast model to study the flexible operation of large coal fired power plants. Minimal boundary conditions were utilised by the model, which included the inlet temperature and mass flow rate to the condensers, the ambient air conditions, and the fuel composition. The model was used to investigate dynamic operations, namely a turn down of 100% to 27.5% load, and a corresponding ramp up scenario from 27.5% to 100% load, with results showing that the model can simulate the dynamic behaviour of the plant with high accuracy.

Similarly, Oko and Wang [26] modelled a 500 MW_e sub-critical boiler using a one-dimensional approach that was validated for 70%, 80%, 94,4% and 100% steady-state operating conditions. The furnace used a zero-dimensional model, with only the radiation heat transfer component being considered. The full boiler including all the mechanical components and control logic was incorporated in the model. The software package gPROMS^® was utilised to capture the whole boiler model. The model was validated for the 100% steady state load case using measured plant data. Further dynamic simulations were conducted comparing a step change and ramp change control strategy for implementing load changes, with the latter being best suited for the control strategy since minimal fluctuations were observed.

Rousseau and Laubscher [27] developed a novel thermofluid model for the membrane walls of a pulverised fuel boiler using an equivalent one-dimensional thermal resistance network to resolve the unknown temperatures through the membrane wall. Detailed CFD results were used to validate the model with the overall UA value being within 1%. A further parametric study was conducted to investigate different geometries and operating conditions that cover the full range of operations for both subcritical and supercritical boilers. Trojan [28] developed a one-dimensional process model of a furnace combustion chamber assuming that the flue gas temperature in the combustion chamber was uniform and the specific heat of the flue gas was a function of temperature. The combustion chamber approach proposed in the study is analogous to the Gurvich/Blokh method. The downstream heat exchanger geometry was considered during the modelling setup and the flue gas temperatures at the exit of these heat exchangers were determined using the effectiveness NTU method. The model was verified against measured data with very good agreement between the calculated and measured values at certain locations in the boiler.

Due to the computational efficiency of one-dimensional process models the investigation of various control strategies and boiler retrofitting can be investigated in an efficient manner. Chandrasekharan et al. [29] utilised a one-dimensional mathematical model to optimise the control design and operational strategy of a 250 MW_e coal fired power plant. The economiser and steam drum of a subcritical boiler were modelled using mass and energy balances derived from first principles. The focus was on the open loop control of the steam drum, to retrofit the boiler with bypass valves after the primary and secondary superheater to increase the steam drum temperatures. A retrofitting study conducted by Chantasiriwan [30] focused on determining the optimum position of an air heater in a conventional utility scale boiler. The Gurvich/Blokh approach was used to resolve the furnace exit gas temperature and the heat load to the furnace, with the downstream heat exchangers being modelled using conventional mass and energy balances. Three cases were simulated using the boiler process model with additional heat exchanger models used to resolve the air heater and fuel dryer. A cost-benefit analysis was conducted, but the process model was not explicitly validated.

Gradziel [31] developed a one-dimensional model to investigate the evaporator of a natural circulating boiler including the phase transition phenomena. An in-house program code was developed for the study. The model incorporated the thermal response and flow phenomena with only the evaporator walls being considered. The furnace radiative heat transfer distribution was determined using the zone method, with the distinct nose shaped heat flux profile being adopted to represent the complex thermal loading arising from the combustion chamber. Experimental results were obtained to verify the evaporator model for steady state operation and during start-up and shutdown, with satisfactory results being obtained. Similarly, Wang et al. [32] considered a unique annular furnace structure for a 600 MW_e circulating fluidised bed boiler. The water wall circulation was analysed using the appropriate process modelling techniques that capture the fundamental mass, energy, and momentum balance equations. The model was used to investigate the flow distribution, the thermal-hydraulic response, pressure drop and surface temperatures for the various boiler loads, namely 100%, 75% and 30% of maximum continuous rating. The results show that a vertical water wall design with smooth tubes offers the optimal performance across the simulated load ranges.

As with CFD simulations, researchers have used one-dimensional process models to evaluate various types of boiler configurations. One-dimensional process models have been used extensively in the modelling of circulating fluidised beds, due to the efficient solution times and the use of minimal computational resources. Guoli et al. [33] investigated the thermal-hydraulic characteristics of an ultra-supercritical steam cycle using circulating fluidised bed combustion technology to burn low-grade fuels. The mathematical model of the evaporator system was established based on energy and mass conservation equations and incorporated an effective pressure balance for the parallel connecting tubes of the furnace walls. In addition, a semi-empirical bed-to-wall heat flux model and a set of empirical correlations available for heat transfer and hydraulic resistance calculation were embedded in the thermal-hydraulic model. The model was run for operating conditions of 100%, 75% and 25% of nominal load. The results showed that the conceptual design of the evaporative system is acceptable for a 660 MW_e ultra-supercritical circulating fluidised bed boiler. Similarly, Alobaid et al. [34] developed a sophisticated dynamic process model of a 1 MW_th circulating fluidised bed test facility. The model incorporates the details of the air supply, the fluidised bed, the flue gas path, the water-cooling system, and the implemented control structures. The model was validated and tuned using experimental data. Further load change studies were implemented, and the validated model allowed the researchers to numerically obtain parameters that cannot be fully measured.

Alobaid et al. [18] conducted a comprehensive review in 2016 of dynamic simulation models, its development and application to various types of thermal power plants, including coal-fired, nuclear and combined cycles, as well as municipal waste incineration. The fundamental dynamic balance equations are presented together with its application to various components, from basic pipe models and complex models including two-phase flow, to the essential electrical and automation components. It presents an overview of work completed in dynamic analysis of coal fired power plants in relation to disturbances, start up, flexibility and the oxy fuel concept. It showcases the use of one-dimensional process models for dynamic modelling applications, on-line calculations, whole boiler operation simulations, and the efficient resolution of the water and steam side components. However, it does not provide all the details necessary to enable the development of integrated whole-boiler process models for biomass and coal-fired boilers.

1.3. Coupled Modelling Approaches

Coupled modelling approaches utilise the advantages of CFD and one-dimensional process models to resolve the interactions between the fireside and waterside components for various boiler types [35]. CFD allows for the detailed resolution of the combustion and radiation heat transfer processes in the furnaces, while one-dimensional process models can adequately resolve the complex water/steam network in an efficient manner. Park et al. [36] developed a coupled CFD and one-dimensional water-steam process model for a tangentially fired boiler. ANSYS CFX^® was used to model the fireside effects of the furnace and the main superheaters and PROATES was used to capture the water/steam side components. The coupling was achieved via a direct method whereby the walls were used as a data-transfer interface between the models. The model was validated against a 100% steady-state load case using measured plant data. The coupled model was further used to provide operational guidance and allowed the plant engineers to optimise burner settings and the heat exchanger performance.

Schuhbauer et al. [37] show the full development of a detailed tower type boiler simulation using a coupled CFD and one-dimensional process model. APROS^® is used to model the waterside, while the combustion and three-dimensional effects are captured in ANSYS Fluent^®. The coupling interface is the heat exchanger walls where the CFD model sends the heat flux as a boundary condition to APROS^®, which would update the outer tube temperatures. These temperatures would then be transferred back to the CFD model to update the boundary conditions. MATLAB^® was used as an intermediate data transfer link, since APROS^® does not allow for direct coupling. The convective components in the CFD model were modelled using the porous media approach, since the detailed modelling of the tube bundles was deemed computationally expensive. The issue with porous media model in modelling a tower-type boiler is that no visible surface is given, which resulted in very low radiation exchange between the combustion chamber and the tube bundles. A simple model was subsequently developed to account for the missing radiation flow with reasonable results being obtained in accordance with the design data. Validation was performed for the full load case, and from there, part load cases were simulated. The simulation results are within 6% of the design data.

The coupled simulation technique has allowed researchers to commence detailed studies for the various convective and radiative heat exchangers typically found in utility-scale boilers, since the CFD component can provide detailed distributions of the radiative and convective heat fluxes distributed on the tube surfaces. This, combined with the one-dimensional process models of the waterside, allows researchers to investigate a range of impacts. Laubscher and Rousseau [38] developed a coupled model utilising a one-dimensional discretised two-phase model of the water/steam side of the evaporator and radiative superheaters and a CFD model capturing the furnace combustion and heat transfer. The model was used to investigate the impact of particle radiative effects. The study highlighted the importance of using conversion-dependent particle radiative properties for high-ash content fuels, with the conversion-dependent properties providing better results corresponding to the experimental data.

Chen et al. [39] compared a coupled and decoupled approach of the supercritical steam in the pendant superheater of a 1000 MW_e tangentially fired boiler. The coupled approach utilises a CFD model and a one-dimensional hydrodynamic model which accepts the temperature distribution and energy fluxes to the pendant walls, while the decoupled model primarily uses a CFD with a uniform temperature profile applied to the pendant walls. The coal combustion and flue gas flow distribution were modelled using CFD, and only the pendant superheater was modelled in a one-dimensional process model. It was found that the coupled method offers the best approach to capture the energy flux and wall temperature distributions. However, the model lacks validation for the given operating conditions.

Yang et al. [40] modelled a 20 MW_th supercritical CO₂ boiler by combining a one-dimensional analytical model for the supercritical CO₂ side and a CFD model for the combustion and heat transfer processes. It employed ANSYS Fluent^® user-defined functions to couple the CFD and one-dimensional model. The 1D model was developed using the governing equations and the appropriate Nusselt number correlations. The model was able to capture the non-uniform nature of the heat fluxes in the superheater and reheater tube walls, providing a better understanding of the heat transfer processes and safe operation of the 20 MW_th supercritical CO₂ boiler.

Fei et al. [41] utilised a coupled CFD and a one-dimensional process model in a boiler retrofitting study. CFD was used to capture the coal combustion and heat transfer to the waterwalls and heat exchangers. A reduced order model (ROM) was developed for the integration of the CFD model into the full plant process model. CFD was used to generate and develop a ROM using kriging interpolation. This ROM was integrated with the full plant process model to offer faster simulation times in the retrofitting study. A range of air-coal and oxy-coal conditions with different power loads (350–500 MW_e) were simulated. The results show that it is possible to retrofit the air-coal firing to oxy-coal firing while maintaining the original design performance of air-coal firing.

Rousseau and Laubscher [42] investigated the impact of using very high ash coal versus the design coal on the heat transfer distribution in a utility-scale boiler. A coupled simulation was employed, using CFD to model the fireside and a one-dimensional thermofluid network model capturing the water/steam side. The results show that the current use of low-quality coal impacts the furnace heat uptake and steam generation rate when compared to the design fuel case. The main reason is the increased attenuation induced by the increased particle load to the domain.

Coupled simulation techniques provide a valuable investigative tool for researchers and engineers. However, the CFD modelled component limits the application of the models in the analysis of transient events and the development of online monitoring systems [43] due to the computational burden and long simulation times. The use of one-dimensional process models of the whole boiler system can counteract these shortfalls provided that the correct correlations and modelling approaches are incorporated.

2. Modelling Methodology

This section outlines the integrated boiler process modelling methodology that is built on a one-dimensional thermofluid network approach. Since we are interested in analysing the integrated performance of the system rather than resolving the details associated with each component, we approximate the various components as one-dimensional elements having one inlet and one outlet, which are connected at nodes to form a thermofluid network. The analysis of each component consists of the following generic elements:

• The simultaneous solution of the fundamental balance equations of mass, energy and momentum written between the inlet and the outlet.
• The performance characteristics of the specific component, which may include the rate of heat transfer to or from the fluid, the rate of work performed on or by the fluid, and the pressure drop or pressure rise in the fluid as it moves through the component. We will call the set of equations needed to describe this performance the ‘component characteristic’ equations.
• The specific fluid property relationships or state equations.
• The boundary values that may be specified at the inlet and/or outlet.

2.1. Balance Equations

The generic form of the balance equations of mass, energy and momentum for steady-state incompressible flow through any thermofluid component may be written between the inlets ($i$) and outlets ($e$) as:

$$\begin{array}{cc}\hfill \sum {\dot{m}}_i=& \sum {\dot{m}}_e\hfill \\ \hfill \sum {\dot{m}}_i\left(h_i+\frac{1}{2}{v}_i^2+g{z}_i\right)+\dot{Q}=& \sum {\dot{m}}_e\left(h_e+\frac{1}{2}{v}_e^2+g{z}_e\right)+\dot{W}\hfill \\ \hfill p_i+\frac{1}{2}\rho v_i^2+\rho g{z}_i+\Delta p_M=& p_e+\frac{1}{2}\rho v_e^2+\rho g{z}_e+\Delta p_L\hfill \end{array}$$

(1)

In Equation (1), $\dot{m}$ (kg/s) is the mass flow rate, $h$ (J/kgK) is the static enthalpy, $p$ (Pa) is the static pressure, $\rho$ (kg/m³ ) is the fluid density, $v$ (m/s) is the velocity, $g$ (m/s² ) is the gravitational acceleration, $z$ (m) is the elevation, $\dot{Q}$ (W) is the rate of heat transfer to the fluid, $\dot{W}$ (W) is the rate of work performed by the fluid, $\Delta p_M$ (Pa) is the total pressure rise due to machine work, and $\Delta p_L$ (Pa) is the total pressure loss.

In the context of boiler heat transfer analysis, we usually exclude the detail analysis of pumps and fans and simply specify the required mass flow rates as boundary values. We may also neglect the kinetic energy terms and, if we are not interested in solving the natural convection flow rate and static pressure head, we may also neglect the potential energy terms. This then leads to the generic set of balance equations that may be applied between the outlet and inlet of each of the components in the boiler, namely:

$$\begin{array}{cc}\hfill \sum {\dot{m}}_e=& \sum {\dot{m}}_i\hfill \\ \hfill \sum {\dot{m}}_e{h}_e=& \sum {\dot{m}}_i{h}_i+\dot{Q}\hfill \\ \hfill p_e=& p_i-\Delta p_L\hfill \end{array}$$

(2)

Although there are an infinite number of possible boiler layouts and several different firing systems using many different fuels, it is possible to identify a collection of generic components that represent those found in almost all boilers. It is also possible to identify a simplified high-level layout that represents the connectivity between these different components within a typical boiler. Figure 1 shows a process flow diagram (PFD) that contains these generic components connected in such a simplified layout. We will make use of this PFD to facilitate the discussion of the modelling methodology.

There are three distinct networks, namely the combustion air flow network (the green lines with air nodes labelled from a,0 to a,2), the flue gas flow network (the red lines with gas nodes labelled from fg,0 to fg,7), and the steam (or water) flow network (the blue lines with steam nodes labelled from st,0 to st,5). In addition, we have the attemperator spray water flow (the blue arrow with node labelled “AT”), the fuel feed line (the black arrow with node labelled “F”), and the bottom ash (ba) extracted from the furnace. There could also be a second steam flow network, namely the reheat steam path, but that will not be included in the current discussion since the same principles can be applied there.

The combustion air enters at air node a,0 and is pre-heated in one or more heat exchangers in series, which in this illustration are represented by only two air heater components namely AH1 and AH2. It leaves AH2 at air node a,2 and flows towards the combustion zone within the furnace. The combustion zone is represented by the adiabatic flame temperature (AFT) component. The combustion air is mixed (and/or pre-mixed) with the fuel, which results in combustion and heat release, whereafter the hot products of combustion leave the combustion zone at gas node fg,0. At the same time, a fraction of the incombustibles in the fuel leaves as bottom ash while the rest, known as fly ash (fa), is carried along with the other products of combustion at gas node fg,0. From there, it passes through the furnace (FUR) component and exits at gas node fg,1. The products of combustion consist of a suspension made up of a hot gas mixture and ash particles, which is represented by the flue gas flow network.

Within the FUR component, heat is transferred from the flue gas suspension to the surrounding water walls, primarily via thermal radiation. The partially cooled suspension then enters a train of radiative and convective superheater and reheater heat exchangers. In this illustration, the heat exchanger train is represented by only two superheater components, namely SH1 and SH2, and these are separated by a cavity (CAV) component. Within these heat exchangers, the cavity heat is transferred from the flue gas suspension to the outside surfaces of the heat exchanger tubes and surrounding water walls via a combination of convection and radiation heat transfer. The flue gas leaves the superheater train at gas node fg,4.

From gas node fg,4, the flue gas enters an economiser (EC) heat exchanger. Here, it rejects more heat to the surfaces of the heat exchanger tubes, primarily via convection, since at this point, the flue gas temperatures are significantly lower than in the superheaters and reheaters. The mixture leaves at gas node fg,5 from where it enters the train of air heater heat exchangers. Here, it again rejects heat to the heat exchanger surfaces, predominantly via convection. The mixture leaves AH1 at gas node fg,7, from where it moves on to the flue gas cleaning systems and is ultimately released to the atmosphere via the stack.

The feed water enters at steam node st,0 and is first pre-heated in the EC heat exchanger, from where it leaves in a subcooled liquid state. In natural circulation boilers, the subcooled water is mixed with saturated liquid water taken from the bottom of the steam drum (SD). It then flows via a downcomer to the bottom of the boiler where it enters the membrane-type water walls that surround the furnace. From there, it flows upwards through the upper water walls that also surround the heat exchanger train and various cavities and then exits into the steam drum as a mixture of saturated liquid and vapour, where separation takes place. The saturated liquid leaves at the bottom of the SD and flows into the downcomer. The saturated vapour leaves at the top of the SD at steam node st,2 and flows towards the superheaters. The difference in the steam/water density between the downcomer and the water walls results in continuous re-circulation via natural convection.

From a process modelling perspective, the evaporator (EV) heat exchanger component consists of all the water walls and the SD combined, and it may also include the roof tubes at the top of the boiler if these are connected in series with the water walls. Therefore, pre-heated feed water enters the EV component as subcooled liquid at steam node st,1, and it leaves as saturated vapour at steam node st,2.

In once-through forced flow boilers, the water from the EC also flows to the bottom of the furnace, through the EV water walls to the superheaters. However, in once-through mode, the water only passes through the water walls once; then, it passes through separators and a collector vessel and on to the superheaters. Therefore, there is no re-circulation. At low loads (typically below 40% MCR), water from the collector vessel is also circulated through the EC and EV via a circulation pump. For this discussion, we will only consider the natural circulation case, but the modelling principles can equally be extended to once-through boilers.

After passing through SH1, the superheated steam leaves at steam node st,3 and is typically attemperated before it enters SH2, using high-pressure water from the feedwater pump. This controls the outlet steam temperature of SH2 at steam node st,5 to the required main steam set point value. The attemperation process is represented by the attemperator (AT) component that accounts for the mixing of the superheated steam and subcooled water, with the attemperated steam leaving at steam node st,4. The steam then absorbs more heat in SH2 and leaves at steam node st,5 at the desired main steam temperature, from where it flows through the turbine inlet valve to the high-pressure turbine.

The connectivity between the different components in the boiler is captured in the model via the set of mass, energy and momentum balance equations that can be written for the air, steam and flue gas flow networks, respectively, based on Equation (2).

2.1.1. Air Flow Network

Although the air flow network in our representative PFD is very simple, it is useful to describe the approach for setting up the network equations in some detail. The same approach can then be expanded to the more complex steam and flue gas flow networks.

The mass balance equations for the air flow network are given by

$$\begin{array}{l}{\dot{m}}_{a,0}={\dot{m}}_{a,ca}\\ {\dot{m}}_{a,1}-{\dot{m}}_{a,0}=0\\ {\dot{m}}_{a,2}-{\dot{m}}_{a,1}=0\end{array}$$

(3)

where ${\dot{m}}_{a,j}$ is the air mass flow rate at node $j$ and ${\dot{m}}_{a,ca}$ is the total combustion air flow rate provided by the primary and secondary air fans combined.

Equation (3) can be solved using a suitable direct or iterative method to yield the values of the air mass flow rate at each of the nodes. From Equation (3), it should be clear that the left-hand side effectively specifies the connectivity between the various components, i.e., flow into and out of each component, while the right-hand side contains the inlet boundary value.

The energy balance equations are given by

$$\begin{array}{l}{\dot{m}}_{a,0}{h}_{a,0}={\dot{m}}_{a,0}{h}_{a,ca}\\ {\dot{m}}_{a,1}{h}_{a,1}-{\dot{m}}_{a,0}{h}_{a,0}={\dot{Q}}_{a,AH1}\\ {\dot{m}}_{a,2}{h}_{a,2}-{\dot{m}}_{a,1}{h}_{a,1}={\dot{Q}}_{a,AH2}\end{array}$$

(4)

where $h_{a,j}$ is the static enthalpy of the air at node $j$, $h_{a,ca}$ is the static enthalpy of the incoming combustion air supplied by the primary and secondary air fans combined, and ${\dot{Q}}_{a,AH1}$ and ${\dot{Q}}_{a,AH2}$ are the rates of heat transfer to the air within AH1 and AH2, respectively.

Equation (4) can be solved to provide the values of static enthalpy of the air at each of the nodes. Note that again, the left-hand side specifies the connectivity and the right-hand side contains the inlet boundary value—and in this case also the rates of heat transfer (energy sources) to the air within each of the components.

The momentum balance equations are given by

$$\begin{array}{l}{p}_{a,0}=p_{a,ca}\\ p_{a,1}-p_{a,0}=-\Delta p_{a,AH1}\\ p_{a,2}-p_{a,1}=-\Delta p_{a,AH2}\end{array}$$

(5)

where $p_{a,j}$ is the static pressure of the air at node $j$, $p_{a,ca}$ is the static pressure of the combustion air supplied by the primary and secondary air fans combined, and $\Delta p_{a,AH1}$ and $\Delta p_{a,AH2}$ are the total pressure losses within AH1 and AH2, respectively.

The simultaneous solution of Equations (3)–(5) will yield the values for the mass flow rate, static enthalpy and static pressure at each of the nodes in the air flow network.

However, in order to solve the equations, we need to know the boundary values ${\dot{m}}_{a,ca}$, $h_{a,ca}$ and $p_{a,ca}$, as well as the source term values, namely the rate of heat transfer ${\dot{Q}}_a$ to the air within each of the components and the pressure drop $\Delta p_a$ within each of the components. Furthermore, an iterative approach will be needed where each set of equations are solved successively together with the source term values, which will have to be updated in between the iterations.

The values of the source terms must be determined using the component characteristic equations referred to earlier. To solve the component characteristic equations for each component, we will need to know the values of other secondary fluid properties besides enthalpy and pressure, such as temperature, density, and conductivity. For this, we need to solve the specific fluid property relationships that relate the value of each of these secondary properties to the enthalpy and pressure. The average values of enthalpy and pressure within each component are usually assumed to be the arithmetic mean of the known values at the inlet and outlet that are obtained by solving the balance equations.

We will leave the detailed discussion of the component characteristic equations and the fluid property relationships for later and will now go on to write the balance equations for the steam and flue gas flow networks. Since the focus of this paper is on heat transfer, we will only address the mass and energy balance equations while leaving the momentum balance equations for the reader to explore.

2.1.2. Steam Flow Network

The mass balance equations for the steam flow network are given by

$$\begin{array}{l}{\dot{m}}_{st,0}={\dot{m}}_{st,fw}\\ {\dot{m}}_{st,1}-{\dot{m}}_{st,0}=0\\ {\dot{m}}_{st,2}-{\dot{m}}_{st,1}=0\\ {\dot{m}}_{st,3}-{\dot{m}}_{st,2}=0\\ {\dot{m}}_{st,4}-{\dot{m}}_{st,3}={\dot{m}}_{AT}\\ {\dot{m}}_{st,5}-{\dot{m}}_{st,4}=0\end{array}$$

(6)

where ${\dot{m}}_{st,j}$ is the steam (or water) mass flow rate at node $j$, ${\dot{m}}_{st,fw}$ is the feed water flow rate provided by the feed water pumps and ${\dot{m}}_{AT}$ is the attemperator spray water flow rate. From the perspective of the mass balance equations, ${\dot{m}}_{st,fw}$ is a boundary value and ${\dot{m}}_{AT}$ is a source term. Note that we therefore have six equations and six unknowns, namely the mass flow rates at each of the steam nodes numbered 0 to 5.

The energy balance equations for the steam flow network are given by

$$\begin{array}{l}{\dot{m}}_{st,0}{h}_{st,0}={\dot{m}}_{st,0}{h}_{st,fw}\\ {\dot{m}}_{st,1}{h}_{st,1}-{\dot{m}}_{st,0}{h}_{st,0}={\dot{Q}}_{st,EC}\\ {\dot{m}}_{st,2}{h}_{st,2}-{\dot{m}}_{st,1}{h}_{st,1}={\dot{Q}}_{st,EV}\\ {\dot{m}}_{st,3}{h}_{st,3}-{\dot{m}}_{st,2}{h}_{st,2}={\dot{Q}}_{st,SH1}\\ {\dot{m}}_{st,4}{h}_{st,4}-{\dot{m}}_{st,3}{h}_{st,3}-{\dot{m}}_{AT}{h}_{AT}=0\\ {\dot{m}}_{st,5}{h}_{st,5}-{\dot{m}}_{st,4}{h}_{st,4}={\dot{Q}}_{st,SH2}\\ {\dot{m}}_{st,5}{h}_{st,5}={\dot{m}}_{st,5}{h}_{st,ms}\end{array}$$

(7)

where $h_{st,j}$ is the static enthalpy of the steam at node $j$, $h_{st,fw}$ is the enthalpy of the feed water from the feed water pump, ${\dot{Q}}_{st}$ is the rate of heat transfer to the water/steam within each of the components, $h_{AT}$ is the enthalpy of the attemperator spray water, and $h_{st,ms}$ is the enthalpy of the steam at the desired main steam temperature.

Some observations are necessary regarding Equation (7). Although the node numbering in the steam network is only from 0 to 5 (i.e., six nodes), there are actually seven nodes in the network, since it also includes the AT node. Therefore, there are seven equations and seven unknown variables. The seven unknown variables are the enthalpies at nodes 0 to 5 and the required attemperator spray water mass flow rate ${\dot{m}}_{AT}$ to ensure that we obtain the desired main steam temperature.

Note that this allows us to solve for the attemperator spray water flow rate needed to obtain the desired main steam enthalpy. This implies that we are building the spray water controls directly into the model via the solution of the balance equations. The value of ${\dot{m}}_{AT}$ obtained from the solution of the energy balance equations is then provided as a source term value to the mass balance equations Equation (6). The alternative approach of specifying ${\dot{m}}_{AT}$ and solving for $h_{st,ms}$ is also valid if the aim is to study the impact of different spray flow rates on the main steam temperature.

Another feature that requires consideration is the control of the water level in the steam drum. In practice, this is usually studied by varying the feed water flow rate via a control valve at the outlet of the feed water pump or by varying the rotational speed of the pump. From a process modelling perspective, this means that the feed water flow rate is adjusted to ensure that there will always be saturated vapour leaving the SD component. This control feature can also be built into the model as follows:

$${\dot{m}}_{st,fw}=\frac{{\dot{Q}}_{st,EC}+{\dot{Q}}_{st,EV}}{{h_g|}_{p_{SD}}-h_{st,fw}}$$

(8)

where ${h_g|}_{p_{SD}}$ is the saturated vapour enthalpy at the steam drum pressure $p_{SD}$. The value of ${\dot{m}}_{st,fw}$ obtained from the solution of Equation (8) is provided as a boundary value to the mass balance equations Equation (6).

2.1.3. Flue Gas Flow Network

There are two additional complications that must be considered for the flue gas flow network. The first is that it represents the flow of a suspension consisting of the flue gas mixture and the fly ash particles entrained within it. The second is that there is a possibility that air from the surroundings may leak into the gas ducts due to the negative gauge pressure maintained in the boiler. This is referred to as ingress air.

The fly ash particles are accounted for by specifying the fly ash mass flow rate ${\dot{m}}_{fa}$ leaving the combustion zone and entering the FUR at gas node 0. This value is then assumed to remain the same throughout all the components. The mass balance equations for the fly ash are therefore simply ${\dot{m}}_{fa,j}={\dot{m}}_{fa}$ for $j=0\dots 7$.

The ingress air is accounted for by allowing for a specified ingress air mass flow rate ${\dot{m}}_{ia}$ to be included at the inlet to each component, excluding the furnace. Therefore, the mass balance equations for the gas network are given by

$$\begin{array}{l}{\dot{m}}_{fg,0}={\dot{m}}_{fg,aft}\\ {\dot{m}}_{fg,1}-{\dot{m}}_{fg,0}=0\\ {\dot{m}}_{fg,2}-{\dot{m}}_{fg,1}={\dot{m}}_{ia,SH2}\\ {\dot{m}}_{fg,3}-{\dot{m}}_{fg,2}={\dot{m}}_{ia,CAV}\\ {\dot{m}}_{fg,4}-{\dot{m}}_{fg,3}={\dot{m}}_{ia,SH1}\\ {\dot{m}}_{fg,5}-{\dot{m}}_{fg,4}={\dot{m}}_{ia,EC}\\ {\dot{m}}_{fg,6}-{\dot{m}}_{fg,5}={\dot{m}}_{ia,AH1}\\ {\dot{m}}_{fg,7}-{\dot{m}}_{fg,6}={\dot{m}}_{ia,AH2}\end{array}$$

(9)

${\dot{m}}_{fg,j}$ is the mass flow rate of the gas mixture (excluding the fly ash) at each node $j$ and ${\dot{m}}_{fg,aft}$ is the mass flow rate of the gaseous combustion products leaving the combustion zone at the adiabatic flame temperature, which is a boundary value.

Note that ${\dot{m}}_{fg,aft}$ is a boundary value that will be calculated within the AFT component, while the ingress air flow rates at each component are known mass source terms.

The energy balance equations for the gas flow network include the effect of the ingress air as well as the energy associated with the fly ash within the mixture and is given by

$$\begin{array}{l}{\dot{m}}_{fg,0}{h}_{fg,0}={\dot{m}}_{fg,0}{h}_{fg,aft}\\ {\dot{m}}_{fg,1}{h}_{fg,1}={\dot{m}}_{fg,1}{h}_{fg,fet}\\ {\dot{m}}_{fg,2}{h}_{fg,2}-{\dot{m}}_{fg,1}{h}_{fg,1}={\dot{m}}_{ia,SH2}{h}_{ia}+{\dot{Q}}_{fg,SH2}-{\dot{m}}_{fa}{c}_{p,ash}\left(T_{fg,2}-T_{fg,1}\right)\\ {\dot{m}}_{fg,3}{h}_{fg,3}-{\dot{m}}_{fg,2}{h}_{fg,2}={\dot{m}}_{ia,CAV}{h}_{ia}+{\dot{Q}}_{fg,CAV}-{\dot{m}}_{fa}{c}_{p,ash}\left(T_{fg,3}-T_{fg,2}\right)\\ {\dot{m}}_{fg,4}{h}_{fg,4}-{\dot{m}}_{fg,3}{h}_{fg,3}={\dot{m}}_{ia,SH1}{h}_{ia}+{\dot{Q}}_{fg,SH1}-{\dot{m}}_{fa}{c}_{p,ash}\left(T_{fg,4}-T_{fg,3}\right)\\ {\dot{m}}_{fg,5}{h}_{fg,5}-{\dot{m}}_{fg,4}{h}_{fg,4}={\dot{m}}_{ia,EC}{h}_{ia}+{\dot{Q}}_{fg,EC}-{\dot{m}}_{fa}{c}_{p,ash}\left(T_{fg,5}-T_{fg,4}\right)\\ {\dot{m}}_{fg,6}{h}_{fg,6}-{\dot{m}}_{fg,5}{h}_{fg,5}={\dot{m}}_{ia,AH2}{h}_{ia}+{\dot{Q}}_{fg,AH2}-{\dot{m}}_{fa}{c}_{p,ash}\left(T_{fg,6}-T_{fg,5}\right)\\ {\dot{m}}_{fg,7}{h}_{fg,7}-{\dot{m}}_{fg,6}{h}_{fg,6}={\dot{m}}_{ia,AH1}{h}_{ia}+{\dot{Q}}_{fg,AH1}-{\dot{m}}_{fa}{c}_{p,ash}\left(T_{fg,7}-T_{fg,6}\right)\end{array}$$

(10)

In Equation (10), $h_{fg,j}$ is the enthalpy of the gas mixture (excluding the fly ash) at each node $j$, $h_{fg,aft}$ is the gas mixture enthalpy at the adiabatic flame temperature, and $h_{fg,fet}$ is the gas mixture enthalpy at the furnace exit temperature. $h_{fg,aft}$ will be calculated within the AFT component, and $h_{fg,fet}$ will be calculated within the FUR component. These calculations will be addressed later as part of the component characteristic equations. $h_{ia}$ is the enthalpy of the ingress air, ${\dot{Q}}_{fg}$ is the rate of heat transfer to the flue gas mixture in each component, $c_{p,ash}$ (J/kgK) is the ash specific heat capacity, and $T_{fg,j}$ (K) is the flue gas mixture temperature at each node $j$.

In Equation (10), the left-hand side again specifies the connectivity, while the right-hand side contains the boundary values, the rates of heat transfer to the gas mixture, and now also the energy sources/sinks in the gas mixture due to the cooling/heating of the entrained fly ash particles.

The flue gas mixture (excluding the fly ash particles) contains H₂O, CO₂, N₂, SO₂ and O₂. Therefore, besides the overall mass balances given in Equation (9), it is also necessary to write the mass balances for each of the species while taking into account the ingress air, which contains H₂O, N₂ and O₂. This can be written for each component as follows:

$$\begin{array}{l}{Y}_{fg,H2O,e}=\left({\dot{m}}_{fg,i}{Y}_{fg,H2O,i}+{\dot{m}}_{ia}{Y}_{ai,H2O}\right)/{\dot{m}}_{fg,e}\\ Y_{fg,CO2,e}={\dot{m}}_{fg,i}{Y}_{fg,CO2,i}/{\dot{m}}_{fg,e}\\ Y_{fg,N2,e}=\left({\dot{m}}_{fg,i}{Y}_{fg,N2,i}+{\dot{m}}_{ia}{Y}_{ai,N2}\right)/{\dot{m}}_{fg,e}\\ Y_{fg,SO2,e}={\dot{m}}_{fg,i}{Y}_{fg,SO2,i}/{\dot{m}}_{fg,e}\\ Y_{fg,O2,e}=\left({\dot{m}}_{fg,i}{Y}_{fg,O2,i}+{\dot{m}}_{ia}{Y}_{ai,O2}\right)/{\dot{m}}_{fg,i}\end{array}$$

(11)

with $Y_{fg,j,i}$ and $Y_{fg,j,e}$ representing the mass fraction of each species $j$ within the flue gas mixture at the inlet and outlet respectively, and $Y_{ia,j}$ representing the mass fraction of each species $j$ within the ingress air.

2.2. Adiabatic Flame Temperature Model

The adiabatic flame temperature ($T_{fg,aft}$) is defined as the maximum theoretically possible temperature of the flue gas mixture (flame ball) at the outlet of the combustion zone. Combustion is a chemical reaction during which the carbon, hydrogen and sulphur in the fuel are oxidised and energy is released. The composition of any fuel can be described by a unique combination of elemental carbon, hydrogen, oxygen, nitrogen and sulphur, together with moisture, and non-combustibles. Mass balance requires that the total mass of each element be conserved during the chemical reaction, which in turn implies that the number of moles of each element is conserved.

In the context of steady-state integrated boiler process models, we typically assume instantaneous (i.e., not accounting for the kinetics) and complete combustion, except in the case of solid fuels such as coal and biomass, where a fraction of unburned carbon is assumed. This unburned carbon is conveyed out of the boiler via the flue gas flow stream, but while assuming it has a negligible impact on the energy balance downstream of the combustion zone.

The mass fractions of C, H, O, N, S, moisture, and ash in the fuel specification are usually provided on an as-received (AR) basis with units kilogram per kilogram of fuel (kg/kg_f). However, for boiler process analysis, it is convenient to utilise the mass balance on an ash-free basis.

The total mass flow rate of ash is equal to

$${\dot{m}}_{ash}=Y_{ash}{\dot{m}}_{AR,f}$$

(12)

with ${\dot{m}}_{AR,f}$ representing the fuel mass flow rate on an as-received basis and $Y_{ash}$ representing the ash mass fraction. The fly ash flow rate that is carried along with the flue gas is assumed to be a fraction of the total ash flow rate, therefore

$${\dot{m}}_{fa}=f_{fa}{\dot{m}}_{ash}$$

(13)

with $f_{fa}$ representing the fraction of fly-ash. Remember that ${\dot{m}}_{fa}$ is an input to the source terms of Equation (10).

The bottom ash flow rate is simply given by

$${\dot{m}}_{ba}={\dot{m}}_{ash}-{\dot{m}}_{fa}$$

(14)

The fuel mass flow rate ${\dot{m}}_f$ on an ash-free basis is given by

$${\dot{m}}_f=\left(1-Y_{ash}\right){\dot{m}}_{AR,f}$$

(15)

The ash-free mass fraction $Y_{f,j}$ of each element $j$ in the fuel mixture (including the moisture $Y_{f,M}$) can be determined from

$$Y_{f,j}=\frac{Y_{AR,j}}{1-Y_{ash}}$$

(16)

where $Y_{AR,j}$ is the mass fraction of each element in the as-received fuel mixture (including the moisture $Y_{AR,M}$).

It is also convenient in the combustion model to include only the fraction of carbon that is participating in the combustion process in the value of $Y_{f,C}$, which then excludes the unburned carbon. It is therefore good practice to reduce the participating carbon mass fraction at the outset as follows:

$$Y_{f,C}=Y_{f,C}-Y_{UC}$$

(17)

where $Y_{UC}$ is the mass fraction of unburned carbon in the ash-free fuel mixture in kg/kg_f. Note that the fuel mass fractions with units kg/kg_f are normalised with the total ash-free fuel mass flow rate that includes the unburned carbon but which excludes the ash.

The mean molecular mass of the ash-free fuel mixture is given by

$$M_f={\left({\displaystyle \sum _j\frac{Y_{f,j}}{M_j}}\right)}^{-1}$$

(18)

with $M_j$ (kg/kmol) representing the molecular mass of each element $j$.

We now define $n_{fC}$, $n_{fH}$, $n_{fO}$, $n_{fN}$, $n_{fS}$ and $n_{fM}$ as the fuel flow-based mole fractions of C, H, O, N, S and moisture (M) in the fuel, respectively, with units of kilomole per kilomole of fuel (kmol/kmol_f). The fuel flow-based mole fraction for each element is determined from

$$n_{f,j}=Y_{f,j}\frac{M_f}{M_j}$$

(19)

The actual mass flow rate ${\dot{m}}_j$ of any element $j$ can be obtained directly from the fuel flow-based mole fraction as follows:

$${\dot{m}}_j=n_j\frac{M_j}{M_f}{\dot{m}}_f$$

(20)

The minimum oxygen needed for complete combustion is called the stoichiometric or theoretical oxygen. It is made up of the O₂ required for the combustion of the participating carbon, hydrogen and sulphur, minus the elemental oxygen provided within the fuel. Therefore, it is given by

$$n_{stO2}=n_{fC}+\frac{1}{4}{n}_{fH}-\frac{1}{2}{n}_{fO}+n_{fS}$$

(21)

Now, let the fuel flow-based mole fractions of H₂O, O₂ and N₂ in the combustion air be $n_{rH2O}$, $n_{rO2}$ and $n_{rN2}$, respectively, again with units of kmol/kmol_f. The combustion air will include a given fraction of excess air indicated by $\alpha$, where $\alpha =1.0$ implies zero excess air and $\alpha =1.2$ implies 20% excess air. Given this, the fuel flow-based mole fraction of O₂ in the combustion air for a given value of $\alpha$ will be

$$n_{rO2}=\alpha ·n_{stO2}$$

(22)

Since dry air can be approximated as a mixture consisting of 21% oxygen and 79% nitrogen on a volume/molar basis, the fuel flow-based mole fraction of N₂ in the combustion air will be equal to

$$n_{rN2}=3.7619\alpha ·n_{stO2}$$

(23)

The moisture content of the air is expressed in psychrometry terms as a ratio called the specific humidity (or humidity ratio), which is measured in kilogram of moisture per kilogram of dry air (kg/kg_a) and indicated with the symbol $w$. Therefore, the fuel flow-based mole fraction of H₂O in the combustion air will be equal to

$$n_{rH2O}=4.7619\alpha w\frac{M_{air}}{M_{H2O}}{n}_{stO2}$$

(24)

where $M_{air}=28.84\mathrm{kg}/\mathrm{kmol}$ is the mean molecular mass of air.

Based on Equations (20) to (24), the required total mass flow rate of humid combustion air is given by

$${\dot{m}}_{a,ca}=\left(n_{rO2}{M}_{O2}+n_{rN2}{M}_{N2}+n_{rH2O}{M}_{H2O}\right)\frac{{\dot{m}}_f}{M_f}$$

(25)

Remember that ${\dot{m}}_{a,ca}$ is the boundary value required in Equation (3).

Let the fuel flow-based mole fractions of CO₂, N₂, SO₂, H₂O and O₂ in the combustion products be $n_{pCO2}$, $n_{pN2}$, $n_{pSO2}$, $n_{pH2O}$ and $n_{pO2}$, respectively, with units of kmol/kmol_f. All the participating carbon in the fuel is oxidised to form CO₂ and therefore

$$n_{pCO2}=n_{fC}$$

(26)

Some of the nitrogen in the fuel and combustion air forms various NOx variants in the products. However, in the current context, it is assumed that all the nitrogen is converted to stable diatomic nitrogen in the products. Therefore

$$n_{pN2}=\frac{1}{2}{n}_{fN}+n_{rN2}$$

(27)

The first term on the right-hand side of Equation (27) represents the elemental nitrogen in the fuel, and the second terms represents the nitrogen introduced via the combustion air.

We also assume that all the elemental sulphur in the fuel burns to SO₂, therefore

$$n_{pSO2}=n_{fS}$$

(28)

The flue gas mixture will include the moisture introduced via the fuel, the water formed during the oxidation of the hydrogen, and the additional moisture introduced via the combustion air. Therefore

$$n_{pH2O}=n_{fM}+\frac{1}{2}{n}_{fH}+n_{rH2O}$$

(29)

The additional oxygen that enters via the excess air does not take part in combustion and therefore forms part of the flue gas mixture as follows:

$$n_{pO2}=\left(\alpha -1\right)n_{stO2}$$

(30)

The total mass flow rate of the gaseous combustion products leaving the combustion zone can be determined from

$${\dot{m}}_{fg,aft}=\left(n_{pCO2}{M}_{CO2}+n_{pN2}{M}_{N2}+n_{pSO2}{M}_{SO2}+n_{pH2O}{M}_{H2O}+n_{pO2}{M}_{O2}\right)\frac{{\dot{m}}_f}{M_f}$$

(31)

Remember that ${\dot{m}}_{fg,aft}$ is the boundary value required in Equation (9) and that a constant flow of fly ash at a rate of ${\dot{m}}_{fa}$ will be carried along with it.

Based on Equations (21) to (30), the fuel flow-based molar balance in kmol/kmol_f can conveniently be written as

$$\begin{array}{l}{C}_{n_{fC}}{H}_{n_{fH}}{O}_{n_{fO}}{N}_{n_{fN}}{S}_{n_{fS}}{M}_{n_{fM}}+n_{rH2O}{H}_{2}O+n_{rO2}{O}_2+n_{rN2}{N}_2\\ \to n_{pCO2}C{O}_2+n_{pN2}{N}_2+n_{pSO2}S{O}_2+n_{pH2O}{H}_{2}O+n_{pO2}{O}_2\end{array}$$

(32)

Equation (20) can be employed to convert each term in Equation (32) to the mass flow rate of the individual species.

The energy balance equation for the combustion zone is given by

$${\dot{m}}_{fg,aft}{h}_{fg,aft}={\dot{m}}_f{h}_f+{\dot{m}}_{a,ca}{h}_{a,ca}+{\dot{m}}_f\left(HH{V}_f-Y_{UC}HH{V}_C\right)-{\dot{m}}_{fa}{h}_{ash,aft}-{\dot{m}}_{ba}{h}_{ash,bat}$$

(33)

$HH{V}_f$ (J/kg) is the higher heating value of the fuel on an ash-free basis, $HH{V}_C=32763\mathrm{J}/\mathrm{kg}$ is the higher heating value of carbon, and $h_f$ is the sensible enthalpy of the fuel entering the combustion zone. As stated before, $h_{fg,aft}$ is the gas mixture enthalpy at $T_{fg,aft}$, which can be determined from

$$h_{fg,aft}={\varphi }_{fg}\left(Y_{fg},T_{fg,aft},p_{fg,aft}\right)$$

(34)

with ${\varphi }_{fg}$ representing the appropriate fluid property relationship for the flue gas mixture, which accounts for the local flue gas species mass fractions as well as the local temperature and pressure. $h_{ash,aft}$ is the sensible enthalpy of the fly-ash particles at the adiabatic flame temperature, and $h_{ash,bat}$ is the sensible enthalpy of the bottom ash particles at the bottom ash temperature ($T_{bat}$). Therefore

$$h_{ash,aft}={\varphi }_{ash}\left(T_{fg,aft}\right)$$

(35)

and

$$h_{ash,bat}={\varphi }_{ash}\left(T_{bat}\right)$$

(36)

with ${\varphi }_{ash}$ representing the sensible enthalpy versus temperature relationship for the solid fly ash particles. Note that a suitable value must be assumed for $T_{bat}$. For the fuel, we have

$$h_f={\varphi }_f\left(T_f\right)$$

(37)

with ${\varphi }_f$ representing the sensible enthalpy versus temperature relationship for the solid fuel particles, and $T_f$ representing the temperature of the fuel at the inlet to the combustion zone.

It is important to note that Equation (33) utilises the HHV of the fuel as opposed to the LHV (lower heating value), which is often employed in conventional combustion calculations. This difference will be addressed in the next section, which explains the fluid property relations.

Since $T_{fg,aft}$ is not known beforehand, an iterative approach is required to solve simultaneously for Equations (33) to (35), which then yields the values of $T_{fg,aft}$, $h_{fg,aft}$ and $h_{ash,aft}$. Remember that $h_{fg,aft}$ is a boundary value required for Equation (10).

This brings us to the point where we need to consider what the appropriate property relationships would be for ${\varphi }_{fg}$, ${\varphi }_{ash}$ and ${\varphi }_f$.

2.3. Fluid and Solid Particle Property Relationships

Values of enthalpy are always provided relative to a given reference state $h_{ref}$ where ${h_{ref}|}_{p_{ref}}^{T_{ref}}=0$ with $T_{ref}$ and $p_{ref}$ representing the chosen reference temperature and pressure, respectively. Therefore, the values of enthalpy that are commonly used in the energy balance equations are really enthalpy differences, i.e., $\Delta h=h-h_{ref}$. In combustion analysis, the composition of the fluids at the outlet of the process is no longer the same as at the inlet of the process. This implies that there must be a common reference state for all the substances taking part, which includes the combustion air, the fuel, the ash, and the flue gas mixture. This chosen common reference state is at 25 °C and 101.325 kPa, which is known as the standard reference state. Therefore, for the solid particles

$${\varphi }_{ash}\left(T\right)=c_{p,ash}\left(T-25°\mathrm{C}\right)$$

(38)

and

$${\varphi }_f\left(T\right)=c_{p,f}\left(T-25°\mathrm{C}\right)$$

(39)

where $c_{p,ash}$ and $c_{p,f}$ are the specific heat capacities of the ash and fuel, respectively, which are usually assumed to be constant throughout the boiler.

The combustion and ingress air enthalpies $h_a$, $h_{a,ca}$ and $h_{ia}$ as well as the flue gas enthalpies $h_{fg}$, $h_{fg,aft}$ and $h_{fg,fet}$ are all expressed relative to the standard reference state. The enthalpy of the H₂O vapour in the gas mixture is therefore assumed to be zero at 25 °C and 101.325 kPa. However, at the reference state, the H₂O would be a compressed liquid rather than vapour. This means that the enthalpy of the vapour in the gas mixture should include the latent heat of vaporisation at the reference temperature, i.e., ${h_{fg}|}_{25°\mathrm{C}}\approx 2442.3\mathrm{kJ}/\mathrm{kg}$. In the conventional textbook approach, combustion analysis is usually based on ideal gas enthalpies where the latent heat of vaporisation is not included in the enthalpy of the vapour. It is therefore assumed that the enthalpy of the water vapour is equal to zero at 25 °C. This choice of using the ideal gas approach is rational, since boiler analysis usually involves gas mixtures at low pressures and high temperatures. However, the ideal gas approach necessitates the use of the LHV in the energy balance equation instead of the HHV.

The fact that we have used the HHV in Equation (33) implies that we are using real gas properties for the flue gas mixture rather than ideal gas properties. It can easily be shown that we can replace HHV with LHV provided that the real gas enthalpy $h_{fg,aft}$ is also replaced by its equivalent ideal gas value $h_{fg,aft,ideal}=h_{fg,aft}-Y_{fg,H2O}{h_{fg}|}_{25°\mathrm{C}}$. We accomplish this by replacing HHV in Equation (33) with its definition in terms of LHV, as follows:

$$\begin{array}{c}{\dot{m}}_{fg,aft}{h}_{fg,aft}={\dot{m}}_f{h}_f+{\dot{m}}_{a,ca}{h}_{a,ca}+{\dot{m}}_f\left(HH{V}_f-Y_{UC}HH{V}_C\right)-{\dot{m}}_{fa}{h}_{ash,aft}-{\dot{m}}_{ba}{h}_{ash,bat}\\ {\dot{m}}_{fg,aft}{h}_{fg,aft}={\dot{m}}_f{h}_f+{\dot{m}}_{a,ca}{h}_{a,ca}+{\dot{m}}_f\left(\begin{array}{l}\left(LH{V}_f+\frac{{\dot{m}}_{fg,H2O}}{{\dot{m}}_f}{h_{fg}|}_{25°\mathrm{C}}\right)\\ -Y_{UC}HH{V}_C\end{array}\right)-{\dot{m}}_{fa}{h}_{ash,aft}-{\dot{m}}_{ba}{h}_{ash,bat}\\ {\dot{m}}_{fg,aft}{h}_{fg,aft}-{\dot{m}}_{fg,H2O}{h_{fg}|}_{25°\mathrm{C}}={\dot{m}}_f{h}_f+{\dot{m}}_{a,ca}{h}_{a,ca}+{\dot{m}}_f\left(LH{V}_f-Y_{UC}HH{V}_C\right)-{\dot{m}}_{fa}{h}_{ash,aft}-{\dot{m}}_{ba}{h}_{ash,bat}\\ {\dot{m}}_{fg,aft}{h}_{fg,aft,ideal}={\dot{m}}_f{h}_f+{\dot{m}}_{a,ca}{h}_{a,ca}+{\dot{m}}_f\left(LH{V}_f-Y_{UC}HH{V}_C\right)-{\dot{m}}_{fa}{h}_{ash,aft}-{\dot{m}}_{ba}{h}_{ash,bat}\end{array}$$

(40)

One implication of using real gas properties for the gas mixture is that in the application of the mixing law, the real gas enthalpy of each of the constituents must be evaluated at the mixture temperature and partial pressure of the constituent rather than at the mixture pressure. This is consistent with reality. The real gas enthalpy of the constituent at the reference state is then subtracted from this to obtain the value that is employed in the mixing law. The benefit of this approach is that it will provide consistent results even if the gas mixture temperature falls below the dew point temperature, resulting in condensation of the water vapour. It is important to account for this in the analysis of condensing boilers, and condensation should be avoided in non-condensing boilers.

The mixing law applied for enthalpy and specific heat is usually a simple mass-weighted average, while for viscosity and conductivity, the well-known Wilke method is preferred. The density of the mixture is determined using the ideal gas law with an effective mixture specific gas constant defined as $R_m={\displaystyle {\sum }_j{Y}_j{R}_j}$, where $R$ (J/kgK) is the specific gas constant and $j$ indicates the specific constituent in the mixture.

For the steam network, the enthalpies are also determined using real gas properties. However, it is not necessary to scale the values relative to the standard reference state because the steam does not take part in the combustion process. The reference state for the standard steam properties is usually the triple point of water, i.e., 0.01 °C and 0.6113 kPa.

2.4. Furnace Model

The energy balance equation for the gas and particle suspension in the furnace can be written similar to Equation (10) as:

$${\dot{Q}}_{fg,FUR}={\dot{m}}_{fg,aft}\left(h_{fg,aft}-h_{fg,fet}\right)+{\dot{m}}_{fa}{c}_{p,ash}\left(T_{fg,aft}-T_{fg,fet}\right)$$

(41)

with ${\dot{Q}}_{fg,FUR}$ representing the total rate of heat transfer from the flue gas mixture and

$$h_{fg,fet}={\varphi }_{fg}\left(Y_{fg},T_{fg,fet},p_{fg,fet}\right)$$

(42)

Note that for the furnace model, the rate of heat transfer from the gas is traditionally written as a positive value as opposed to the conventional definition used in Equation (10) where heat transfer from the gas would be a negative value. The traditional convention of positive heat transfer in the furnace is also adopted here.

The so-called projected (Gurvich/Blokh) and direct (Hottel) methods are semi-empirical zero-dimensional models of furnace heat transfer in which all the physical quantities are assumed to be uniform, and the results are averaged over the heating surfaces. It includes various empirical parameters that were originally derived for specific boiler layouts and a limited range of operating conditions. Despite this, it is extensively applied in industry for the design and analysis of all kinds of boilers over a wide range of operating conditions.

These zero-dimensional models assume that the heat transfer in the furnace is dominated by radiation and that the furnace contains a single homogeneous participating medium radiating to a single isothermal absorbing, reflecting, and emitting the surrounding surface, even though it may consist of several different surfaces at different temperatures. A single equation is therefore employed to calculate the total rate of radiant heat (${\dot{Q}}_{rad}$) projected from the gas volume to the surrounding surfaces. Following this, appropriate portions of the total radiation are allocated to the different surfaces, which could include the water walls, refractory, and the virtual furnace exit plane. The water flow path through the water walls and the gas flow path through the furnace volume are again assumed to be one-dimensional with the balance equations applied between the outlet and inlet, as described in Section 2.1.2 and Section 2.1.3. The difference between the projected method and the direct method is in the calculation of ${\dot{Q}}_{rad}$.

2.4.1. Projected Method

The projected method essentially assumes that the surrounding wall temperature is low enough so that emission from the wall can be neglected. Therefore, the total rate of radiant heat projected from the gas volume to the surrounding surfaces is given by [44]

$${\dot{Q}}_{rad}=\psi {\epsilon }_f\sigma A_{rad}{T}_g^4$$

(43)

with $T_g$ representing the effective gas temperature, $A_{rad}$ (m²) representing the total surrounding projected area, $\sigma$ (W/m²K⁴) representing the Stefan-Boltzmann constant, and ${\epsilon }_f$ representing the furnace effective emissivity. $\psi$ is an area-weighted furnace efficiency factor that accounts for the effects of fouling and angular coefficient (view factor) of the different surfaces that make up the surrounding surface, including the virtual furnace exit plane. The efficiency factor is defined as the ratio of the heat flux actually absorbed by the wall to the radiation heat flux incident on the wall (i.e., the irradiation). This area-weighted efficiency factor is found from

$$\psi A_{rad}={\displaystyle \sum _i{\psi }_i{A}_{rad,i}}={\displaystyle \sum _i{x}_i{\xi }_i{A}_{rad,i}}$$

(44)

where $x_i$, ${\xi }_i$ and $A_{rad,i}$ are the angular coefficient, wall fouling factor and projected radiation surface of each surface $i$. The angular coefficient value used for membrane-type water walls is a function of the tube spacing and diameter and is typically around 0.9–1.0. For water walls covered by refractory, it is 1.0, and for the furnace exit plane, it is also 1.0. The wall fouling factor used for gas and oil-fired furnaces is 0.65; for heavy oil, it is 0.55; for coal, it is 0.4; for grate fuels, it is 0.6; for other fuels, it is 0.2; for refractory, it is 0.1; and for the furnace exit, it is 1.0.

The furnace effective emissivity is given by

$${\epsilon }_f={\left(1+\psi \left(\frac{1}{{\epsilon }_{gp}}-1\right)\right)}^{-1}$$

(45)

with ${\epsilon }_{gp}$ representing the effective emissivity of the flame ball, i.e., the hot flue gas mixture and particle suspension, at the gas temperature $T_g$ (the calculation of ${\epsilon }_{gp}$ will be addressed later). Equation (45) can be derived from fundamental principles by approximating the furnace wall and the gas as two infinitely large parallel surfaces. Then, we write the energy balance equations for the radiation heat transfer network between the two in terms of emissive power, irradiation and radiosity while including the emissivity and reflectivity of both sides. The ratio of heat flux absorbed by the wall to the radiation heat flux incident on the wall (which is equal to the radiosity of the gas) is replaced by $\psi$. Finally, the radiosity of the gas is written as ${\epsilon }_f\sigma T_g^4$, i.e., the effective rate at which radiant energy leaves the gas, which then leads to Equation (45).

The heat losses from the furnace are accounted for via the heat preservation coefficient $\phi$ as follows:

$${\dot{Q}}_{rad}=\phi {\dot{Q}}_{fg,FUR}$$

(46)

Since $\phi <1$, it implies that the heat transfer from the flue gas to the surrounding surfaces is less than the heat extracted from the flue gas, thereby accounting for the losses from the furnace.

The energy balance equation for the furnace can also be written as

$${\dot{Q}}_{fg.FUR}={\dot{m}}_f\overline{vc}\left(T_{fg,aft}-T_{fg,fet}\right)$$

(47)

with $\overline{vc}$ (J/kg_fK) defined as the mean overall heat capacity based on the fuel flow rate. Based on Equations (41) and (47), $\overline{vc}$ can be determined from

$$\overline{vc}=\frac{{\dot{m}}_{fg,aft}\left(h_{fg,aft}-h_{fg,fet}\right)+{\dot{m}}_{fa}{c}_{p,ash}\left(T_{fg,aft}-T_{fg,fet}\right)}{{\dot{m}}_f\left(T_{fg,aft}-T_{fg,fet}\right)}$$

(48)

Combining Equations (43), (46) and (47) yields

$$\phi {\dot{m}}_f\overline{vc}\left(T_{fg,aft}-T_{fg,fet}\right)=\psi {\epsilon }_f\sigma A_{rad}{T}_g^4$$

(49)

which provides the relationship between $T_{fg,aft}$, $T_{fg,fet}$ and $T_g$. Assuming that $T_{fg,aft}$ is known from the solution of the AFT model, an additional empirical closure relationship is required to solve for $T_{fg,fet}$ and $T_g$. The projected method provides the following empirical equation, which is also provided in references [20,44,45]:

$$\frac{T_{fg,fet}}{T_{fg,aft}}={\left(M{\left(\frac{{\epsilon }_f}{Bo}\right)}^{0.6}+1\right)}^{-1}$$

(50)

$M$ is the flame centre modification factor (or M-factor), and $Bo$ is the Boltzmann number defined as

$$Bo=\frac{\phi {\dot{m}}_f\overline{vc}}{\psi \sigma A_{rad}{T}_{fg,aft}^3}=\frac{T_g}{T_{fg,aft}-T_{fg,fet}}$$

(51)

The M-factor is derived from the empirical relationship

$$M=A-B\left(X_r+\Delta X\right)$$

(52)

where $X_r$ is the normalised burner height, $A$ and $B$ are empirical constants depending on the kind of fuel, and $\Delta X$ is a correction factor for the actual flame position, which depends on the burner type and tilt.

Knowing $T_{fg,aft}$ and $h_{fg,aft}$ from the solution of the AFT model, the simultaneous solution of Equations (48), (51), (50) and (42) yields the values of $\overline{vc}$, $Bo$, $T_{fg,fet}$ and $h_{fg,fet}$. Following this, ${\dot{Q}}_{fg,FUR}$ can be determined from Equation (41) or (47), and ${\dot{Q}}_{rad}$ can be determined from Equation (46).

Following this, the radiant heat transfer to each individual surface $i$, excluding the water walls, can be determined from

$${\dot{Q}}_{rad,i\ne ww}={\eta }_i\frac{{\psi }_i{A}_{rad,i}}{\psi A_{rad}}{\dot{Q}}_{rad}$$

(53)

where ${\eta }_i$ is a non-uniformity factor that depends on the assumed heat flux distribution along the height of the furnace and the position of surface $i$. Note that a low value of ${\psi }_i$ will result in a small portion of the heat transfer being allocated to surface $i$. Therefore, for refractory covered walls, we typically have $\xi =0.1$ [46], which is equivalent to specifying a large thermal resistance, resulting in a high surface temperature and low heat absorption rate. The value of $\xi$ therefore effectively compensates for the fact that the projected method does not account for the surrounding wall temperature in Equation (43).

For the direct radiation onto the furnace exit plane, an additional re-radiation factor $\beta$ (or furnace–platen superheater interaction factor) is also included as follows

$${\dot{Q}}_{rad,fe}=\beta {\eta }_{fe}\frac{{\psi }_{fe}{A}_{rad,fe}}{\psi A_{rad}}{\dot{Q}}_{rad}$$

(54)

where $\beta$ is a function of $T_{fg,fet}$ and the type of fuel, and its value varies between 1.0 and 0.5, and $A_{rad,fe}$ is the area of the furnace exit plane. An attractive attribute of the projected method is that it allows for the calculation of this direct radiation onto the furnace exit plane, which is then assumed to impinge onto the heat exchangers situated at the furnace exit without knowing the details of the downstream heat exchanger.

Finally, to ensure an overall energy balance for the furnace, the remainder of the total radiant heat is allocated to the water wall, i.e.,

$${\dot{Q}}_{rad,ww}={\dot{Q}}_{rad}-{\displaystyle \sum _{i\ne ww}{\dot{Q}}_{rad,i}}$$

(55)

Combining Equations (49) and (50) also provides the relationship for $T_g$ as a function of $T_{fg,fet}$ and $T_{fg,aft}$ as follows:

$${\left(\frac{T_g}{T_{fg,aft}}\right)}^4={\left(M\frac{T_{fg,fet}}{T_{fg,aft}}\right)}^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}{\left(1-\frac{T_{fg,fet}}{T_{fg,aft}}\right)}^{-\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$$

(56)

or $T_g$ can simply be determined implicitly from Equation (43) or (49).

2.4.2. Projected Method for Large Boilers

According to Zhang et al. [44], the projected method described above was developed based on experimental data from a boiler with output of 200–300 t/h (160–240 MW), and the influence of non-uniformity in the furnace was ignored during the data analysis and modelling. Therefore, it should not provide very accurate results for larger boilers in the sense that the calculated furnace outlet temperature is typically lower than the measured values. This in effect implies that the total heat transfer is overestimated. Therefore, for larger boilers, they provide an alternative empirical formulation to Equation (50), namely:

$$\frac{T_{fg,fet}}{T_{fg,aft}}=1-M{\left(\frac{{\epsilon }_f\psi T_{fg,aft}^2}{10800·q_H}\right)}^{0.6}$$

(57)

where $q_H$ is the heat release rate per unit of radiative heating surface area in (kW/m²), which can be approximated as

$$q_H\approx {\dot{m}}_f\left(HH{V}_f-Y_{UC}HH{V}_C\right)/A_{rad}$$

(58)

2.4.3. Direct Method

In contrast to the projected method of Equation (43), the direct method accounts for the temperature and emissivity of the surface surrounding the furnace. It states that [44,47]

$${\dot{Q}}_{rad}=\frac{{\epsilon }_w}{{\alpha }_{gp}+{\epsilon }_w-{\alpha }_{gp}{\epsilon }_w}\sigma A_{rad}\left({\epsilon }_{gp}{T}_g^4-{\alpha }_{gp}{T}_w^4\right)$$

(59)

with $T_w$ representing the effective surrounding wall/surface temperature, ${\epsilon }_w$ representing the wall emissivity and ${\alpha }_{gp}$ representing the effective gas and particle suspension absorptivity determined at the wall temperature $T_w$ (the calculation of ${\epsilon }_{gp}$ and ${\alpha }_{gp}$ will be addressed later). The value of the wall emissivity is usually around 0.8–0.85.

It is important to understand where Equation (59) originated. It can be derived by writing the energy balance equations for radiation heat transfer between multiple diffuse-grey surfaces surrounding a single participating medium (gas). It accounts for the view factors between the various surfaces, the emissivity of the walls, and the absorption and emission due to the participating medium, but it neglects scattering. This results in a set of equations representing a radiation heat transfer network, which can be solved simultaneously to obtain the net rate of heat transfer from the gas and from each of the surfaces, as a function of the effective gas temperature $T_g$ and the various wall temperatures $T_{w,i}$. If it is now assumed that there is only a single surrounding surface with effective temperate $T_w$, we arrive at Equation (59).

The guideline provided by [47] regarding the representative gas temperature $T_g$ is that it should be taken as equal to the furnace exit temperature $T_{fg,fet}$. (This may not be completely consistent with the derivation of the equation from fundamental principles, since the weighted average temperature of the gas volume will be higher than the furnace exit temperature).

The direct radiation leaving the furnace exit plane and impinging onto the heat exchangers situated behind the furnace exit can then be determined from

$${\dot{Q}}_{rad,fe}=\frac{{\epsilon }_{w,hx}}{{\alpha }_{gp,hx}+{\epsilon }_{w,hx}-{\alpha }_{gp,hx}{\epsilon }_{w,hx}}\sigma A_{rad,fe}\left({\epsilon }_{gp}{T}_g^4-{\alpha }_{gp,hx}{T}_{w,hx}^4\right)$$

(60)

where $T_{w,hx}$ represents the heat exchanger surface temperature, ${\epsilon }_{w,hx}$ represents the heat exchanger wall emissivity and ${\alpha }_{gp,hx}$ represents the gas absorptivity determined at the heat exchanger wall temperature $T_{w,hx}$. This implies that some detail of the downstream heat exchanger must be known to solve the furnace model. Note that the direct radiation eventually absorbed by the downstream heat exchanger will depend on the view factor between the virtual furnace exit plane and the tube bundle, which will be addressed later.

The direct method accounts for the heat losses from the furnace via a radiation loss factor ($f_{radloss}<1$) as follows:

$${\dot{Q}}_{rad}+{\dot{Q}}_{rad,fe}=\left(1-f_{radloss}\right){\dot{Q}}_{fg,FUR}$$

(61)

One benefit of accounting for the wall temperature in Equation (59) is that one can now explicitly include the conductive/convective thermal resistances of the different surfaces, such as the water wall and refractory that may cover portions of the water wall. This means that the heat fluxes through the different surfaces to the water/steam can be estimated. The energy balance at the surrounding wall is given by

$${\dot{Q}}_{rad}=\left(T_w-T_{sat}\right)/R_w$$

(62)

with $T_{sat}$ representing the saturation temperature of the steam within the water wall tubes. $R_w$ (K/W) is the effective wall thermal resistance determined from $R_w={\left({\displaystyle \sum _{i}1/R_{w,i}}\right)}^{-1}$, with $R_{w,i}$ representing the thermal resistance of the physical wall associated with each surface. (In this regard, the convective heat transfer coefficient for the evaporating water/steam mixture within the water wall is often simply assumed to be very large).

Knowing $T_{fg,aft}$ from the solution of the AFT model, simultaneous solution of Equations (41), (42), (59) and (60)–(62) yield the values of $T_g=T_{fg,fet}$, $h_{fg,fet}$, ${\dot{Q}}_{fg,FUR}$, ${\dot{Q}}_{rad}$, ${\dot{Q}}_{rad,fe}$ and $T_w$.

Following this, the heat transfer to each individual surface $i$, excluding the furnace exit, can be determined from

$${\dot{Q}}_{rad,i\ne fe}=\left(T_w-T_{sat}\right)/R_{w,i}$$

(63)

2.5. Thermal Radiation of Gas–Solid Suspension

There are three possible approaches to the calculation of the effective emissivity and absorptivity of the flue gas and particle suspension, which will be referred to as the ‘standard model’, the ‘low particle load model’, and the ‘high particle load model’.

In all the models, the mean beam length ($S$) characterises the geometry of the heat transfer surface involved. For the furnace, it is determined as

$$S_{fur}=\frac{3.6V_{fur}}{A_{fur}}$$

(64)

with $V_{fur}$ (m³) representing the furnace volume and $A_{fur}$ (m²) representing the total surface area surrounding the furnace volume, which includes the virtual furnace exit area.

2.5.1. Standard Model

The standard model is presented by Zhang [44] as the former Soviet Union’s proposed 1973 standard for boiler thermal calculation, as well as by Lin [45] and Basu [46]. It only addresses the calculation of the emissivity of the flue gas and particle suspension and not the absorptivity.

The effective emissivity of the flue gas and particle suspension is given by

$${\epsilon }_{gp}=1-e^{-k_{fur}{p}_{fur}{S}_{fur}}$$

(65)

with $k_{fur}$ representing the effective furnace pressure-based absorption/emission coefficient with units (m∙MPa)⁻¹ and $p_{fur}$ representing the furnace pressure in MPa. (Note that the absorption coefficient $\kappa$ with units m⁻¹ is defined in terms of the pressure-based absorption coefficient as $\kappa =k\times p$).

The effective furnace pressure-based absorption coefficient is given by

$$k_{fur}=k_g\left(X_{fg,H2O}+X_{fg,CO2}+X_{fg,SO2}\right)+k_{fa}{\mu }_{fa}+k_{coke}{x}_1{x}_2$$

(66)

with $k_g$ representing the effective pressure-based absorption coefficient of the triatomic gasses in the flue gas mixture with units (m∙MPa)⁻¹, $X_{fg,H2O}$, $X_{fg,CO2}$ and $X_{fg,SO2}$ are the mole fractions of H₂O, CO₂ and SO₂, respectively, in the flue gas mixture, $k_{fa}$ representing the fly ash pressure-based absorption coefficient with units (m∙MPa)⁻¹, ${\mu }_{fa}={\dot{m}}_{fa}/\left({\dot{m}}_{fg}+{\dot{m}}_{fa}\right)$ representing the mass fraction of the fly ash in the flue gas and particle suspension, and $k_{coke}=10.2$ (m∙MPa)⁻¹ representing the pressure-based absorption coefficient of coke particles. $x_1$ is an empirical constant associated with the kind of coal, which is equal to 1.0 for anthracite coal with low volatile matter and 0.5 for other coals with high volatile matter. $x_2$ is an empirical constant accounting for the influence of the burning equipment, which is equal to 0.03 for grate firing and 0.1 for suspension firing.

The pressure-based absorption coefficient of the triatomic gasses at the gas temperature $T_g$ is given by Zhang [44], Lin [45] and Basu [46] as

$$k_g=\left(\frac{7.8+16X_{fg,H2O}}{3.16\sqrt{\left(X_{fg,H2O}+X_{fg,CO2}+X_{fg,SO2}\right)p_{fur}{S}_{fur}}}-1\right)\left(1-0.37\frac{T_g}{1000}\right)$$

(67)

with the pressure and temperature given in MPa and Kelvin, respectively.

The fly ash pressure-based absorption coefficient at the gas temperature is given by Zhang [44] as

$$k_{fa}=\frac{48350{\rho }_{fg,n}}{{\left(T_g^2{d}_p^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}$$

(68)

with ${\rho }_{fg,n}$ representing the flue gas density in kg/m³ at normal conditions (101.325 kPa, 0 °C) and $d_p$ representing the average diameter of the fly ash particles in μm. For tubular ball mills, $d_p$ is typically 13 μm, for hammer mills, it is 16 μm, and for grate firing, it is 20 μm.

With regard to Equation (68), it is important to note that we have presented the version applied in the example calculation provided in Zhang [44]. In Chapter 2 of the same reference [44], the equation is given as $k_{fa}=\frac{43000{\rho }_{fg,n}}{{\left(T_g^2{d}_p^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}$ (m∙MPa)⁻¹, which will provide a value for $k_{fa}$ that is about 89% of the value obtained with Equation (68). Furthermore, in Basu [46], the equation is given as $k_{fa}=\frac{5990}{{\left(T_g^2{d}_p^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}$ (m∙MPa)⁻¹, which was taken directly from Lin [45].

Using this would result in a value of $k_{fa}$ that is about 9.2% of the value obtained when using Equation (68). Having these different versions of what is supposedly the same equation creates uncertainty about the standard method for calculation of the fly ash pressure-based absorption coefficient. This will be addressed again in the case studies presented in Section 3.

2.5.2. Low Particle Load Model

The emissivity of the gas mixture (excluding the particles) is given by

$${\epsilon }_g=1-e^{-{\kappa }_{\epsilon ,g}{S}_{fur}}$$

(69)

with ${\kappa }_{\epsilon ,g}$ representing the emission coefficient of the gas mixture with unit m⁻¹ (as opposed to the pressure-based coefficient $k$ with units (m∙MPa)⁻¹).

The low particle load model is presented by Brummel [48] and employs the weighted sum of grey gases model for the emissivity of the flue gas mixture, which was originally presented by Hottel [49] and given in Vortmeyer [50] while referencing Johnson and Beér [51]. For a mixture of H₂O and CO₂ at a total pressure of $p=1.0bar$, $\frac{p_{H2O}}{p_{CO2}}=1$ and $1100K<T_g<1800K$ and $0.2m<S<6m$, the emissivity of the flue gas mixture at the gas temperature $T_g$ is determined as follows:

$${\epsilon }_g\left(T_g\right)={\displaystyle \sum _{i=1}^3\left[\left(b_{1,i}+b_{2,i}\frac{T_g}{1000K}\right)\left(1-e^{-k_i\left(X_{fg,H2O}+X_{fg,CO2}\right)pS}\right)\right]}$$

(70)

with the pressure given in bar and the coefficients provided in

Table 1. Coefficients for Equation (70).

i	b1,i	b2,i	ki (m∙bar)−1
1	0.130	0.265	0.0
2	0.595	−0.150	0.824
3	0.275	−0.115	25.91

.

According to Vortmeyer [50], the corresponding absorptivity ${\alpha }_g$ of the gas mixture at the wall temperature $T_w$ can also be calculated using Equation (70) and Table 1 by replacing the gas temperature with the wall temperature.

Brummel [48] presents the low particle load model while referencing Biermann [52]. Here, the emissivity of the particle cloud is determined via a simple absorption model while neglecting any scattering of radiation by the dispersed particles. This simplification is only valid for low particle loads given by $L_p=\frac{{\dot{m}}_{fa}}{{\dot{m}}_{fg}/{\rho }_{fg}}<0.005\frac{k{g}_{fa}}{m_{fg}^3}$. The effective emissivity of the particles is then given by

$${\epsilon }_p=1-e^{-Q_{abs}{A}_p{L}_p{S}_{fur}}$$

(71)

with $A_p=\frac{3}{2}\frac{1}{{\rho }_p{10}^{-6}{d}_p}$ (m²) the projected area of the fly ash particles per unit mass ($d_p$ has units μm). $Q_{abs}$ is the non-dimensional particle absorption efficiency, which will be addressed later.

The effective emissivity of the flue gas and particle suspension at the gas temperature is given by

$${\epsilon }_{gp}\left(T_g\right)={\epsilon }_g\left(T_g\right)+{\epsilon }_p-{\epsilon }_g\left(T_g\right)·{\epsilon }_p$$

(72)

Similarly, the effective absorptivity of the flue gas and particle suspension at the wall temperature is given by

$${\alpha }_{gp}\left(T_w\right)={\alpha }_g\left(T_w\right)+{\epsilon }_p-{\alpha }_g\left(T_w\right)·{\epsilon }_p$$

(73)

2.5.3. High Particle Load Model

The high particle load model is also presented by Brummel [48]. It employs the same weighted sum of grey gases model for the emissivity of the flue gas mixture as before. Having determined ${\epsilon }_g$ from Equation (70), the emission coefficient of the gas mixture at $T_g$ can be determined from:

$${\kappa }_{\epsilon ,g}\left(T_g\right)=-\frac{\ln \left(1-{\epsilon }_g\left(T_g\right)\right)}{S_{fur}}$$

(74)

According to Brummel [48], the low-particle load model presented before tends to overestimate the heat transfer from clouds of small ash particles with higher loads, and the root cause for this is the scattering of the radiation by the particles. The total scattering is divided into forward and backward scattering (backscattering). While forward scattering is assumed to have no impact on the radiation, the backscattering attenuates the emissivity by redirecting part of the radiation in the opposite direction. This backscattering becomes significant when the wavelength of the radiation and the size of the particles are of the same magnitude.

In the high particle load model, the effective emissivity of the flue gas and particle suspension is given by [48]

$${\epsilon }_{gp}\left(T_g\right)=\left(1-\beta \right)\left(\frac{1-e^{-{\mathsf{\Phi }}_{\epsilon ,gp}\left(T_g\right)}}{1+\beta e^{-{\mathsf{\Phi }}_{\epsilon ,gp}\left(T_g\right)}}\right)$$

(75)

with $\beta =\frac{\gamma -1}{\gamma +1}$ and $\gamma =\sqrt{1+\frac{2Q_{bsc}}{Q_{abs}}}$. ${\mathsf{\Phi }}_{\epsilon ,gp}$ is the optical thickness for the gas particle mixture which incorporates the emission coefficient of the gas mixture from Equation (74) and is given by

$${\mathsf{\Phi }}_{\epsilon ,gp}\left(T_g\right)=\left({\kappa }_{\epsilon ,g}\left(T_g\right)+Q_{abs}{A}_p{L}_p+k_{coke}{x}_1{x}_2{p}_{fur}\right)S_{fur}\gamma$$

(76)

$Q_{abs}$ is the non-dimensional particle absorption efficiency (also encountered in Equation (71)) and $Q_{bsc}$ is the non-dimensional particle scattering efficiency. Brummel [48] provides graphs for $Q_{abs}$ and $Q_{bsc}$ as a function of $d_p$ for representative coal ashes. To facilitate computer-based calculations, Laubscher and Rousseau [53] presents the following correlations derived from the graphs (with $d_p$ given in μm):

$$Q_{abs}=0.275d_p^{0.298}-0.305$$

(77)

$$Q_{bsc}=\frac{6.2188\times {10}^{-3}-1.0492\times {10}^{-2}{d}_p+7.287\times {10}^{-3}{d}_p^2-2.1925\times {10}^{-5}{d}_p^3}{1.851\times {10}^{-1}-2.0405\times {10}^{-3}{d}_p^2+6.254\times {10}^{-4}{d}_p^3}$$

(78)

The corresponding absorptivity ${\alpha }_g$ of the gas mixture at the wall temperature $T_w$ can be calculated using Equation (70) and Table 1 by replacing the gas temperature with the wall temperature. The absorption coefficient of the gas mixture can then be determined from:

$${\kappa }_{\alpha ,g}\left(T_w\right)=-\frac{\ln \left(1-{\alpha }_g\left(T_w\right)\right)}{S_{fur}}$$

(79)

The associated optical thickness ${\mathsf{\Phi }}_{\alpha ,gp}$ of the gas particle mixture is then given by

$${\mathsf{\Phi }}_{\alpha ,gp}\left(T_w\right)=\left({\kappa }_{\alpha ,g}\left(T_w\right)+Q_{abs}{A}_p{L}_p+k_{coke}{x}_1{x}_2{p}_{fur}\right)S_{fur}\gamma$$

(80)

Finally, the effective absorptivity of the flue gas and particle suspension at the wall temperature is then given by

$${\alpha }_{gp}\left(T_w\right)=\left(1-\beta \right)\left(\frac{1-e^{-{\mathsf{\Phi }}_{\alpha ,gp}\left(T_w\right)}}{1+\beta e^{-{\mathsf{\Phi }}_{\alpha ,gp}\left(T_w\right)}}\right)$$

(81)

2.6. Radiative-Convective Heat Exchanger Model

Superheaters and reheaters are typically tubular heat exchanger banks inserted in the flue gas stream downstream from the furnace to heat the steam to the required turbine inlet conditions. These heat exchangers are usually classified as either the radiant type, where thermal radiation is the predominant heat transfer mechanism, or the convective type, where convection is the predominant mechanism. Radiant-type superheaters and reheaters are located either in the furnace zone (platen zone) or in the convective pass zone closer to the furnace, whereas the convective superheaters and reheaters are found either in the convective pass zone further away from the furnace or in the back-pass zone of the boiler.

However, all the heat exchangers transfer heat via a combination of radiation and convection, where the ratio between the two phenomena varies according to the temperature of the flue gas mixture surrounding the tubes and the flue gas velocity due to the spacing between the tubes. For this reason, it is convenient to formulate a generic one-dimensional radiative–convective heat exchanger model where the ratio between the mechanisms is determined automatically based on the local operating conditions. This model effectively provides the ‘component characteristics’ of the heat exchanger that were referred to earlier. The principles encapsulated in this model can be applied to superheaters, reheaters, economisers and air heaters while correctly accounting for the differences in geometry.

Figure 2 provides a schematic showing all the relevant heat transfer phenomena together with the fluid flow paths in such a generic radiative/convective heat exchanger model.

The steam flow path (in blue) with inlet and outlet enthalpies $h_{st,i}$ and $h_{st,e}$, respectively, is separated by the heat exchanger wall from the flue gas flow path (in red) with inlet and outlet enthalpies $h_{fg,i}$ and $h_{fg,e}$, respectively. A fixed mass flow rate of fly ash (in black) with inlet and outlet enthalpies $h_{fa,i}$ and $h_{fa,e}$, respectively, is carried along with the flue gas, while a specified mass flow rate of ingress air (in green) with enthalpy $h_{ia}$ may also leak into the flue gas flow path. The effects of the ingress air on the overall mass balances and flue gas mixture compositions were addressed in Equations (9) and (11).

Heat is extracted from the flue gas in several ways. The first is via a combination of convective and radiative heat transfer to the heat exchanger surfaces from the flue gas directly surrounding it, as indicated with ${\dot{Q}}_{hx}$. The radiation heat transfer involved here is simply referred to as gas and particle radiation (as opposed to direct radiation). The second is via a combination of convective and radiative heat transfer to the furnace water walls that surround the cavity in which the heat exchanger is situated, as indicated with ${\dot{Q}}_{wall}$. The third is via a combination of convective and radiative heat transfer to the furnace roof that may form part of the enclosure surrounding the cavity in which the heat exchanger is situated, as indicated with ${\dot{Q}}_{roof}$. The fourth is the direct radiation from the hot flue gas inside the heat exchanger cavity radiating out the back towards the next heat exchanger surface, which is indicated with ${\dot{Q}}_{back}$.

In addition to this, there is also direct radiation impinging on the heat exchanger surface from the previous heat exchanger situated in the flue gas flow path or from the furnace in the case of the platen heat exchanger situated behind the furnace exit. This is indicated by ${\dot{Q}}_{ri}$. A portion of this incoming direct radiation bypasses the heat exchanger, as indicated by ${\dot{Q}}_{bp}$, while the rest is absorbed by the heat exchanger walls, as indicated by ${\dot{Q}}_{abs}$. ${\dot{Q}}_{bp}$ together with ${\dot{Q}}_{back}$ makes up the total direct radiation impinging onto the next heat exchanger surface, as indicated by ${\dot{Q}}_{re}$.

These heat transfer rates are all formulated to provide positive values. The total rate of heat transfer to the flue gas needed as a source term in Equation (10) will have a negative value and is given by the following energy balance

$${\dot{Q}}_{fg}=-\left({\dot{Q}}_{hx}+{\dot{Q}}_{wall}+{\dot{Q}}_{roof}+{\dot{Q}}_{back}\right)$$

(82)

The other energy balances internal to the heat exchanger model are

$${\dot{Q}}_{abs}={\dot{Q}}_{ri}-{\dot{Q}}_{bp}$$

(83)

$${\dot{Q}}_{re}={\dot{Q}}_{bp}+{\dot{Q}}_{back}$$

(84)

and

$${\dot{Q}}_{st}={\dot{Q}}_{hx}+{\dot{Q}}_{abs}$$

(85)

${\dot{Q}}_{st}$ provides the source term required for the energy balance of the steam given in Equation (7).

We will now consider the calculation of each of the heat transfer rates in more detail.

2.6.1. Convection and Gas Radiation Heat Transfer

For an increment of a heat exchanger surface, for which the fluid properties may be assumed to be constant, the rate of heat transfer can be characterised by the overall UA value. The overall UA value is equal to the inverse of the thermal resistance between the two fluid streams. Considering fouling and/or scaling on the outside and inside of the tubes, the UA (W/K) value can be written as

$$\frac{1}{UA}=\frac{1}{{\eta }_i{h}_i{A}_i}+\frac{R_{fi}}{{\eta }_i{A}_i}+R_w+\frac{R_{fo}}{{\eta }_o{A}_o}+\frac{1}{{\eta }_o{h}_o{A}_o}$$

(86)

$h_i$ (W/m²K) and $h_o$ (W/m²K) are the convective heat transfer coefficients on the inside and outside, respectively, $A_i$ (m²) and $A_o$ (m²) are the total heat transfer area on the inside and outside, respectively, ${\eta }_i$ and ${\eta }_o$ are the overall surface efficiency on the inside and outside, respectively, $R_{fi}$ (m²K/W) and $R_{fo}$ (m²K/W) are the fouling factors on the inside and outside, respectively, and $R_w$ is the thermal resistance of the tube wall material.

In the case of platen heat exchangers, the outside heat transfer area is taken as the projected flat surface areas formed by the tubes that are stacked closely together alongside one another in the flow direction. For all other heat exchangers, it is the total wetted outside surface area of all the tubes that are in contact with the flue gas.

To account for the gas radiation on the flue gas side, the outside convective heat transfer coefficient is replaced with an effective convective/radiative heat transfer coefficient

$$h_o=\xi \left(h_{o,c}+h_{o,r}\right)$$

(87)

where $h_{o,c}$ (W/m²K) is the convective heat transfer coefficient and $h_{o,r}$ (W/m²K) is the effective gas radiative heat transfer coefficient. A correction factor $\xi$ is also included to compensate for the non-uniform sweeping of the flow across the outer surface, with typical values of 1.0 for cross-current flow and 0.95 for mixed-current flow. In the case of a platen heat exchanger, a multiplication factor is included to account for the fact that the convection heat transfer area is the total wetted outside surface area in contact with the flue gas and not only the projected flat surface area. Therefore, for platen heat exchangers, we have

$$h_o=\xi \left(\left(\frac{\pi d}{2S_L}\right)h_{o,c}+h_{o,r}\right)$$

(88)

with $d$ (m) representing the tube outer diameter and $S_L$ (m) representing the longitudinal pitch (spacing in the flow direction) between the tube rows.

The value of $h_{o,r}$ is obtained by again considering the radiation heat transfer given in Equation (59). According to Zhang [44], the coefficient $\frac{{\epsilon }_w}{{\alpha }_{gp}+{\epsilon }_w-{\alpha }_{gp}{\epsilon }_w}$ can be approximated with $\frac{1+{\epsilon }_w}{2}$. Further assuming that the flue gas directly surrounding the heat exchanger surface is an isothermal and diffuse-grey medium leads to ${\epsilon }_{gp}={\alpha }_{gp}$. Therefore, the radiative heat transfer between the gas and the tube surfaces is given by

$${\dot{Q}}_{rad}=\frac{1+{\epsilon }_w}{2}{\epsilon }_{gp}\sigma A_{rad}\left(T_g^4-T_w^4\right)$$

(89)

The effective heat transfer coefficient due to radiation can then be obtained from ${\dot{Q}}_{rad}=\left(\frac{1+{\epsilon }_w}{2}{\epsilon }_{gp}\sigma \frac{T_g^4-T_w^4}{T_g-T_w}\right)A\left(T_g-T_w\right)=h_{o,r}A\left(T_g-T_w\right)$, which implies that

$$h_{o,r}=\frac{1+{\epsilon }_w}{2}{\epsilon }_{gp}\sigma \frac{T_g^4-T_w^4}{T_g-T_w}$$

(90)

The procedures described in Section 2.3 are also used here to determine the flue gas and particle suspension emissivity ${\epsilon }_{gp}$ in the heat exchanger cavity. For this purpose, the mean beam length for platen heat exchangers is given by

$$S_{hx}=\frac{1.8}{\frac{1}{H}+\frac{1}{b_L}+\frac{1}{S_T}}$$

(91)

$H$ (m) is the average height of the heat exchanger, $b_L$ (m) is the depth in the flow direction determined from $b_L=\left(N_L-1\right)S_L$ with $N_L$ representing the number of tube rows in the flow direction, and $S_T$ (m) is the transvers pitch (spacing perpendicular to the flow direction) between the tube rows. For all other heat exchanger tube banks, it is given by

$$S_{hx}=0.9d\left(\frac{4S_T{S}_L}{\pi d^2}-1\right)$$

(92)

with $d$ (m) representing the outer tube diameter.

Figure 3 shows a single increment of a convective/radiative heat exchanger with the various thermal resistances between the flue gas and steam flow paths.

Note that $h_{o,r}$ in Equation (90) is a function of the outer wall temperature of the heat exchanger, which according to Figure 3 will be the outer surface temperature of the fouling layer on the outside $T_{fo}$. Given this, it is convenient to assume a uniform surface temperature and to divide the total thermal resistance into two thermal resistances in series given by

$${\frac{1}{UA}|}_o=\frac{1}{{\eta }_o{h}_o{A}_o}$$

(93)

and

$${\frac{1}{UA}|}_i=\frac{1}{{\eta }_i{h}_i{A}_i}+\frac{R_{fi}}{{\eta }_i{A}_i}+R_w+\frac{R_{fo}}{{\eta }_o{A}_o}$$

(94)

The heat transfer from the flue gas to the outside of the outer fouling layer can then be written in terms of the well-known effectiveness-NTU method as

$${\dot{Q}}_{hx}={\epsilon }_{hx,fg}{\dot{m}}_{fg}\left(T_{fg,i}-T_{fo}\right)$$

(95)

where ${\epsilon }_{hx,fg}$ is the effectiveness on the flue gas side defined as

$${\epsilon }_{hx,fg}=1-e^{-\left(\frac{{UA|}_o}{{\dot{m}}_{fg}{c}_{p,fg}}\right)}$$

(96)

with $c_{p,fg}$ (J/kgK) is the specific heat capacity of the flue gas.

Similarly, the heat transfer from the outside surface of the outer fouling layer to the steam can be written in terms of the effectiveness-NTU method as

$${\dot{Q}}_{st}={\epsilon }_{hx,st}{\dot{m}}_{st}\left(T_{fo}-T_{st}\right)$$

(97)

where ${\epsilon }_{hx,st}$ is the effectiveness on the steam side defined as

$${\epsilon }_{hx,st}=1-e^{-\left(\frac{{UA|}_i}{{\dot{m}}_{st}{c}_{p,st}}\right)}$$

(98)

with $c_{p,st}$ representing the specific heat capacity of the steam.

Knowing the value of ${\dot{Q}}_{abs}$ from Equation (83), the simultaneous solution of Equations (85), (95) and (97) yields the values of ${\dot{Q}}_{hx}$, ${\dot{Q}}_{st}$ and $T_{fo}$.

Note that the LMTD method may be employed here instead of the effectiveness-NTU method, since these are exactly equivalent. However, the authors prefer the latter, since it is less problematic for computer-based implementation.

There are several correlations available for calculating the relevant convective heat transfer coefficients based on the geometry of the specific heat exchanger surface. For instance, for flow over tube bundles, the reader is referred to [54], and for flow within tubes, the reader is referred to [55].

The heat transfer to the water walls (${\dot{Q}}_{wall}$) and roof (${\dot{Q}}_{roof}$) surrounding the heat exchanger is usually simply calculated as

$${\dot{Q}}_{wall}=h_o{A}_{wall}\left({\overline{T}}_{fg}-T_{wall}\right)$$

(99)

and

$${\dot{Q}}_{roof}=h_o{A}_{roof}\left({\overline{T}}_{fg}-T_{roof}\right)$$

(100)

with $A_{wall}$ and $A_{roof}$ representing the wetted wall and roof heat transfer areas, respectively, $T_{wall}$ and $T_{roof}$ representing the wall and roof temperature, respectively, ${\overline{T}}_{fg}$ representing the average flue gas temperature and $h_o$ being obtained from Equation (88) with $0.5<\xi <1.0$.

2.6.2. Direct Radiation towards the Next Heat Exchanger Surface

The direct radiation from the hot flue gas inside the heat exchanger cavity radiating out the back towards the next heat exchanger surface, namely ${\dot{Q}}_{back}$, is required in the energy balance of Equation (84) and can be determined from [22,46]

$${\dot{Q}}_{back}={\zeta }_r{\epsilon }_{gp}\sigma A_e{\overline{T}}_{fg}^4$$

(101)

$A_e$ is the outlet area behind the heat exchanger, which is also the inlet area to the next heat exchanger, ${\overline{T}}_{fg}$ is the average flue gas temperature within the heat exchanger cavity, and ${\zeta }_r$ is a correction factor for the type of fuel, which is equal to 0.5 for coal and liquid fuels.

2.6.3. Direct Radiation Bypassing the Heat Exchanger

The direct radiation that bypasses the heat exchanger surfaces, namely ${\dot{Q}}_{bp}$, is required in the energy balances of Equations (83) and (84), and it is given by [22,48]

$${\dot{Q}}_{bp}=\frac{{\dot{Q}}_{ri}\left(1-{\epsilon }_{gp}\right)x_e}{\beta }$$

(102)

$\beta$ is the same as in Equation (54) for the platen heat exchanger at the outlet of the furnace, while $\beta =1$ for all other heat exchangers. $x_e$ is the angular coefficient and is given by

$$x_e={\left[{\left(\frac{b_L}{S_T}\right)}^2+1\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-\frac{b_L}{S_T}$$

(103)

If the heat exchanger is situated behind an empty cavity in the flue gas flow path, there will be additional direct radiation from the hot gas inside the cavity to the surface of the heat exchanger. A simplified approach to include this is to apply a correction to the effective gas radiative heat transfer coefficient as follows

$$h_{o,r}^{’}=h_{o,r}\left(1+A{\left(\frac{T_{fg,cav}}{1000}\right)}^{0.25}{\left(\frac{b_C}{b_L}\right)}^{0.07}\right)$$

(104)

with $T_{fg,cav}$ representing the flue gas temperature inside the cavity, $b_C$ (m) representing the depth of the cavity in the flow direction and $b_L$ (m) representing the depth of the heat exchanger tube bundle as before. $A$ is an empirical constant which is equal to 0.3 for fuel oil and gas, 0.4 for bituminous coal and anthracite, and 0.5 for lignite.

An alternative to this simplified approach is to apply a cavity model as presented below.

2.7. Cavity Model

The convective heat transfer between the flue gas and the walls surrounding a cavity is typically assumed to be negligible. However, the radiative heat transfer to the wall may be significant, depending on the flue gas temperature, and there will be direct radiation entering and leaving the cavity. Therefore, the gas radiation to the wall surrounding the cavity wall is given by ${\dot{Q}}_{hx}$ in Equation (95) but with $h_{o,c}=0$. Furthermore, since the primary heat transfer surface consists of the surrounding wall and roof, we have ${\dot{Q}}_{wall}={\dot{Q}}_{roof}=0$. All the other components of heat transfer are calculated in the same way as for the radiative-convective heat exchanger model, including ${\dot{Q}}_{bp}$, ${\dot{Q}}_{abs}$, ${\dot{Q}}_{back}$, ${\dot{Q}}_{re}$, ${\dot{Q}}_{fg}$ and ${\dot{Q}}_{st}$.

2.8. Evaporator Model

The EV component consists of all the water walls in the furnace and the SD combined, plus all water walls surrounding the flue gas duct, and the roof tubes at the top of the boiler if these are connected in series with the water walls. Therefore, the total heat transfer to the steam in the evaporator is given by

$${\dot{Q}}_{st,EV}={\dot{Q}}_{rad,ww}+\sum {\dot{Q}}_{wall}+\sum {\dot{Q}}_{roof}+{\displaystyle \sum _{cav}{\dot{Q}}_{hx}}$$

(105)

Remember that ${\dot{Q}}_{st,EV}$ is needed in the steam energy balance given by Equation (7) and in the steam drum level control given by Equation (8).

2.9. Membrane Water Wall Model

In Section 2.4, we have seen that the projected method for furnace heat transfer does not account for the temperature of the water wall surface. The reason for this is that even if it was considered, the calculated radiation heat transfer to the water walls is relatively insensitive to the wall temperature. However, in some cases, it may be important to determine the temperature distribution within the water wall to ensure that the maximum metal temperature is within the allowed range. The maximum metal temperature occurs in the fin section and could be much higher than the average wall temperature, depending on the specific geometry. Rousseau and Laubscher [27] proposed a thermofluid network-based model for membrane walls of boiler furnaces that allows the prediction of the temperature distribution within the water wall tubes and fins. A schematic of the proposed network model consisting of equivalent 1D thermal resistances is shown in Figure 4.

The model takes as input the uniform heat flux onto the planar projected inner surface facing the furnace, which would be equal to ${\dot{Q}}_{rad,ww}/A_{rad,ww}$ in the furnace model, the bulk steam/water temperature within the water wall tubes, and the inside convective heat transfer coefficient. It then provides as outputs the area-weighted average temperature of the total planar projected surface area exposed to the furnace, which would be $T_w$ in the furnace model, as well as the overall $UA$ value of the water wall, which in the furnace model would be the inverse of the effective wall thermal resistance $R_w$. Given this, the membrane wall model can be connected to the furnace model and solved simultaneously to provide the wall temperature and heat flux. Finally, the membrane model then provides the maximum temperature at the centre of the fin facing the furnace. The reader is referred to the original paper in [27] for further details.

3. Application Case Studies

This section presents the results of two case studies prepared by the authors with the results compared to site measurements to demonstrate the application of the modelling methodology described in Section 2.

3.1. Utility-Scale Coal-Fired Boiler

The first case study that will demonstrate the application of the modelling methodology is a utility-scale 620 MW_e pulverised coal-fired subcritical drum-type boiler located in Southern Africa.

3.1.1. Boiler Layout and Operating Conditions

The unit has a maximum continuous rating (MCR) of 473 kg/s of steam at an exit pressure and temperature of 165 bar(a) and 540 °C, respectively. In addition to the high-pressure steam circuit, the boiler has a medium pressure reheat circuit which delivers steam at 37 bar(a) and 540 °C and a flow rate of 457 kg/s. The boiler has a two-pass configuration, with fully welded membrane walls surrounding the combustion chamber, superheaters, reheaters and economiser. The combustion chamber is fed by 36 wall-mounted swirl burners, which direct the heated primary air (PA) and secondary air (SA) along with the fuel into the furnace. The general arrangement of the boiler is shown in Figure 5 (left).

Combustion gases generated in the furnace travel up to the platen superheater (SH2), after which they travel to the final superheater (SH3), then the final reheater (RH2), primary superheater (SH1), primary reheater (RH1), and finally the economiser (EC). Once the gas exits the EC, it is used to preheat the incoming combustion air in the air heater (not shown). The combustion chamber is approximately 64 m high and 24 m wide. Both SH2 and SH3 are widely pitched (SH2-1143 mm, SH3-828 mm) platen-type superheaters used to cool the flue gas down before it reaches the other more tightly pitched heat exchangers (RH2-229 mm, SH1-304 mm, RH1-153 mm and EC-77.8 mm).

The as-received ultimate analysis of the coal burnt in the boiler is provided in

Table 2. Coal composition (from actual fuel analysis).

Constituent	C	H	O	N	S	Ash	Moisture
Units	wt %	wt %	wt %	wt %	wt %	wt %	wt %
Value	0.4156	0.0222	0.079	0.0097	0.0094	0.409	0.055

. From the composition, it is noted that the fuel has a high ash concentration (>35 wt %), which together with the applied excess air ratios (

Table 3. Coal-fired boiler case study model inputs (from site measurements).

	Load Case Values
Parameter	100% MCR	80% MCR	Units
Fuel
HHV of fuel	15,070	15,070	kJ/kg
Unburned carbon	0.021	0.021	kg/kg
Radiation loss percentage	0.0056	0.0056	W/W
Air
Ambient air temperature	306	306	K
Ambient air pressure	≈87,100	≈87,100	Pa
Relative humidity	≈50	≈50	%
AH outlet air temperature	561	559	K
Excess air ratio	1.15	1.2	-
Water
Inlet water temperature	520.8	510	K
Inlet water pressure	17.2	16.7	MPa
Steam drum pressure	16.9	16.5	kPa(a)
ATT 1 flow rate	38.9	29.5	kg/s
ATT 2 flow rate	3.1	2.3	kg/s
Final steam flow rate	472.1	386.5	kg/s
Reheat water
Inlet water temperature	607	601	K
Inlet water pressure	3.7	3	MPa
ATT RH flow rate	13.5	13.5	kg/s
ATT water temperature	429.9	430	K
Final steam flow rate	457.9	374.5	kg/s

) leads to a flue gas particle loading of $0.013k{g}_{ash}/m_{fg}^3$. This is well above the limiting threshold of $0.005k{g}_{ash}/m_{fg}^3$ sited for the low particle load model. The HHV of the fuel is provided in Table 3 along with the assumed furnace radiation loss and unburned carbon percentages.

For the case study boiler under consideration, numerous CFD studies have been performed previously by the authors to investigate phenomena such as the effect of swirl direction on furnace heat uptake, heat uptake of radiant superheaters for varying loads and effect of high ash on furnace heat transfer [2,13,38]. The CFD model was developed using ANSYS Fluent^® 2019 and includes the evaporator waterwalls and the radiant superheaters directly downstream of the furnace outlet plane. The model is based on the steady-state Eulerian–Lagrangian approach to approximate the flue gas and fuel particle flow and heat transfer. The conservation equations are solved using the respective Reynolds averaged forms, with particle tracking accounting for the momentum and energy interaction between the phases. A separate transport equation is solved for each chemical constituent in the reaction mechanism. The turbulent quantities in the transport equations are closed using the realisable k epsilon turbulence model, and the radiation transport is modelled using the Discrete Ordinates Method (DO). The diffusion-kinetic model was applied to model the char combustion reaction rates and the eddy dissipation model for the gas phase chemical reaction rates. For a more detailed description of the model, the reader is referred to [2].

During these studies, the CFD model has been thoroughly validated using various site measurements at different load conditions. Figure 5 (right) shows a plot of the flue gas temperature results for the 100% MCR case along with the furnace exit temperature of 1678 K.

For the present case study, two load conditions were simulated using the process modelling methodology described in Section 2, namely 100% and 80% MCR. The measurements collected earlier by the authors are again employed in the current case study. However, there are no flue gas temperature measurements available. Therefore, selected gas temperatures calculated with the process model will be compared to the CFD results from the above-mentioned studies. The various input boundary conditions taken from site measurements are provided in Table 3. The excess air ratios were determined using O₂ and CO₂ measurements at the exit of the EC.

3.1.2. Thermofluid Network Model

The PFD of the thermofluid network for the process model is shown in Figure 6. It consists of eight heat exchanger elements (numbered in square blocks) and various nodes (numbered in circles) per fluid stream between these heat exchangers (air, high-pressure water, reheater water and flue gas). Note that the high-pressure water circuit has two attemperators, namely at the inlet of SH2 and at the inlet of SH3, which are both fed with incoming feed water piped from the inlet of the EC. The reheat circuit has a single attemperator at the inlet of RH2, which is fed with water taken from the pre-boiler plant.

In this boiler, the roof tubes convey saturated steam from the top of the steam drum to SH1 while picking up more heat from the flue gas. The total EV heat load consists of the heat transfer rate to the furnace water walls together with the walls and roof surrounding heat exchangers 2–7. Therefore, the steam at the outlet of the EV is slightly superheated rather than saturated. In the diagram, the steam drum is simply shown for illustrative purposes, and node 3 represents the outlet of the evaporator (EV) circuit, which is also the inlet to SH1.

The network model consists of 460 variables and 460 equations that must be solved simultaneously. These equations include all the mass, energy and momentum (momentum for the water side only) balance equations similar to those described in Section 2.1 as well as all the component characteristic equations as described in the rest of Section 2. The equations were custom programmed and solved using the Engineering Equation Solver (EES) software package (https://fchartsoftware.com/ees/, accessed on 13 December 2022). It is a general equation-solving program that can numerically solve large sets of coupled non-linear algebraic and differential equations by applying a modified Newton–Raphson method. It also includes a built-in high accuracy thermodynamic and transport property database. The equations are solved down to a maximum relative residual of 10⁻³, which is defined as the difference between the left-hand and right-hand sides of the equation divided by the magnitude of the left-hand side of the equation. The time required for the solution is approximately 30 s on a laptop computer with an Intel(R) i7 7500CPU processor and 16 GB of RAM.

The basic projected method described in Section 2.4.1 was applied in the furnace heat transfer calculations, and the generic radiative-convective heat exchanger model described in Section 2.6 was applied for each of the heat exchangers. To demonstrate the role of the flue gas and particulate emissivity model, the analysis was conducted using both the high ash loading model presented in Section 2.5.3 (M1) and the standard model presented in Section 2.5.1 (M2).

In both models, the final steam flow rate was fixed at the mean measured value, and the attemperator spray water flow rates were set equal to the measured values. To ensure that the correct heat is transferred to the water flowing through the EV circuit, the fuel flow rate was then adjusted in M1 until the predicted temperature at node 3 equalled the mean measured value. In this manner, the appropriate firing rate was obtained while enforcing the other boundary conditions and inputs shown in Table 2 and Table 3. For direct comparison, the same fuel flow rate was then also applied in M2.

The appropriate excess air ratio was also found by adjusting it in the process model until the calculated dry volume fraction of O₂ in the flue gas equalled the value measured on site.

3.1.3. Results and Discussion

The site measurements were taken over a period of approximately five hours, and the mean, minimum and maximum values were extracted from the data, as provided in

Table 4. Comparison between measured and calculated water/steam values for utility-scale coal-fired boiler model for 100% MCR load case.

			Measurements			Model
Parameter	Node	Units	Mean	Min	Max	M1	M2
High-pressure water
Feed water temperature *	1	K	520.8	520.8	520.8	520.8	520.8
Feed water pressure *	1	kPa(a)	17,200	17,000	17,400	17,200	17,200
EC outlet temperature	2	K	566.3	558.7	572.2	560.7	561.4
Steam drum pressure *	3	kPa(a)	16,900	16,600	17,100	16,900	16,900
SH1 inlet temperature	3	K	625.1	624.3	626.6	626.1	638
SH1 outlet temperature	4	K	682.6	679	688.2	678	700
SH2 outlet temperature	6	K	757.2	746.4	766.4	740.6	753
SH3 outlet temperature	8	K	808.9	806.9	810.7	803	807.6
SH3 outlet pressure	8	kPa(a)	16,200	16,000	16,500	16,100	16,010
Final steam flow rate *	8	kg/s	472.1	469.6	475.8	472.1	472.1
Reheat water
RH1 inlet temperature *	1	K	607.2	606.2	608.4	607	607
RH1 inlet pressure *	1	kPa(a)	3700	3700	3700	3700	3700
RH1 outlet temperature	2	K	709	705.1	714.3	687	689
RH2 inlet temperature	3	K	673.8	662.8	684.8	654	657
RH2 outlet temperature	5	K	808	804.7	811.5	808	798
Final steam flow rate *	5	kg/s	457.9	455.5	461.5	457.9	457.9

* Values are inputs to model.

. Using the M1 model, a fuel flow rate of 99.5 kg/s was obtained for the 100 MCR case, with a dry volume fraction of O₂ in the flue gas of 2.8 V/V%, while the measured value is 2.7 V/V%. Table 4 also shows the various values predicted with the M1 and M2 process models for the 100% MCR case.

Table 4 shows that the M1 and M2 models predict similar EC exit water temperatures and that both model predictions fall within the min/max range of the measurements. For the inlet water temperature of SH1, the M1 prediction compares well with the mean measured value, but M2 significantly overpredicts the SH1 inlet steam temperature. This implies that the rate of heat transfer to the EV circuit is overestimated, resulting is an excessive degree of superheating at the inlet of SH1. The overprediction of the inlet temperature of SH1 by M2 also leads to an overprediction of the outlet water temperature of SH1 by approximately 2.6%, whereas M1 slightly underpredicts, but by less than 1%. Both M1 and M2 provide good results for the outlet temperatures of SH2 and SH3. For the reheat circuit, both models underpredict the RH1 outlet water temperature by approximately 3%, but M1 performs better when considering the outlet water temperature of RH2.

Figure 7 (left) shows the heat transfer rates per heat exchanger calculated from the measured values along with the values predicted using M1 and M2. Figure 7 (right) shows the Q-T diagram, which has the flue gas temperatures on the vertical axis and the cumulative heat transfer rates on the horizontal axis.

The results show that M1 and M2 overpredict the EV heat load, i.e., the furnace heat transfer rate, by 3.5% and 13.7%, respectively. The furnace exit gas temperature calculated using the CFD model (Figure 5 right) is 1678 K, whereas the furnace exit gas temperatures calculated using M1 and M2 are 1572 K and 1496 K, respectively. This again shows that both models overpredict the furnace heat transfer rate, thereby underpredicting the furnace exit gas temperature, and much more so in the case of M2. These results show the importance of utilising the high particle load model for the effective emissivity and absorptivity of the flue gas and particle suspension, rather than the standard model, especially in the case of a high ash fuel. It also confirms the statement by Brummel [48] that at higher values of mean beam length (such as in the furnace) and at very high particle loads, “a thermal insulation of the hot core zone of a gas–solid mixture can occur. The radiation of the core zone cannot reach the enclosure walls anymore, as the layer of the gas solids mixture close to the walls turns out to be opaque, insulating gas as well as particle radiation from the core zone”.

In Section 2.5.1, it was pointed out that Equation (68) was taken from Zhang [44] and that a different version is provided by Lin [45] and Basu [46]. If we replace the version by Zhang [44] with the version of Lin [45] and Basu [46] in the M2 model, it will result in the EV heat transfer being about 4% lower for this case study. Although this brings it closer to the measured value, it would still overpredict the EV heat load by 8.7%. The uncertainty about Equation (68) as well as the overprediction of the EV heat loads lead the authors to prefer the use of the high ash loading model presented in Section 2.5.3 (M1) rather than the standard model presented in Section 2.5.1 (M2).

For the other heat exchangers, we see that M1 and M2 perform similar with the relative error percentages ranging from 3.3% (SH2) to 26.7% (RH1) for M1 and from 6.1% (RH2) to 22.4% (RH1) for M2.

Figure 8 shows the contribution of the different heat transfer mechanisms for each heat exchanger predicted using the M1 and M2 models. It shows that hardly any of the direct radiation from the furnace bypasses SH3. Therefore, the platen-type superheaters (SH2 and SH3) are achieving the goal of shielding and cooling the gas temperature below the ash deformation temperature (1250 K) before it enters the closely pitched convective type heat exchangers (RH2, SH1, RH1 and EC). The higher furnace exit gas temperature in M1 results in higher gas and particulate radiation in SH2 and SH3 compared to M2. Figure 7 showed that M1 overpredicts the steam heat transfer rate in RH2 water. This is due to the higher gas and particulate heat transfer rate shown in Figure 8, which in turn is due to the higher calculated inlet gas temperatures for RH2.

A study was also performed at 80% MCR to demonstrate the applicability to part load conditions using the M1 model only. The model inputs are also provided in Table 3. The fuel flow rate was again adjusted until the predicted SH1 inlet water temperature matched the mean measured value. Using this strategy, a fuel flow rate of approximately 82 kg/s was obtained. Figure 9 shows that the M1 model results compare well with the site data with a maximum relative error of 24% (RH1) and minimum relative error of 1% (EV). The predicted furnace exit gas temperature is 1510 K, which is within 5.5% of the value of 1599 K from the CFD prediction [2].

In addition to the above analyses, the direct method for furnace analysis presented in Section 2.4.3 was compared to the results obtained with the projected method for the 100% MCR load case, both employing the high particle load model. The results are provided in

Table 5. Direct method and projected method heat transfer rate results for 100% MCR coal-fired boiler case study.

Heat Exchanger	EV	SH1	SH2	SH3	RH1	RH2	EC
Direct method	581	153	185	91	83	159	83
Projected method	559	161	190	95	85	163	84
Experimental values	540	175	196	78	108	141	97

together with the measured values. The results show that the direct method overpredicts the EV heat load by 7.6%, compared to the projected method, which overpredicts by 3.5%. The projected method also fares better around the platen superheater situated behind the furnace exit (SH2), with an underprediction within 3% of the measured value as opposed to the direct method with an underprediction within 5.6% of the measured value.

3.2. Industrial Biomass-Fired Boiler

The second case study is an industrial-scale sugarcane bagasse fired boiler located in Southern Africa. The sugarcane burnt in the combustion chamber is pre-processed in a cane diffuser; therefore, the raw fuel has a small fibrous texture and burns in semi-suspension. In other words, a portion of the fuel burns on the stoker grate and the remainder burns in the freeboard above the grate.

3.2.1. Boiler Layout and Operating Conditions

The selected boiler has a maximum continuous rating of 28.8 kg/s of steam flow at a temperature and pressure of 3 MPa(a) and 400 °C, respectively. The steam is used to generate electricity and process heat for a sugar mill. The general arrangement of the boiler is shown in Figure 10 (left).

The boiler consists of a moving-grate stoker, membrane-type furnace walls, two bare tube superheaters (SH1 and SH2), spray-type attemperator, condenser, an evaporator tube bank (EV) which forms part of the steam generation circuit along with the furnace, single-pass tubular air heater (AH), and an extended surface economiser (EC). The combustion chamber is approximately 22 m high and 6.5 m wide. The two inline superheaters are both transversely pitched at 200 mm, and the downstream EV tube bank is pitched at 100 mm. SH2 consists of four longitudinal rows, SH1 of 10 rows and the EV bank of 23 rows. For the air heater, flue gas travels inside the tubes and ambient combustion air passes across the tube bank. The AH has a staggered arrangement with transverse and longitudinal pitches of 90 mm and 70 mm, respectively, and with tube diameters of 63.5 mm. The economiser consists of tubes with outside diameters of 50.8 mm, which are pitched transversely and longitudinally at 125 mm. The EC fins are 2.8 mm thick square fins that are 120 mm wide.

The as-received ultimate analysis of the biomass fuel can be seen in

Table 6. Biomass fuel composition (from actual fuel analysis).

Constituent	C	H	O	N	S	Ash	Moisture
Units	wt %	wt %	wt %	wt %	wt %	wt %	wt %
Value	0.2127	0.027	0.2094	0.0016	0.0002	0.047	0.496

. Important differences to note between the biomass fuel and the coal in the previous case study are that for the biomass fuel, the moisture content is significantly higher (49.6% vs. 5.5%), and there is a significantly lower ash mass fraction (4.7% vs. 40.9%). Using the excess air ratios listed in

Table 7. Biomass-fired boiler case study model inputs (from site measurements).

	Load Case Values
Parameter	100% MCR	62% MCR	Units
Fuel
HHV of fuel	8894	8894	kJ/kg
Unburned carbon	0.007	0.007	Kg/kg
Radiation loss percentage	0.00514	0.00514	W/W
Air
Ambient air temperature	298	305	K
Ambient air pressure	101,325	101,325	Pa
Relative humidity	95	85	%
Excess air ratio	1.24	1.24	-
AH outlet air temperature	493	493	K
Water
Inlet water temperature	371	370	K
Spray water temperature	487	509	K
Final steam pressure	2773	3095	kPa(a)
Desired final steam temperature	675	675	K
Final steam flow rate	28.8	18.8	kg/s

results in a particle loading of $0.002k{g}_{ash}/m_{fg}^3$, which is well below the threshold value of $0.005k{g}_{ash}/m_{fg}^3$ for low particle loading. The HHV of the biomass fuel along with the assumed radiation loss from the furnace and the assumed unburned carbon are shown in Table 7.

As with the coal-fired boiler model in Section 3.1, a thoroughly validated CFD model of the biomass boiler was created as part of a previous study by the authors [3]. The modelling approach used is the same as for the coal-fired boiler except for the fuel firing system. Whereas the pulverised coal suspension-fired boiler has wall-mounted swirl burners, the biomass boiler has pneumatic fuel spreaders and a moving grate with the primary air forced upwards through the grate. The grate is therefore modelled as a porous zone and special scripts are used to account for the grate-particle interaction. This includes the translation motion of the grate and the particle drag and gravitational forces. For a more detailed description of the model, the reader is referred to [3]. Figure 10 (right) shows the side view of the temperature contour results for the 100% MCR load case along with the predicted furnace exit temperature of 1264 K. For the present case study, two load conditions were simulated using the process modelling methodology described in Section 2, namely 100% and 62% MCR. Site measurements collected earlier by the authors are again employed in the current case study. However, there are only flue gas temperature measurements available at the exit of the EV bank, AH and EC. Therefore, the CFD results will be used to evaluate the furnace exit gas temperature predictions obtained with the process models.

The input boundary conditions obtained from the site measurements are also provided in Table 7. The excess air ratios were again determined from O₂ and CO₂ measurements taken at the exit of the EC.

3.2.2. Thermofluid Network Model

The PFD of the thermofluid network for the process model is shown in Figure 11. The network consists of seven heat exchanger components and various fluid nodes between the heat exchangers. Feed water enters the economiser (EC) and is heated by the flue gas. The feed water then flows through a shell-and-tube condenser and is further pre-heated while condensing steam that is extracted from the steam drum. The pre-heated feedwater (node 3) is then fed into the steam generation circuit which comprises the furnace water walls and the membrane walls surrounding SH1 and SH2 in parallel with the EV tube bank. The generated steam is then collected in the steam drum from where it is piped to the inlet of SH1. From there, the steam is heated to its final conditions (node 8) in SH2. The liquid water generated in the condenser (node 10) is used as attemperator spray water at the inlet of SH2.

The air heater (AH) is explicitly modelled in this case. The AH has a bypass control mechanism to ensure that the temperature of the air entering the combustion chamber does not exceed 493 K. This is included in the process model between nodes 1, 2 and 3 of the air circuit. Furthermore, in the boiler, only 85% of the combustion air passes through the AH, which makes up the primary air. The primary air is fed from below upwards through the stoker grate bars. The remaining 15% of the air is fed into the furnace at atmospheric conditions as overfire air and spreader air. The latter is used to physically fling the fibrous fuel onto the boiler grate.

The network model consists of 392 variables and 392 simultaneous equations that were again solved using the EES software package.

As in the first case study, the basic projected method was applied in the furnace heat transfer calculations and the generic radiative-convective heat exchanger model for each of the heat exchangers. The analysis was made using both the high particle loading model (M1) and the standard model (M2).

In both models, the final steam flow rate was fixed at the mean measured value. The fuel mass flow rate was then adjusted in M1 to ensure that only saturated steam leaves the evaporator circuit at the top of the steam drum (node 5). Since the attemperator flow rate was not measured, it was adjusted in M1 to ensure that the final steam temperature (node 8) is equal to the desired final steam temperature shown in Table 7.

3.2.3. Results and Discussion

The site measurements were taken over a period of approximately two hours, and the mean, minimum and maximum values were extracted from the data, as provided in

Table 8. Comparison between measured and calculated water, flue gas and air values for industrial-scale biomass-fired boiler model for 100% load case.

			Measurements			Model
Parameter	Node	Units	Mean	Min	Max	M1	M2
Water side
SH1 outlet temperature	6	K	661.6	651.3	669	640	639.8
ATT outlet temperature	7	K	566.1	555	576.6	567	568.6
Spray water temperature *	10	K	487	472	506	487	487
SH2 outlet pressure *	8	kPa(a)	2773	2228.7	3054.4	2773	2773
SH2 outlet temperature	8	K	674.4	668.4	687	675.2	676.6
Final steam flow rate *	8	kg/s	28.8	25.3	31	28.8	28.8
Steam drum pressure	5	kPa(a)	3115.2	2660.7	3357.9	3053	3053
EC inlet water temperature *	1	K	370.2	262.2	376.2	371.2	371.2
EC outlet water temperature	2	K	460.2	456.2	465.2	462	458.6
Flue gas side
EV outlet temperature	5	K	701.2	695.2	712.2	705.2	693.9
AH outlet temperature	6	K	570.9	565.2	576.2	603.2	595.2
EC outlet temperature	7	K	445.2	441.2	448.2	424	422
Air side
AH inlet temperature *	1	K	298.5	298.2	299.2	298.5	298.2
AH outlet temperature	3	K	492.2	488.2	497.2	492.4	486

* Values are inputs to model.

. Using the M1 model, a fuel flow rate of 12.92 kg/s was obtained for the 100% MCR case. The measured O₂ fraction was 4.2% (minimum 3.1% and maximum 6.2%), and the measured CO₂ was 16% (minimum 12.5% and maximum 17.5%), while the values obtained with the M1 model after adjustment of the excess air ratio are 4.2% for O₂ and 16.6% for CO₂. Table 8 also shows the other values predicted with the M1 and M2 process models for the 100% MCR case.

The results show that both M1 and M2 provide good results for the steam temperatures that are within 0.4% at the EC outlet and within 3.3% at the SH1 outlet. M1 only slightly overpredicts the flue gas temperature at the EV exit by 0.6%, while M2 underpredicts by 1.0%. The flue gas temperature at the outlet of the AH is overpredicted by both M1 (5.7%) and M2 (4.3%), while the outlet temperature of the EC is underpredicted by both M1 (4.8%) and M2 (5.3%). The air outlet temperature of the AH obtained with M1 is within 0.04% of the measured value, while for M2, it is within 1.3%. The differences in the results of M1 and M2 are much smaller than in the case of the coal-fired boiler.

Figure 12 (left) shows the heat transfer rates per heat exchanger calculated from the measured values along with the values predicted using M1 and M2, while Figure 12 (right) shows the Q-T diagram of the flue gas temperatures and cumulative heat transfer rates for both model results, excluding the condenser. The results show that M1 and M2 produce very similar results that are mostly within the min/max range of the measured values. For M1, the deviations from the measured values range between 0.9% (EC) and 4.2% (SH1), and for M2, they range between 1% (EC) and 10% (SH1).

M1 overpredicts the EV heat load, i.e., the furnace heat transfer rate, by only 1.7%, while M2 overpredicts by 4.0%. The furnace exit gas temperatures calculated using M1 and M2 are 1206 K and 1151 K, respectively, compared to the CFD result of 1264 K. M1 therefore underpredicts by 4.6%, while M2 underpredicts by 8.9%, which is consistent with the overpredictions of the evaporator load. However, the difference between M1 and M2 is much less in this case (2.2%) than in the case of the high-ash coal-fired boiler (9.8%). These results show that for the biomass case with its low particle loading, comparable results are obtained with the high-particle load model and the standard model. The high-particle load model is therefore also applicable in cases of low particle loads. This confirms the statement by Brummel [48] that “It is generally recommended to apply the more sophisticated coupling calculation model (referring to the high particle load model) even for operating conditions with low particle loads…” It also confirms the authors’ preference for using the high ash loading model presented in Section 2.5.3 (M1) rather than the standard model presented in Section 2.5.1 (M2), even in cases with low ash loading.

The breakdown of heat transfer phenomena for each heat exchanger is shown in Figure 13. For SH2, which is directly exposed to the furnace, the direct radiation is significant, and the gas and particle radiation is of similar magnitude as the convection heat transfer. From the EV onwards, convection is by far the dominant mechanism. This is different from what was found in the coal-fired boiler, where gas and particle radiation continued to play a notable role, even far downstream from the furnace. The reasons for this are that in the biomass case, the furnace exit gas temperature is much lower (and therefore, the flue gas temperatures are lower in general), there are much lower entrained particle concentrations, and the heat exchangers tubes are more tightly pitched, thereby enhancing convection.

To showcase the capability of the process model to capture low load performance of the boiler, an additional simulation was performed and compared to site data. This load case was for the boiler operating at 62% MCR. Again, the fuel flow rate along with the attemperator flow rate were adjusted to find the correct steam generation rate and final steam temperature. The fuel flow rate and attemperator flow rate are found to be 8.2 kg/s and 1.75 kg/s respectively. Figure 14 shows the predicted heat transfer rates per heat exchanger using model M1, along with the measured values. The results show that the results from the process model compare well with the site data at part load, with a minimum error of 0.5% (EC) and maximum error of 22% (SH2). The predicted furnace exit gas temperature is 1097 K, which is within 0.7% of the value of 1089 K from the CFD prediction [3].

The direct method for furnace analysis presented in Section 2.4.3 was also compared to the results obtained with the projected method for the 100% MCR case, both employing the high particle load model. The results are provided in

Table 9. Direct method and projected method heat transfer rate results for 100% MCR biomass boiler case study.

Heat Exchanger	Furnace+ EV	SH1	SH2	AH	EC
Direct method	51	11.2	8.3	6.9	11.8
Projected method	53.5	9.9	7.8	6.6	11.2
Experimental values	52.6	10.3	7.6	6.7	11.1

together with the measured values. For the biomass case, the direct method underpredicts the furnace and EV bank heat load by 3%, while the projected method slightly overpredicts the furnace and EV bank heat load by 1.7%. The projected method fares better downstream with SH1 underpredicted by 3.9% as opposed to an overprediction of 8.7% for the direct method and SH2 overpredicted by 2.6% as opposed to an overprediction of 9.2% for the direct method. The overpredictions by the direct method is due to the underprediction of the furnace heat transfer rate, which leads to higher furnace exit gas temperatures and consequently higher predicted heat absorption rates by SH1 and SH2.

4. Summary and Conclusions

The integrated boiler process modelling methodology described here is built on a generic framework that is applicable to each component within the boiler. It consists of the simultaneous solution of the fundamental balance equations of mass, energy and momentum written between the inlet and the outlet, the component characteristic equations that provide values for the source terms, the appropriate fluid property relationships, and the boundary values that may be specified at the inlet and/or outlet. The key component models are the adiabatic flame temperature model, which includes the combustion analysis, the furnace model, and the generic convective-radiative heat exchanger model.

The adiabatic flame temperature model assumes instantaneous and complete combustion, except in the case of solid fuels such as coal and biomass, where a fraction of unburned carbon is assumed. It also employs real fluid properties in the gas mixture analysis referred to the standard reference state rather than the ideal gas assumption. This requires the HHV to be used in the energy balance equation rather than the LHV.

There are two zero-dimensional models typically employed in the furnace analysis, namely the projected method and the direct method. These provide the radiation heat transfer to the water wall and the direct radiation leaving the furnace exit plane and impinging onto the heat exchangers situated behind the furnace exit. For the coal-fired boiler 100% MCR case study, the projected method overpredicts the evaporator heat load by 3.5%, while the direct method overpredicts by 7.6%. For the 100% MCR biomass case, the projected method slightly overpredicts the furnace and evaporator bank heat load by 1.7%, while the direct method underpredicts the furnace and evaporator bank heat load by 3%. The projected method therefore provides better results overall.

An important aspect throughout the analysis is the calculation of the effective emissivity and absorptivity of the flue gas and particle suspension. For this, there are three possible approaches, namely the ‘standard model’, the ‘low particle load model’, and the ‘high particle load model’. The results of the case studies show the importance of utilising the high particle load model rather than the standard model in the case of a high ash fuel such as the current coal-fired boiler case study. However, it can also be used with confidence in the case of low particle loading, such as the current biomass-fired boiler case study. The uncertainty about the equation for calculating the fly ash pressure-based absorption coefficient in the standard model, as well as the resultant overprediction of the EV heat loads, leads the authors to prefer the use of the high ash loading model rather than the standard model even in cases with low ash loading.

A generic radiative–convective heat exchanger model was presented. It accounts for the combination of convective and radiative heat transfer to the heat exchanger surfaces from the flue gas directly surrounding it, the direct radiation from the hot flue gas inside the heat exchanger cavity radiating out the back towards the next heat exchanger surface, the direct radiation impinging on the heat exchanger surface from the previous heat exchanger situated in the flue gas flow path, the portion of this incoming direct radiation that bypasses the heat exchanger, as well as the rest that is absorbed by the heat exchanger walls. The interplay between the different heat transfer mechanisms is determined automatically based on the local operating conditions. It can therefore be applied to superheaters, reheaters, economisers and air heaters, provided that the differences in geometry are correctly considered.

The case studies showed that the results obtained with the process models using the projected method and the high particle load model compare well with the values obtained from site measurements and the detail CFD models for full load and at part load, even though the two boilers have very different geometries, firing systems and fuels.