Validation of a thermal non-equilibrium Eulerian-Eulerian multiphase model of a 620 MWe pulverized fuel power boiler

The use of a thermal non-equilibrium Eulerian-Eulerian model for the simulation of a 620 MWe power boiler is proposed for capturing the combustion and radiative heat transfer found in the pulverized fuel systems. The models eliminates the use of a Lagrangian reference frame in tracking solid fuel particles thereby reducing the computational expense and time. The model solves the scalar transport for the particle mass, energy and radiation interactions between the pseudo-particle and continuous phases. The goal is to apply the modelling approach to generate a simulation database for different load cases and firing conditions which in turn will be used to study flexible operation. The model is validated against both numerical and applicable site data measurements. It is shown that the model is able to adequately resolve the furnace and superheater wall heat fluxes. Additionally the resolution of the flow field, combustion dynamics and wall fluxes are demonstrated for both an 80% and 60% operational loads. Moreover, it is shown that the Eulerian-Eulerian model results in approximately a 30% computational resource reduction when compared to traditional modelling approaches.


Introduction
The mid-merit operation of Coal-Fired Power Plants (CFPP) arises from the intermittent supply of electricity produced by renewable sources. Research into the flexible/mid-merit operation of power plants is needed to better understand the dynamic operations arising from these off-design conditions if renewable energy sources are to be integrated onto electrical grids.
Three-dimensional computational fluid dynamics (CFD) can capture the thermal performance of CFPPs at different load conditions with the necessary accuracy. These 3D CFD models include the modelling of combustion, radiation, and fluid flow phenomena. The traditional method for modelling solid-fuel combustion systems, such as CFPPs, is acknowledged as the Eulerian-Lagrangian (EL) approach. Successful implementations, using the EL approach include the following studies [1], [2] and [3], where the latter focused on the effect that the burner swirl direction has on the heat uptake to the furnace and superheater walls of a 620 MWe boiler.
Another approach is to use a Eulerian multi-fluid model as used by the following studies [6], [4] and [5]. Knaus et al [7] compared the Eulerian-Eulerian (EE) and EL approaches for a 550 MWe coal-fired utility boiler. The predicted combustion products yield only small differences, with carbon monoxide being better predicted when using the Eulerian-Lagrangian approach. The computational effort of the Eulerian approach proved more efficient and faster convergence times were noted. It should be stated, that these studies were conducted using in-house CFD codes which is difficult to apply by engineers and scientists outside the research group.
The current work proposes the use of an EE modelling strategy that incorporates thermal non-equilibrium between the gas and particle phases which allows for adequate particle concentration tracking and heat transfer interactions, with a reduced computational cost when compared to the EL simulation approach. The model is validated for a full scale utility power boiler for various load cases. The validated modelling technique is to be used for future studies in developing a CFD data-driven surrogate model of a 620 MWe CFPP's furnace and radiant superheaters.

Numerical modelling methodology
An outline of the modelling methodology, utilized by the proposed EE model, is provided in this section highlighting the general conservation equations, pseudo-solid phase tracking, and the combustion and radiation modelling techniques employed. The basis of the model assumes that the gas and solid phases behave as inter penetrating continua with both phases being described in a Eulerian reference frame. A further assumption is that the fluid is in dynamical equilibrium with the solid phase, (similar to the work of [6]), meaning that velocities of both phases are locally the same.

Transport equations
The time-averaged steady-state continuity, momentum, k-turbulence and energy equations for the gas-phase can all expressed using the following generalised transport equation for a variable φ expressed as: where ρ is the gas density and Γ φ is the diffusion coefficient for variable φ.
In the present work the flue gas phase is modelled using the species-transport approach where a transport equation is solved for each constituent mass fraction in the gas mixture. The governing equations for the gas phase can be written as follows; Mass conservation: 2 MATEC Web of Conferences 347, 00004 (2021) SACAM2020 https://doi.org/10.1051/matecconf/202134700004 Momentum conservation: The momentum equation is closed by approximating the Reynolds stresses using the Boussinesq relationship as shown in equation (4).
This approximation is used in conjunction with the realizable k-turbulence model of Shih et al [8] Energy conservation: where E in the total energy term, P the pressure, λ is the effective thermal conductivity and T g is the gas temperature.
Species transport: where k is the number of species present in the domain and J k is the species diffusion flux. r R j,k is the net rate of production of species j due to r reactions.
To correctly account for the particle inertial effects on the gas phase convection an effective density is defined as follows and is used in the above equations.
In equation (7), ρ p is the particulate matter density and φ mp [kg p /kg] is the local mass concentration of particulate matter present in the gas volume.

Solid phase modelling
The pseudo particles transported into the domain are modelled using a scalar field and transported through the computational domain according to equation (1). The pseudo particle scalar fields are used to define the fuel characteristics based on the proximate analysis composition (refer to Table 2), namely consisting of moisture, volatile matter, fixed carbon and ash. Each of the scalar field equations solved for are given in Table 1. Original/initial mass of particles Enthalpy of particle Equation (8) The energy transport of the pseudo particle, is transported by defining the particle enthalpy using the following equation: The equation accounts for all the processes associated with energy transport to the particle, namely convection Q conv , radiation Q rad , latent heat dM evap dt h f g and near surface char oxidation f heat dM c dt h rxn . This gives the model the ability to track the particle temperature in the domain, moving the model away from the thermal equilibrium approach incorporated by previous studies using an EE approach ([6], [7]). The particle temperature is important in describing the sequential steps found in modelling combustion processes, especially at low boiler loads where mixing and ignition become problematic.

Combustion modelling
The combustion modelling of coal comprises four sequential steps, namely the heating/cooling of the particles, evaporation, devolatilization and char oxidation. The process starts with the pseudo-particles entering the domain and heating up until the moisture's boiling temperature of 373 K is reached [9]. The realization of the boiling temperature would release moisture based on a boiling rate which is defined as: Following this is the devolatilization process whereby the volatiles present in the fuel are liberated. The initialization temperature of the devolatilization process was found to vary between 553 − 623 K for a variety of coals by Ranade and Gupta [10], thus particle heating would need to ensue after the moisture evaporation to bring the particles to the initialization temperature. A single-step Arrhenius kinetic rate model is used to approximate the devolatilization process since these models can predict the devolatilization of coal and are simpler to implement, according to Sankar et al [11]). The rate is expressed in equation (10) as: where and E a,vol = 6.7E7 [J/kmol] from the works of Sheng et al [12]. Here m vol is the mass of the volatile released during combustion and m 0,vol is the initial volatile mass present before combustion.
The final step of the solid fuel combustion process is realized when all the volatiles have been driven off from the pseudo-particles. The product species of the char oxidation process is CO, with 100% of the resultant heat being absorbed by the pseudo-particle contributing to the energy rise of the particle matter in a cell. The char oxidation rate is derived from the works of Baum and Street [13] and Spalding [14]. The model is commonly referred to as the diffusion-kinetics limited model, written as: The interested reader is directed to [15] for an in-depth derivation of the above equation. The kinetic parameter R c = A c exp(−E a,c /RT p ), where A c = 0.0053 [kg/m 3 ] and E a,c = 8.37E7 [J/kmol] from the works of Sheng et al [12] and Hanjalić et al [16]. Here A p is the particle surface area, R di f f = (5E − 12)(0.5(T p + T g )) 0.75 )/d p is the diffusion rate coefficient, which is a function of the gas and pseudo-particle temperature and average diameter, and P O 2 is the partial pressure of oxygen.
According to Sankar et al [11], the gas phase reactions in boilers are fast, with the chemical time scales being orders of magnitude faster than turbulent mixing time scales. This leads to a mixed-is burnt assumption [17], which assumes the chemistry is infinitely fast and irreversible, being adequate for species and temperature distribution predictions. The gasphase volatile and CO reactions are approximated using the two-step global reaction shown in equation (12).
This study makes use of the two-step finite-rate/eddy-dissipation model (FR/EDM) in ANSYS Fluent® 19.3 [15] to calculate the gaseous reaction rates. The model calculates three rates and uses the minimum of the three values for the source term calculations for the species transport equations and subsequently for the energy release source terms. Chemical reaction rate: Rate of turbulent production eddies dissipation: Rate of reaction eddies dissipation: In the above equations Y p is the product species mass fraction, Y R is the mass fraction of any reaction species, A and B are model constants, M w,k is the molecular weight of j th species, ϑ k,r is the stoichiometric coefficient of reactant k in reaction r, [C l,r ] m,r is the molar concentration of species l in reaction r with rate exponent m, A r = 2.56E11 [s −1 ] and E a,r = 1.081E8 [J/kmol] taken from the study of He et al [18].

Radiation modelling
The heat transfer to the furnace, platen and final superheater walls of the validation case are mainly due to radiation. The radiation transport for the current work was resolved using the P1 radiation model which includes the effects of particle absorption (α p ) and scattering (σ p ) as well as gas mixture absorption (α g ). The P1 model transport variable is the incident radiation (G -W/m 2 ), and can be written for a particle laden domain as: The gas mixture absorption coefficient was calculated using the weighted sum of gray gases model (WSGGM) which accounts for the tri-atomic gases CO 2 and H 2 O as proposed by Smith et al [19]. The solid phase radiative properties (namely α p , σ p and E p ) were determined using an equivalent Eulerian description for each expressed as: (19) where N p = ρφ mp0 V/ ρ p πd p 3 /6 , is the number of particles present based on the average diameter pseudo particle sized p , A pro j = N p πd 2 p /4V is the projected area, σ S B = 5.67E − 8; [W/m 2 K 4 ], T p is the particle temperature and V is the cell volume. The particles emissivity ( p ) and scattering factor ( f p ) were implemented using a variable property correlation as used in the work of Laubscher and Rousseau [20].

Numerical setup and model inputs
The modelling methodology, described previously, was successfully incorporated in the modelling of a 620 MWe subcritical power boiler. Figure 1 illustrates the computational domain used in this validation study highlighting the various heat receiving surfaces and inlets. The burner feeds the combustion chamber with primary and secondary air through the burners inner and outer annuli respectively. At 100% maximum continuous rating (MCR) the primary air (PA) supplies the fuel and air at a temperature of 373 K and secondary air (SA) enters at a temperature of 577 K. The fuel characteristics used for CFD modelling purposes are given in Table 2. The plant in designed to operate at a 100% MCR with a excess air ratio of 15.5%. As the plant is turned down from 100% MCR the excess air ratio is increased to ensure sufficient air and fuel mixing occurs. At MCR ratings of 81% and 60% the excess ratios were calculated as 20.9% and 26.3% respectively, based on the design schedule data. The inlet mass flow rates and temperatures for the plant's MCR ratings of 100%, 81% and 60% are listed in Table 3.

Numerical solution strategy
The current simulations were performed using the ANSYS Fluent ® 19.2 pressure-based solver. Pressure-momentum coupling was set to the SIMPLE technique, with momentum, energy and species equations being discretized using the second-order upwind method. The pressure equation was discretized using PRESTO!. This setup holds for both the EE model and the detailed EL model. The difference arises in the use of scalar fields in the EE configuration and using the discrete phase modelling approach in the EL setup. Simulations were performed on a numerical mesh consisting of roughly 6.2 million cells. A mesh independence study was conducted for mesh sizes consisting of 4.2 and 10.2 millions with the 6.2 million cell model deemed acceptable for this study. To ensure an accurate numerical solution the mesh aspect ratio was kept below 15 and the minimum orthogonal quality above 0.2. The discrete phase equations, for the EL model, were solved once every 30 fluid phase iterations. The number of particles injected per burner was set to about 7800, totalling to 140 000 particles in the entire domain. To ensure a stable converged solution the spatial discretization for all fields were set to first order upwind (except pressure) and solved for 1500 iterations before increasing the discretization order. The simulations were then run for a further 10000 iterations. For all cases the maximum mass imbalance was 0.046 kg/s for a total gas flow rate of 376 kg/s. The maximum heat imbalance was 2450 kW for a total heat input of 855 MW.

Results and discussion
The discussion herewith compares the results of an EL numerical model with the results of the developed EE modelling approach, (discussed in section 2) for a 100%, 81% and 60% MCR loads cases. Overall a 30% (9 hours to approximately 6 hours) decrease in simulation time was observed across the simulated load cases. Table 4 highlights the relative errors obtained between the EL and EE numerical models. A maximum error of 8.15% occurs at the platen superheater for the 100% load case. The relative errors are deemed acceptable for the decrease in computational time required by the EE model. The graphs shown in figure 2 illustrate the calculated (based off operational plant measurements), EL and EE models results, for the heat loads to the furnace, platen and high temperature superheater (HTSH) walls.  The general accuracy of the EE model is seen to be able to sufficiently capture the overall heat loads when compared to the calculated and EL model results. A notable difference is seen in the 60% load case of figure 2 (c), where both the EL and EE numerical models underpredict the platens heat load. Figure 3 compares the gas velocity fields for the EL and EE model. It can be seen that results are in good agreement. The EE model tends to under predict the lower burner velocity transport resulting in combustion occurring closer to the burner mouth and lower furnace walls. This effect is seen in both the temperature contour plots of figure 4 (d-f) and wall heat flux contours of figure 5 (d-f).
The temperature contour plots of figure 4 (d-f) illustrate the EE model's ability to sufficiently resolve the temperature field when compared to the EL model. The lower burners tend to initiate combustion closer to the burner leading to high temperatures. Due to the lack of velocity of the gas phase in this area, the high temperatures are closer to the walls of the furnace leading to a higher heat flux observed as seen in figure 5 (d-f). In general the temperature and velocity fields are deemed sufficiently accurate for the purposes of surrogate model development, since the parameter of interest is the resolution of the wall flux field. The heat flux profiles are given in figure 5 (a-f). The EE model over predicts the wall fluxes in the lower half of the boiler, namely due to the high temperatures experienced in these zones in and around the bottom burners. The platen and HTSH radiant superheaters show similar flux profiles with the top half of the furnace illustrating sufficient agreement between the EL and EE models. The spatial distribution of the particles (expressed in kg/m 3 ) shown in figure 6 (a-f) highlights the EE model's ability to resolve the particle concentration throughout the domain for the various load cases. The general trend of the EE model is with sufficient accuracy when compared to the EL model.

Conclusion
The validation of a full scale 620 MWe boiler was conducted using a thermal nonequilibrium Eulerian-Eulerian model for three load cases, namely a 100%, 81% and 60% steady-state loads. The validation cases include a comparison with results obtained by a numerical model using a Eulerian-Lagrangian approach and where applicable measured site data. The developed EE model demonstrated adequate performance in predicting the flow field, wall heat fluxes and fluid property distributions throughout the domain.The relative accuracy of the EE model ranged from 2 to 8% for key parameters used in the study. The computational speedup, of 30%, by using the EE model is a beneficial attribute when considering the EE models intended future use in the development of a CFD data-driven surrogate model.