## 1. Introduction

Currently we observe an increase in natural gas utilization, deployment of renewable energy sources, and a need to improve coal plant efficiency, which tend to decrease coal consumption in OECD. In spite of the above, coal will long remain a key energy fuel for electricity generation in a number of developed countries. Performance optimization of large-scale pulverized coal (PC) utility boilers has become more and more relevant through the recent years for the utility industry. Optimization efforts are focused on increasing thermal efficiency, extending their lifetime, and lowering pollutant emissions.

The coal-fired power plants face a big challenge of improving combustion process enforced by environmental concerns. The new units designed for high steam parameters must meet strict limits for emissions of gaseous pollutants. Meeting the requirements of CO_{2} emissions while maintaining highly efficient production process is still the subject of many research [1–6]. Regulations for reducing NO_{x} and SO_{x} emissions become more strict and meeting them must be achieved by both new and old production units. It is currently possible to achieve the required quality of exhaust gasses with primary and secondary measures. Primary measures of NOx reduction are used inside the boiler combustion chamber and consist of activities such as proper selection of the excess air ratio and temperature, as well as modification of combustion techniques (reburning, exhaust gas recirculation, air staging, cooling the flame, and burners redesign). Secondary measures of NOx reduction include activities and auxiliary installations located behind the combustion zone of boiler like selective non-catalytic reduction (SNCR) and selective catalytic reduction (SRC) methods. The efficiency of the NOx reduction of primary measure is about 35% (reduction from 540 to 350 mg/m^{3}). Secondary measures are more effective. SNCR achieve the NOx reduction efficiency of up to 50% and efficiency of SCR method can reach about 95%.

However, in the case of older units, applying these methods in the wrong way can greatly affect the quality of the combustion process and consequently reduce production efficiency. To keep the thermal efficiency of the boiler constantly at a high level, it is necessary to perform optimization work. In recent years, the increasing number of activities using optimization methods based on computer simulations is observed. The optimization activities are carried out to improve performance of coal and biomass cofiring process [7], and decrease the carbon content in ash [8]. Many investigations were devoted to the modeling and prediction of NOx emission [9–11] as well as to optimization methods for NOx reduction [12–14]. Current research on optimization methods are based on data from available measurement systems. To further optimize pulverized coal boiler performance, dedicated systems for monitoring and control of boiler operating conditions [15, 16] can be required. This optimization process can be done by optimizing the amount and distribution of fuel and air supplied to the boiler [17]. Flame shape and its symmetry is a critical factor influencing combustion process performance. Optimization of combustion process can improve thermal efficiency of boilers up to 0.84% [18]. Homogenous temperature distribution of flue-gas promotes lower emissions of NO_{x}, CO, and minimizes total organic carbon (TOC) content in ash. Acoustic gas temperature and mass flow rate of pulverized coal measurement systems allow for online monitoring of the most important combustion parameters. A temperature measurement can also be applied to verify the results of computational fluid dynamics (CFD) modeling [19].

The next big challenge for the electricity producers will be the flexible operation with regard to increasing share of renewable sources of energy in the domestic markets. The increase of electricity produced from renewable energy sources means that coal units have to adapt to market needs and work with increased load variations during the day. Working under changing load requires constant monitoring and controlling the operating parameters of the boiler in order to maintain high efficiency and required levels of gas emissions. Conducted research also shows that the time required to load changes can be much shorter than the one that is currently used [20]. As a result, in the near future, coal-fired units will work more often in transient conditions which will require the use of new methods presented in this paper, to control and optimize the combustion process also at transient conditions.

In recent years the development of computer modeling techniques, in particular the computational fluid dynamics, has allowed for accurate analysis of coal combustion process in tangentially fired [21–26] as well as front-fired furnaces [27–29]. Results of CFD modeling [30] can be used at the stage of the combustion process optimization to provide complete information about the process. A comprehensive large-scale furnace CFD model should be capable of properly predicting trends in NO_{x} reduction by means of primary [31] and secondary [32] methods.

## 2. Coal combustion process monitoring and control

Air flow as well as fuel flow distribution in a power plant has to be measured and controlled in order to achieve optimum combustion conditions and improve the boiler efficiency. In low-NOx combustion system, the air and fuel flow has to be precisely controlled. Both the total air distribution (primary and secondary air to the boiler) and the mass flow rate through the individual burner. **Figure 1** presents dependence of the key combustion parameters (NO_{x}, O_{2}, CO, LOI, and boiler efficiency) on the air/fuel ratio. To optimize the combustion process the accurate measurements and control systems are needed.

The combustion stoichiometry can be easily controlled by air distribution. For low air-fuel ratio value the O_{2} versus CO dependency is high. As long as the amount of air is not optimized in reference to fuel the CO can reach very high values. While increasing O_{2} the CO declines—the combustion is getting complete. The optimum zone is defined as low CO value for the lowest possible O_{2} supply, while the excess of O_{2} has to be maintained for combustion stabilization. For optimum zone the NOx value also has to be as low as possible. O_{2} has crucial influence on NOx generation. The NOx level is proportional to the amount of O_{2}, depending on NOx generation stoichiometry. If the excess of air will reach the highest values the NOx will also grow. For those conditions we start to lose the combustion efficiency. The output losses will grow rapidly according to O_{2} value. This can be observed by highest FEGT.

### 2.1. Low NOx emission technique of pulverized coal combustion

The purpose of improving NO_{x} emissions resulted in the use of low-emission combustion technology, which is characterized by staged combustion process with extended low oxygen combustion zone. In pulverized coal boiler, the low-NO_{x} combustion technology is most often implemented using low-emission burners and spatial distribution of air supplied to the furnace. Furnace installations can also be equipped with a system of fuel staged supply to the boiler. This process involves providing coal-air mixture with better fineness in the higher regions of furnace. General scheme of low-emission combustion system in PC boiler OP-650 (EDF Poland, Rybnik Power Plant) is presented in **Figure 2**.

OP-650 is a front wall fired, pulverized coal boiler with steam drum and natural circulation. Maximum continuous rating is 650 t/h of live steam generation (700 t/h in peak load). Boiler is equipped with six ring-ball mill units (A, B, C, D, E, and F) that supply 24 burners. Mills are equipped with static classifiers with movable blades.

Low-emission combustion system consists of two rows of main swirl burners located at the lower part of windbox (six burners in each row). Additional diluted coal-air mixture is provided through 12 drop tubes located at the upper part of windbox (**Figure 3**). Boiler also consists of two levels of overfired air nozzles (OFA) located at the front wall (6 OFA ports) and rear boiler wall (10 OFA ports). All burners are fitted in the common windbox located on the front wall of the boiler. The burners’ arrangement was designed to reach the low-NO_{x} combustion conditions with air and fuel staging method. The first two rows of burners (low-NO_{x} burners) (burners I and II) are supplied with a high concentration mixture and the third row (burners III—without swirl) is supplied with low concentration mixture. Depending on load demand, the mill units operation configuration is different (some of the mills are out of operation).

The main swirl burner design with marked flow zone of coal-air mixture and secondary air is presented in **Figure 4**. Combustion air is separated into core, primary, secondary, and tertiary as is shown in **Figure 4**. Additional diluted coal-air mixture is provided through 12 burners located in a single row at the upper part of windbox. These burners are made of drop tubes with additional six nozzles for optional biomass injection.

### 2.2. Monitoring and control of fuel flow rate

Pulverized coal separation into low and high concentration mixture is realized by means of louver splitters installed behind mill outlet. Each mill has two outlets (multiturret), where each outlet is equipped in splitting box. Behind each splitter there are two separate coal pipes which transport the low and high concentration mixture to the burners (**Figure 5**). The fines of high concentration mixture is approximately 30% and 5% residue for R88 and R200 mesh sieve, respectively, and low concentration mixture fineness is approximately 15% (R88) and 1% (R200).

Fuel distribution is controlled in open loop by unit operator. Each louver has individual actuator which is adjustable from distributed control system (DCS) in specific range from 0% to 100% which corresponds to 30–70% share of mixture density. **Figure 6** depicts the main idea of fuel splitter operation.

The use of louver separator allows to control the pulverized fuel (PF) distribution in share range 60–80% for high concentration mixture and 20–40% for low concentration mixture. According to this PF distribution, the low emission combustion condition can be met by means of fuel staging method.

The set point for louver position depends on current live and reheated steam temperatures and current PF distribution which is related to mill unit operation configuration. The fuel splitters settings allow to control the fire ball (flame core) location along the boiler height. This feature is often used to achieve the nominal live and reheated steam temperature (540°C). At low unit load, where the temperature is too low, the fire ball is moved up in the combustion chamber—the high concentration mixture is driven into burners placed on third level (burners III). In case of very high temperature, the fire ball is moved down—the high concentration mixture is driven for first and second burners level (burners I and II). The more complicated situation is where the mill unit operation is changed. The PF distribution differs, which cause the imbalance between left and right side and generates nonuniformity of gas temperature. For PF distribution and velocity measurement the dedicated online measuring system is used. Each coal pipe is equipped with sensors which are connected to the measuring portable cabinets (each of four coal pipes, one for mill unit). Depending on location constraints, different section lengths (~1.2–2.4 m) of pipes with internal diameter (ID) of 0.46–0.51 m were chosen for the system installation. The general rule is to install the sensors in vertical pipes at a location of three internal diameters (3 × ID) for the inlet and one internal diameter (1 × ID) for the outlet from bends and curves. If there is no possibility to use the system for vertical pipes the horizontal sections are chosen.

Additionally, online optimization of the PF distribution is performed by unit operation, while the periodic optimization is carried out by process engineer personnel. Due to the fact that the system is portable and each mill has separate measuring cabinet, the boiler/mill group optimization is possible. The optimization strategy is simple. Initially, PF mass flow and velocities are analyzed for each mill separately to meet the following requirements: 20–40% mass flow for burners III (low concentration mixture), and 60–80% mass flow for burners I and II (high concentration mixture) for low-NOx combustion system. Next, PF flow is analyzed for boiler according to the results for each mill and consequently the symmetry between left and right sides of the furnace optimization is reached. Once satisfactory pulverized coal distribution was achieved, the system is moved to another group of mills/boiler.

### 2.3. Monitoring and control of air flow rate

The air distribution and flue gas systems in presented front-fired PC boiler consist mainly of three forced draft (FD) fans and three inducted draft (ID) fans (**Figure 7**). Two FD fans supply primary air to the mills and one FD fan supplies the secondary air to the boiler through the common windbox. The boiler is also equipped with two protection air fans, and one bottom air fan which supply the air to protect the boiler’s evaporator (protective air) and to burn out the fuel near bottom ash hopper (bottom air).

The air distribution control process is supported by dedicated measurement systems of air temperature, mass flow rate, and pressure. These measured values are required for evaluation of set points in control system: fuel demand, combustion air demand, and steam temperatures. In addition, the boiler is also equipped with system for combustion quality monitoring, where the flue gas composition measurement subsystem (O_{2}, NO_{X}, CO, and SO_{2}) is the most important component.

The most commonly used air flow measurement system is based on differential pressure *dP* value. Differential pressure measurement signal is used in primary air flow measuring providing the possibility to avoid the typical problems connected with air purity and thermal condition. The main disadvantage is the pressure drop in air duct which can limit the fan’s maximum capacity. For different air duct profile other similar measurements are used, i.e., venturi, orifice, annubar, veribar, or pitot tubes. New advanced technologies for air flow measurement are based on electrostatic, thermal, or optical principles and are free of difficulties related to complex geometry. They are also more accurate than *dP* measurements; however, pose some limitations depending on measuring technology.

The boiler operation is controlled by distributed control system, in which control algorithms for fuel/air amount are implemented. The presented power unit OP-650 operates mainly in coordinated boiler follow mode. The load demand signal is applied to the turbine valves (turbine master) and combustion demand (boiler master). Boiler master is the main control signal for boiler load algorithms which acts on fuel demand control algorithm for coal feeders’ capacity control. Combustion air (primary and secondary) demand is controlled according to oxygen rate in flue gases at the boiler outlet. Basing on oxygen measurement, the load of FD fans is evaluated and controlled. The air distribution system, secondary air and OFA, can be also controlled by unit operator in manual mode. In the primary and secondary air distribution optimization of advanced air balance control, some requirements have to be met. The primary/secondary air balance should be maintained in the range of 20–30%/70–80% of the total air flow.

The amount of total air flow is measured by means of *dP* measurement on main ducts, which are common for all the three FD fans. The unit controls the load of FD1 and FD3. The FD2 is out of control in this loop. The total amount of air is controlled indirectly by amount of secondary air. For automatic total air amount control the live steam mass flow measurements have been used together with O_{2} content in flue gases correction. The automatic control unit acts on two FD fans (FD1 and FD3) and on subordinate control unit responsible for primary air to the mills control (**Figure 8**), which respectively acts on FD2 load. The leading signal of live steam mass flow is added with correction signal from O_{2} content in flue gas behind the boiler. The controlled O_{2} content is kept on *O*_{2Z} set point level:

where *O*_{2z}, *O*_{2z}(*M*_{P}), and Δ*O*_{2z}(*O*) are oxygen set point at the boiler outlet, oxygen value as a function of steam mass flow rate *M*_{p}, and difference of oxygen content at the left and right side of boiler outlet, respectively.

The difference between those two signals is given on the input of *R*_{O2} (O_{2} controller), which evaluate the correction signal from O_{2} content. This correction is added to subordinated controller. The signal of total air flow to the boiler follows the leading signal of live steam flow recalculated in VM module (**Figure 8**). The correction from O_{2} and power changes is as follows:

where *V*_{C} is the total air flow rate, *V*_{P} and *V*_{L} are air flow rates on the right and left side of boiler, respectively. The difference of those two signals is introduced on two controllers of secondary air (R_{W1} and R_{W2}), which controls the load of FD1 and FD3 fan.

The main controller is proportional-integral type and works accordingly to tracking control idea. The set point signal is the live steam mass flow value, the controlled variable is the total air flow to the boiler. The correction from O_{2} content in flue gas is comprised in subordinated PID controller with O_{2} set point value. The deviation controls directly the adjuster of steam/air share. From practical point of view, the range ±20% for steam/air share is enough to keep the set value of O_{2} content at the boiler outlet.

The amount of primary air is measured by use of dedicated *dP* orifices between the mill fan and mill inlet. The primary air is dependent on the fuel mass flow rate supplied to the burner. The control process of air flow can only be achieved by modifying distribution of the secondary air supplied to both the burners and the OFA nozzles. The amount of secondary air is not measured, but is balanced by use dedicated performance calculations implemented in DCS. The secondary air is controlled by the position of the dampers presented in **Figure 9(a)**. Each burner has the ability to individually control the position of dampers, as opposed to the OFA nozzles. In the case of OFA nozzles on the front wall, the control process is carried out by changing the position of two adjacent nozzle dampers. On the rear wall of the control process is achieved by changing the damper position of the OFA nozzle located on the left and right side of the boiler (**Figure 9b**).

### 2.4. Protective air system

An adverse effect of low-emission combustion technology is the creation of a permanent reducing atmosphere (the reducing zone with reduced O_{2} and increased share of CO) in the boundary layer of screens. If close to the screens of combustion chamber O_{2} content is below 1%, and the CO content rises above 0.2%, the process of corrosion of screen tubes can occur (low-oxygen corrosion). In order to protect screen against negative effects of low-oxygen corrosion in form of intense loss of the tube material, boiler is equipped with protective air nozzles. The air supplying through nozzles located on the rear and the side walls of boiler (**Figure 10**) allows to increase the share of oxygen in the boundary layer.

The share of protective air in total amount of air supplied to the boiler can be up to 20%, and it needs be taken into account in any process related to the optimization of air distribution to the boiler.

### 2.5. Acoustic measurement system of flue-gas temperature

The acoustic gas temperature measurement system (AGAM) allows to asses impact of operating conditions on the temperature distribution inside the combustion chamber. Measurement system is located at the level of 30.2 m and is composed of 8 heads (**Figure 11**), which form the 21 paths of measurements. Based on the 21 paths, the flue-gas temperature can be calculated in 12 zones to create maps of temperature distribution at the outlet of combustion chamber.

The principle of AGAM system is based on the acoustic pyrometry, which allows to measure delay time of an acoustic wave. Time delay depends on the temperature in the environment of wave propagation and this relation is defined by the following formula:

where *v* is the speed of sound (m/s), *R* is the universal gas constant (kJ/(kmol∙K)), *M* is the molar mass of gas (kg/kmol), and *k* is the adiabatic index.

The adiabatic index and molar mass of gas are calculated based on the typical average values of volume fraction for flue-gas in coal-fired furnaces. For analyzed boiler, the parameters are equal to: *k* = 1.275, *N*_{2} = 76.5 Vol.%, O_{2} = 3 Vol.%, CO_{2} = 13 Vol.%, H_{2}O = 7.5 Vol.%.

Velocity of a wave is determined by the propagation time of an acoustic impulse between symmetrically placed transmitters and receivers. The distances between the transmitter and the receiver are constant and the temperature is calculated using the formula presented below:

where *l* is the distance between the transmitter and the receiver (m), *τ* is the time of delay (s), and *B* is the acoustic constant.

## 3. CFD modeling techniques in pulverized coal-fired boilers

Performance optimization of large-scale pulverized coal utility boilers has become more and more relevant through the recent years for the utility industry. CFD modeling has been extensively applied to provide information on the complex phenomena in tangentially fired [21–27] and in front wall-fired furnaces [27–29], including gas-solid flow, combustion, and heat transfer.

In [21], two tangentially fired furnaces (OP-380 and OP-430) of similar design and thermal power were numerical investigated. OP-380 was retrofitted by replacing traditional jet burners with RI–JET2 (Rapid Ignition JET-burner) swirl burners. Comparison between the boilers combustion characteristics was carried out based on the CFD simulation. In [22], different operation regimes of pulverized coal furnace have been investigated with the 3-D CFD code. Selected parameters have been compared to measurements, showing good agreement. Similar purpose was achieved in [23]. Additionally, a grid refinement study was performed. Pulverized coal ignition behavior in a 40 MW tangentially fired boiler was predicted in [24]. Ignition image was obtained from high-temperature-resistant camera and compared to simulation results. Accuracy of general simulation approach was confirmed by available operating and design data. Yin et al. [25] investigated furnace and part of the rear pass in the tangentially fired boiler. The simulation has been validated with global design parameters including O_{2} at the furnace outlet, heat transfer in the furnace, and furnace exit temperature. Site operation data was used to verify NO_{x} predictions. In [26], Euler-Lagrangian approach was incorporated to investigate numerical flow characteristics in a tangentially fired furnace. Temperature deviation has been analyzed. An example of employing a commercial code Fluent to investigate a 500 MW utility boiler firing medium volatile coal was demonstrated in [27]. Temperature profiles were calculated for different boiler loads. Calculations have been compared with measurements data. Minghou et al. [28] addresses CFD simulations of front-fired boiler for different operating conditions. The model was validated by comparing measured unburned carbon in fly ash, NOx, and total heat transfer to the walls with measured values. The combustion and wall heat flux in a 100 MW boiler under air and oxy-fuel combustion conditions was analyzed by means of computational modeling in [29]. No validation of the numerical approach was done against the experimental data.

Reducing a complex physical problem to a series of models that can be solved numerically requires a number of assumptions to be made. Specifically, for engineering problems momentum and species transport equations must be modeled. Simulations are computed using CFD code, which solves Reynolds averaged Navier-Stokes equations using a low-order finite volume formulation. In the current work, the steady-state solution is calculated using second-order discretization for all equations.

Simulation of the following processes takes place in the furnace: turbulent flow, coal combustion, gas phase combustion, particle transport, and radiative transport. The gas phase is modeled assuming an Eulerian approach, while for the solid phase both the Lagrangian as well as Euler-Euler approach is applied.

Realizable *k*-*ε* model [34] was used as a closure of turbulent Reynolds equations. The realizable *k*-*ε* model is relatively widely used for engineering applications and provides better performance in many industrial turbulent flows than the standard *k*-*ε* model. The flow near the wall is influenced by molecular viscosity rather than by turbulence. The wall function method of [35] uses algebraic formulations to link quantities at the wall to those further away. The *y*+ values have been kept above 20.

### 3.1. Coal devolatilization

The pulverized coal combustion process can be divided into two parts, devolatilization and char combustion. The most commonly used conventional devolatilization model is based on a single kinetic rate [36], which assumes that the rate of devolatilization is dependent on the amount of volatiles remaining in the particle via a first-order reaction:

where *k*, *f*_{v,0}, and *m*_{p,0} denote reaction rate (following Arrhenius rate), volatile fraction, and initial particle mass. These models can be extended by using network devolatilization models like FG-DVC and FLASHCHAIN [37, 38] as a preprocessor. An example of using FG-DVC model with assumed particle heating rate equal to 10^{5} K/s is demonstrated below. It predicts the rate of the production and high temperature yields for the char, tar, volatiles, and the composition of key species during the devolatilization of any coal. The results as well as the proximate and ultimate analysis for the used coal are given in **Table 1**.

During primary devolatilization each volatile species is evolved with different rate and tar undergoes secondary pyrolysis [39]. In the CFD modeling of turbulent flow with combustion it was assumed that volatiles are produced as a single compound that undergoes instantaneous breakup reaction into tar, light hydrocarbons, CO, CO_{2}, and H_{2}O. FG-DVC calculates devolatilization rate of tars and mentioned species. The most significant mass drop of fuel particle occurs when tar is evolved, which is produced in first place (**Figure 12**). For this reason tar release rate is used in devolatilization model.

Knowing the volatile fraction of dry ash free coal (*f*_{volatile}) and assuming that residual char is pure carbon, we can calculate lower heating value of volatiles.

Assuming that lower heating value of light hydrocarbons is approximately equal to that of methane (*LHV*_{gas} = 50 MJ/kg), we can easily calculate lower heating value of tars from the instantaneous breakup reaction of volatiles:

where *y*_{gas}, *y*_{tar}, and *y*_{CO} stand for mass fraction in volatiles.

A novel approach in devolatilization simulations was described as tabulated devolatilization process model (TDP) [40]. The authors indicated that in the previously mentioned modeling approach all the devolatilization parameters used in the CFD code were calculated assuming a constant heating rate for all the particles (with a typically used value of 10^{5} K/s). In the TDP approach a look-up table of the devolatilization rates and yields for a set of heating rates is generated based on the FLASHCHAIN model. In the course of CFD simulations, a particular set of devolatization parameters is used based on the calculated particle heating history.

In nonpremixed reacting flows the local time-dependent mixing and chemical reaction of the species and heat transfer away from the reaction area determine the combustion process. The key gas phase combustion modeling issue is the necessity of source terms calculations in reactive species transport equations, which are the average values of strongly nonlinear reaction rates. Early combustion models have been derived on the assumption of chemical equilibrium. Taking into account detailed kinetics of reactions usually results in much higher computational effort.

#### 3.1.1. Eddy dissipation concept (EDC)

Eddy dissipation concept [41] was used as a general concept for treating interaction between turbulence and chemistry in flames. In this model the total space is subdivided into fine structures and the surrounding fluid. All reactions of the reactive components are assumed to react only in these spaces which are locally treated as perfectly stirred reactors (PSR) with a residence time:

where *ν* is the kinematic viscosity and *ε* denotes turbulent kinetic energy dissipation rate. These parameters are calculated from turbulence model (*k*-*ε*). Mass fraction occupied by fine structures is modeled as:

The reaction rates of all species are calculated on a mass balance for the fine structure reactor. Denoting quantities with asterisk, the conservation equation of species *i* can be defined:

where *i*, *M*_{i} is the molecular weight of the species *i*, and *ω*_{i}*** denotes the chemical reaction rate calculated from Arrhenius equation.

The mean net mass transfer rate of species *i* between the fine structures and the surrounding fluid can be expressed as:

The implementation of the EDC model into CFD code is realized by solving the nonlinear system of equations for the fine structure reactor in each control volume and finding *R*_{i}, which is the source term in species *i* transport equation.

In most of the engineering cases implementation of the detailed reaction mechanism into 3D codes is prohibitive due to their large computational effort. For many purposes the required information can often be obtained with a less complete chemistry description. CFD simulations most often use simplified global reaction mechanism. Four-step global mechanism mainly based on the one demonstrated in [42] was employed. It contains four global reactions: hydrocarbon and tar decomposition reactions, carbon monoxide, and hydrogen oxidation:

Following the above mechanism with known heat of reactions we can further calculate enthalpies of formation of volatiles, light hydrocarbons, and tars. Enthalpy of formation of tars was calculated by assuming zero heat of volatiles instantaneous breakup reaction. Heat of pyrolysis was not included in the analysis.

#### 3.1.2. Mixture fraction/probability density function (PDF) method

Pulverized coal combustion process is a one of example of turbulent nonpremixed combustion systems and can be modeled using the mixture fraction/PDF model. Mixture fraction *f* can be expressed as the local fuel mass fraction [43], where *Y*_{F} and *Y*_{P} is mass fraction of fuel and products, respectively:

The main approach of diffusion models is that the combustion is limited only by mixing of fuel and oxidizer. The release of reaction products takes place if the fuel and oxidizer are locally mixed and does not depend on the reaction rate. In mixture fraction/PDF model transport equations for individual species are not solved. The mass fraction of oxidizer, fuel, and products is calculated based on the mixture fraction value *f*. **Figure 13** presents the graphical interpretation of mixture fraction approach.

If 0 ≤ *f* ≤ *f*_{S} fuel is deficient and the mixture is called fuel lean. The mass fraction of oxidizer and products is represented in the following form:

On the other side, if *f*_{S} < *f* ≤ 1 the mixture is called fuel rich, and the following equation can be used to calculate mass fraction of products and fuel:

where *f* is the mixture fraction, *f*_{S} is the stoichiometric mixture fraction, *Y*_{S} is the mass fraction of products of stoichiometric reaction at *f* = *f*_{S}, *Y*_{i} is the mass fraction linear functions (*Y*_{O}, *Y*_{F}, *Y*_{P}), *Y*_{O} is the local mass fraction of oxidizer, *Y*_{F} is the local mass fraction of fuel, *Y*_{P} is the local mass fraction of products, *Y*_{Ox} is the mass fraction of oxidizer at *f* = 0, and *Y*_{Fuel} is the mass fraction of fuel at *f* = 1.

The reaction is complete if the whole mass of fuel and oxidizer vanish (*Y*_{F} = 0, *Y*_{Ox} = 0) and this state is described by stoichiometric mixture fraction value *f*_{S}. The mixture fraction approach allows to calculate mass fractions at each control volume based on the one value and the modeling of combustion process is simplified to a mixing problem. The governing transport equation of mixture fraction is given by:

where *D* is a molecular diffusion coefficient (m^{2}/s), *u*_{i} is the velocity in direction of *i* (m/s), *x*_{i} is the direction *i* in Cartesian coordinates system (m), and *ρ* is the density (kg/m^{3}).

Equation (14) can be applied taking into account the assumption of equal diffusivities of fuel, oxidizer, and products. In turbulent flows the turbulent convection dominates in comparison with molecular diffusion and assumption of molecular diffusion coefficient equality is acceptable.

When local value of mixture fraction is known, it becomes possible to calculate the local enthalpy *h* of combustion products. The enthalpy of the burnt mixture is a linear piecewise function of *f* and is presented in **Figure 14**. Therefore, the calculation of the gas temperature and an amount of heat generated (heat losses) in the combustion process becomes possible.

For fuel-lean mixture and 0 ≤ *f* ≤ *f*_{S}:

and for fuel-rich mixture and *f*_{S} < *f* ≤ 1:

where *h*_{L} is the enthalpy on the lean side of *f*_{S} (kJ/kg), *h*_{R} is the enthalpy on the rich side of *f*_{S} (kJ/kg), *h*_{Ox} is the enthalpy of oxidizer at *f* = 0 (kJ/kg), *h*_{Fuel} is the enthalpy of fuel at *f* = 1 (kJ/kg), *h*_{S} is the enthalpy of products at the stoichiometric mixture fraction *f*_{S} (kJ/kg).

To determine the thermodynamical properties in turbulent flow, the probability density function of the mixture fraction is needed. The mean enthalpy can be computed from the following equation:

The conservation equation for the mean mixture fraction *β* function [43–45] used in coal combustion modeling has a presumed shape and depends only on the mean mixture fraction and its variance .:

where the normalization factor *B*(*a*, *b*) and is defined as:

and the PDF parameters *a* and *b* can be determined using and :

The *β* functions also have limitations and they cannot describe distributions joining an extreme peak (*f* = 0 or *f* = 1) with a maximum intermediate peak in the range of 0 < *f* < 1. The different approaches can be proposed as the Dirac peak at the boundary, in order to eliminate this inconvenience.

Despite simplifications the mixture fraction/PDF method allows the determination of the basic parameters of the combustion process, also in the case of coal combustion. The mixture fraction/PDF approach is commonly used in computational fluid dynamics and particularly in modeling of the turbulent reactive flows which are the most popular cases occurring in industrial practice.

### 3.2. Char combustion

Char undergoes heterogeneous oxidation to CO. The reaction rate is calculated on the assumption that the process is limited by the diffusion of oxygen to the external surface of the char particle and char reactivity [46]. It is assumed that particle is spherical in shape and that surface reaction rate depends on the ratio of the reacting surface to the surface of the sphere. Diffusion phenomena with the rate:

and kinetic reaction:

occur simultaneously. The total char combustion rate is than expressed as:

where *D*_{p} denotes particle radius, *D*_{0} is the diffusion coefficient, *ρ* is the particle density, *X*_{O2} is the oxygen mole fraction, *M*_{O2} is the oxygen mole fraction, *T*_{p} and *T*_{g} are the particle and bulk temperature, *A*_{f} is the empirical constant dependent on the fuel, and *φ* is the ratio of the reacting surface to the surface of the sphere. The particle burns out with constant diameter and variable density. It was assumed that the particle absorbs all the heat of the char burnout according to [47].

### 3.3. Radiative heat transfer

In case of combustion in furnace problem radiation is not only the dominant energy transport mechanism but also one of the most complex problems. The radiative transfer equation (RTE) [48] governs the radiation heat transfer in participating media. It describes the variation of radiation intensity (*I*) as it travels along a certain path (*s*) in the medium, in the direction (*s, ω*). Considering absorption, emission, and scattering apart from and into the direction *s*, the RTE can be described as follows:

where *k* and *σ* denote absorption and scattering coefficients, *Φ* is the phase function. The spatial integration of RTE was carried out with the discrete ordinates (DO) method [49]. The number of RTE depends on the total number of gray gases and takes into account scattering of the particles. The DO method solves the RTE for a set of directions based on the concept of angular discretization scheme. Each octant of the angular space 4π is discretized into polar and azimuthal solid angle. The continuous integral over the solid angle is approximated by a numerical quadrature scheme, where the equations are solved for a series of directions.

In a typical combustion chamber H_{2}O and CO_{2} are the main gaseous absorbers and emitters of radiant energy. The total emissivity of gas is calculated by a number of gray gases using polynomial correlations for weighing factors and absorption coefficient according to the weighted sum of gray gases method (WSGGM) [50]. Widely employed coefficients for emissivity [51], fitted from the benchmark exponential wide-band model, have been used in this work. The WSGGM represents the entire spectrum with three gray gases having uniform absorption coefficients. The total gas phase absorption coefficient is calculated from the total emissivity with the mean path length calculated from the characteristic cell size.

The gas phase absorption coefficient was corrected according to Taylor-Foster model [52], assuming uniform and constant soot concentration (10^{-3} kg/m^{3}) in the furnace. As noticed in [46], the main source of radiative transfer in two-phase mixture is the particle cloud. Therefore, coal particles emissivity (*ε*_{p}) treatment is crucial in coal combustion simulation. In this work the particle emissivity was assumed to be a function of unburned carbon in particle (*U*_{c}) following the relation [53]:

The particle reflectivity and scattering effects are also included in the calculation of heat transfer.

Thermal boundary conditions at walls have been expressed in terms of surface temperature and emissivity. It was assumed that the evaporator surface temperature is about 60° higher than the saturation temperature corresponding to the pressure of 16 MPa in the boiler drum. Calculations have been carried out for three emissivity values equal to 0.5, 0.7, and 0.9 at fixed wall temperature. All the presented figures correspond to emissivity equal 0.7, which is a typical value found in the literature [54]. We have to emphasize that during real boiler operation temperature and emissivity vary with time and spatially due to water-wall slagging and soot blowers operation. This phenomenon can be included in the simulations by implementing deposition model with thermal properties submodel [55].

## 4. Coal combustion process optimization

The main goal in power generation optimization is related to reduction of operational and maintenance costs. Form practical point of view it is directly dependent on two factors: first, combustion effectiveness (effectiveness in fuel consumption), and second, production effectiveness by optimal operation in steady and transient state (process performance). The combustion effectiveness depends on boiler performance and operational skills which are strictly related to the control system. To increase the combustion process performance the fuel consumption and flue gas emission (NO_{X}, CO, SO_{2}) have to be reduced.

### 4.1. Optimization method

The combustion process optimization has been performed for wall-fired pulverized coal boiler of EDF Polska S.A. in Rybnik Power Plant. The main goal was to minimize the boiler losses (maximize boiler efficiency) and to keep the environmental requirements for NO_{X} and exploitation requirements like live and reheat steam temperatures on optimal level.

The general optimization method used in the EDF Rybnik Power Plant was based on online calculation of manipulated variable (MV) according to controlled variable (CV) changes. The dedicated model has been created that defines the dynamic of particular values changes (CV) according to described control by MV. The defined model is adjusted in online mode, which gives good performance for boiler control units. The general description of optimization methods are presented in **Table 2**. The calculation have been performed by using dedicated software SILO [56, 57], developed and implemented by Transition Technologies Company [58]. SILO collects various data to select the best solution for defined models and for current boiler conditions.

The main goal in boiler operating optimization was to improve the fire ball shape to achieve the nominal steam temperatures as well as low NO_{x}, CO, and O_{2}, which depend on flue gas distribution inside the furnace. The information about temperature distribution provided by the AGAM system has to be integrated with control units. Before optimal fire ball shape evaluation the dedicated tests have been performed to develop the relation between fire shape and MVs values.

The AGAM system gave the 12 independent values for 12 temperature zones in which the AGAM measuring surface was divided. The temperature values from 12 zones have been transferred into shape categories, which were described by use of three parameters listed below:

the position of temperature hot spots indicating the maximum temperature

intensity is evaluated as average temperature for a given temperature profile

the dispersion is calculated as standard deviation of AGAM temperatures for each zone

During the optimization process, the boiler operation conditions and combustion process parameters were analyzed relying on standard measurements (steam temperatures, emissions, etc.). Optimization in this phase was focused on finding the highest performance and optimal fire ball shape. The different shapes were evaluated at low and high boiler load. Finally, the monitored process parameters were related to the evaluated fire ball shapes. **Figure 15** presents the evaluated fire ball shapes at low and high load where the process performance was the best.

### 4.2. Combustion process analysis using CFD modeling

To simulate the OP-650 front-fired boiler operation a commercial CFD code ANSYS Fluent was used. Furnace exit gas temperature measured by the AGAM system, located approximately 4 m under the furnace exit, was one of the key values used for verification studies. The simulation results comparison were performed using average temperature as well as temperature profile.

#### 4.2.1. Numerical model of front-fired boiler

To create the three-dimensional geometry and the numerical mesh, ANSYS Design Modeler and ANSYS Meshing software was used [59]. Quality of mesh is one of the most important factors in the numerical simulation of the large-scale boilers [60]. A body-fitted mesh was created containing mostly hexahedral elements and mesh quality was evaluated based on the two parameters:

Orthogonal Quality—defining to what extent the mesh is not orthogonal. The best cells are close to 1. In the generated mesh minimum orthogonal quality was 0.2, with an average value of 0.95.

Aspect Ratio—defining the ratio of cell sizes in different dimensions. Presented mesh is characterized with maximum aspect ratio of 7 and with the average value equal to 1.4.

Because of the discrepancy of scale between the burner and of the larger volume of the boiler, the geometrical model was subdivided into fine-grid regions around the burners and coarser regions elsewhere. In order to facilitate the meshing process, the circular drop tubes have been replaced by rectangular ones and the platen superheaters have been modeled as zero-thickness horizontal planes. A mesh independence analysis was performed for the initial mesh consisting of 4.0 million control volumes. Reduction of cell number to about 3.3 million did not substantially change the predicted temperature field and this number of cells was selected to obtain compromise between solution accuracy and computational time. Computational domain and numerical mesh close to the swirl burners region is shown in **Figure 16**.

Simulations were conducted for a boiler load equal to 200 MWe (90% of nominal load). Two cases have been considered: before and after optimization process. Operating conditions have been collected from power plant online monitoring system by averaging 2-hour measurements during steady-state boiler operation. Air/fuel distribution, temperatures, and pulverizers’ activity are presented in **Table 3**.

#### 4.2.2. Results of numerical simulation

A tool being able to measure not only the average value, but also the temperature field is needed for a comprehensive CFD model validation. The aggressive environment, ash particles, high temperatures, and large dimensions of the boiler make temperature measurement a complex task [61]. Traditionally used suction pyrometry is extremely accurate and provides information of local temperatures only [62]. In order to account for a dynamic and turbulent environment in the combustion chamber, representative temperature distribution can be obtained only by performing simultaneous measurements. Number of available test ports poses a technical limitation. A turning point in accuracy assessment of CFD furnace models came with Acoustic Pyrometry [63]. This technology can provide average value in the selected cross-section. A detail validation study of the CFD model by comparison of computations against Acoustic Pyrometry measurements was demonstrated in [19].

To show the temperature field variation over the surface, an area-based and mass-based uniformity indexes (*UI*_{area}, *UI*_{mass}) of temperature distribution were used:

where subscript *i* is the denoting facet defining the surface, *N* is the total number of facets (12 zones for AGAM, 2260 facets for CFD), *T*_{average} is the average temperature in the cross-section plane (*K*), *A* is the cross-section area (m), *ρ* is the density (kg/m^{3}), and *w* is the velocity (m/s).

The mass-based average temperature in the cross-section plane is defined as:

For the calculation of mass-based uniformity index, the values of average densities and velocities in 12 AGAM zones were taken from CFD results as no such measurements were available. CFD prediction of temperature distribution parameters are presented in **Figures 17** and **18**.

UI values of one indicate the highest uniformity. The area-based UI calculated from measurements exceeds 90% for both cases, which means that no significant gradients exist. Measured standard deviation was 161°C and 164°C for 200 MWe before and after optimization, respectively. A conclusion is that that a more uniform temperature distribution was achieved after optimization process. The temperature homogeneity is important also for platen superheaters operation, which are located above the combustion chamber. High value of area-based UI indicates relatively uniform distribution of radiative heat flux. The mass-based UI of temperature distribution can be associated with convective heat transfer. Calculated values of 77% and 79% suggest less evenly distributed convective heat flux in comparison to the radiative one. The differences of temperature distribution inside the combustion chamber are presented in **Figure 19**.

After optimization the measured average temperature was 34°C higher than calculated temperature, which corresponds to 3% relative error where before optimization the difference was only about 2°C. Since the heat transfer in the furnace is highly dependent on the emissivity of the furnace wall, the choice of this value has significant influence on the calculated average temperature. All the results presented in this paper correspond to wall emissivity equal to 0.7. More detailed description of thermal wall boundary conditions can be introduced using deposition model with thermal properties submodel [55].

In addition to the verification of calculated temperature distributions, comparisons were also made of the basic values measured during boiler operation. Average values of CO, NO_{x}, and UBC were determined at the outlet plane of the model and compared with measured values before and after the optimization process (**Table 4**).

The results of numerical simulation and measurement confirmed the favorable effect of the combustion optimization process on the values received at the boiler outlet. By changing the operating conditions value of CO, NO_{x}, as well as the unburned carbon in ash was reduced.

### 4.3. Results of optimization process

The optimization performance evaluation was done by analyzing historical data collected from 2 months test. The method of analysis consists of results comparison before and after optimization. The results are presented as mean values of boiler efficiency (**Figure 20**) and UBC content (**Figure 21**) in function of unit load in range from low load (135 MWe) to maximum load (225 MWe). The measuring data has been prepared before analysis and only the representative data has been chosen to evaluate the results.

The results of performed tests indicate that boiler efficiency can be increased approximately by 0.2%, thanks to the optimization of boiler operation conditions. The biggest improvement of boiler performance can be achieved at high loads even by 0.6%. The lowest increase of efficiency is observed for temporary states when the time of boiler operating is shorter than at nominal load. All the effects presented above have been obtained keeping the NOx emission and total organic carbon content in ash at the correct level.

Thermal boiler efficiency depends largely on the steam parameters at the boiler outlet. In the case of exceeding defined value, spray attemperators are launched. This process reduces the overall efficiency of the boiler and the amount of spray water has to be minimized. If the steam temperature is below the required value, parameters of the working fluid at the turbine inlet are lower, reducing the overall process efficiency. The measurement results of live and reheated steam temperatures as a function of unit load are shown in **Figures 21** and **22**.

According to dynamic load changes profile and frequent change of mill unit configuration in operation those temperatures could differ from the reference value (540°C) to over 20–30°C. At low load (135 MWe) those temperatures could drop even below 520°C. At high load (225 MWe) the temperatures reach the 540°C level but with support of significant spray water flow. The results of steam temperature analysis (before and after optimization) confirm the positive effect of optimization. Temperature values of live and reheated steam are close to the required temperature 540°C in the full range of unit load. Also, the temperature variations of steam in transition states, between minimum and maximum load, were eliminated. Optimization of reheated steam temperature has a significant effect on the boiler efficiency. Due to the low pressure of reheated steam in comparison with live steam, temperature has large influence on the total energy amount generated inside the boiler.

The observed parameter during optimization process was temperature of live and reheated steam measured on the left and on the right side of boiler outlet, to underline that symmetry process also requires optimization. Optimization of symmetry process in practice means the minimization of differences between left and right side of the boiler in: coal flow, live and reheat steam temperatures, as well as flue gas temperature measured by AGAM and O_{2}, CO, and NO_{X} contents in flue gases. Thanks to the optimization process, the differences of steam temperature at the boiler outlet between left and right side were significantly reduced. This effect is observed for both live and reheated steam.

## 5. Conclusions

The results of experimental and numerical research of coal combustion process in front-fired pulverized coal boiler are presented in the paper. The experimental part was focused on coal combustion optimization to achieve better performances of electricity production process. The main goal of the work was the increase of thermal boiler efficiency based on the advanced monitoring technique of boiler operation conditions.

The experimental research was conducted on the front-fired pulverized coal boiler OP-650, located in EDF Polska Rybnik Power Plant. In the optimization process, the modern measurement systems of flue-gas temperature distribution and fuel-air distribution were used. The acoustic gas temperature measurement system allows to online monitor temperature distribution at furnace exit. The optimal combustion process was identified by uniform temperature distribution measured in online mode. To modify temperature distribution and properties of fire ball, the control system of coal and secondary air mass flow rate has been used. The comparison of results, before and after optimization process, indicates that the boiler efficiency can be increased approximately by 0.2% in the full range of unit load. The rise of boiler efficiency can reach up to even to 0.6% at the nominal load of boiler. Ensuring symmetry of the combustion process by online measuring of temperature distribution allows to reduce the differences between the outlet steam temperature at the left and right of the boiler.

To verify the impact of operating conditions on the temperature distribution as well as CO, NO_{x}, and TOC content in flue-gas, the numerical investigation using computational fluid dynamics modeling was done. To validate results of modeling, the acoustic gas temperature measurements in the form of temperature distribution at the furnace exit was used. The results of numerical simulations were compared with measured values. Good accuracy between numerical and experimental results was observed.

Results of numerical simulations conducted for the two states, before and after optimization, confirmed combustion process improvement by optimization of boiler operating conditions. Boundary conditions included in numerical models were based on values from power plant online monitoring system, before and after optimization process. Reduction of CO, NO_{x}, and TOC content in flue-gas by changing boundary conditions according to the measured values was confirmed. The numerical modeling provides also new important results of combustion process as, e.g., temperature distribution inside whole chamber of boiler.

Based on the results provided by modern measuring systems and complex CFD analysis, it becomes possible to effectively monitor the combustion process and indicate new guidelines (including changes in procedure of automatic control systems) to ensure the optimal operation of the boiler. The online monitoring systems allow to improve the combustion process through the uniform coal and air supply into the boiler. The effects of the change may increase boiler thermal efficiency keeping the NO_{x} emissions and TOC content at the correct level. The combination of these powerful tools can be used more often by the researchers working on the combustion process improvement as well as by the boiler manufacturers.