Numerical analysis of pulverized biomass combustion

Combustion of pulverized biomass in a laboratory swirl burner is computationally investigated. The two-phase flow is modelled by an Eulerian-Lagrangian approach. The particle size distribution and turbulent particle dispersion are considered. The radiative heat transfer is modelled by the P1 method. For modelling turbulence, different RANS modelling approaches are applied. The pyrolysis of the solid fuel is modelled by a single step mechanism. For the combustion of the volatiles a two-step reaction mechanism is applied. The gas-phase conversion rate is modelled by the Eddy Dissipation Model, combined with kinetics control. The results are compared with measurements.


Introduction
For the generation of power and heat, combustion is being used as the major process since many decades [1][2][3]. For solid fuels, gasification techniques are additionally used [4].
In parallel to the efforts of utilizing renewable energies [5,6], as well as recovery techniques [7] coupled with heat transfer enhancement methods [8,9] combustion continues to play an important role in renewable energies. This is due to the fact that the biomass, which is a renewable fuel, is also converted by combustion [10,11]. Thus, the energetic utilization of biomass via combustion process is the main focus of the present contribution.
In large and medium scale utility boilers for centralized power generation, designed to fire pulverized coal, the common way of burning biomass is co-combustion [12]. In utility boilers, usually only a rather small portion of the energy feed (<50%) is provided by the biomass. This limited use of biomass is caused, on the one hand, by the different fuel properties of biomass, especially with respect to its ash that can lead to increased corrosion problems (which are not addressed in the present study). On the other hand, the logistic reasons play here a role, as the available biomass in the catchment area of the power plant with affordable transport costs limits the extent of biomass usage.
In a computational analysis, in general, validation is a central issue. In combustion problems with very complex physics, leading to a very complex mathematical modelling, the validation is particularly important. This is especially critical for non-gaseous fuels that exhibit an even higher complexity due to the two-phase flow.
Pulverized coal combustion was investigated by many researchers in detail over several decades [13][14][15][16][17][18], so that and numerous validation studies have been performed. On the other hand, compared to pulverized coal, the pulverized combustion of biomass and its modelling have not yet been investigated in that detail. The present work aims to contribute to the validation of computational models for the modelling of pulverized biomass.
Large-scale applications are not very convenient for detailed model validation studies, due to the lack of measurements on the one hand (as they are difficult to perform), and the huge computational effort required on the other hand. Thus, for validation, the measurements obtained at the test facility located at RWTH Aachen University, consisting of single swirl burner [18] are taken as the experimental data base.
Flames measured at this test facility were computationally investigated by different researchers in the past. A 100kWth pulverized coal flame was experimentally and computationally analyzed by Toporov et al. [18] applying the k-ε turbulence model. In a more recent study, Warzecha and Boguslawski [19] investigated a 100kWth pulverized coal flame on the same test rig computationally by applying the k-ε model as well as LES. Further LES and k-ε model based analyses of the same flame were presented also by Sadiki et al. [20], later on. A computational study of the same flame was presented by Gaikwad et al. [21], who applied the SST turbulence model.
The present analysis differs from the abovementioned investigations mainly in two aspects. Firstly, the biomass combustion is investigated, which was not addressed in the previous studies mentioned above. Secondly, flames with 40kW thermal load are considered, whereas larger burners were considered in the previous studies. The larger burners of the same test rig incorporate a tertiary air injection, whereas the tertiary air is missing in the 40kWth burner.

Modelling
The general-purpose Computational Fluid Dynamics (CFD) code ANSYS Fluent 18.0 [22] is used, which utilizes a finite volume method of discretization. The density of the gas mixture is calculated assuming an ideal gas. Temperature dependence of material properties is considered via fourth order polynomials [23]. A coupled solver is used, and a second-order upwind scheme is utilized to discretize the convection terms [22]. The gradient computation technique was least squares cell based. Stabilization was achieved by a standard cell to face slope limiter.

Two-phase flow and convective transport
A Eulerian-Lagrangian approach is adopted, as it is the built-in standard formulation of the employed software, for pulverized coal combustion. The gas phase and particle phase equations are solved alternately, where, the particle iterations are performed after each 40th gaseous phase iteration. The volume occupied by particles, and particle-particle interactions are neglected. Only the drag force on particle is considered, assuming a spherical shape [24].
The size distribution is assumed to follow the Rosin-Rammler distribution [25]. A uniform particle temperature is assumed calculating the convective heat transfer coefficient using the correlation of Ranz and Marshall [26].

Turbulence
The RANS approach is used to model the turbulent gas flow. Within this framework, different turbulence viscosity based turbulence models are considered such as the Standard k-ε model [27] (S-KE), the Realizable k-ε model [28] (R-KE) and the Shear Stress Transport k-ω model [29,30] (SST). The near-wall turbulence is modelled by the standard wall-functions approach [27].
For the turbulent diffusion of the scalar quantities, the gradient-diffusion approximation is used assuming constant Prandtl-Schmidt numbers for the prevailing fully developed turbulent flow (0.9 for the energy, 0.7 for the species transport equations). The effect of gas turbulence on the particle motion is modelled by the socalled "discrete random walk" model [31], whereas the influence of the particle phase on the gas turbulence is neglected.

Radiative Transport
The radiative heat transfer is modelled by the P1 model [22]. The absorption coefficient of the gas mixture is calculated using the Weighted Sum of Gray Gases Model [32], assuming an equivalent path length for the domain.
The particle phase radiation is taken into account assuming the value of 0.9 for particle emissivity and scattering factor. The walls are assumed to reflect diffusely. The wall emissivity is assumed to be 0.7.

Combustion Modelling
The fuel particle experiences an evaporation and pyrolysis with increasing temperature. The residual char burns via heterogeneous reactions, as the combustible volatile matter reacts homogeneously in the gas phase.

Pyrolysis
During the pyrolysis, the swelling of the particles is accounted for, with an assumed swelling coefficient of 1.4. The combustible volatile matter is represented by a molecule CxHyOz, assuming a molar mass of 30 kg/kmol, where x, y, z depend on the elementary analysis of the fuel.
Following Badzioch and Hawskley [33], a firstorder, single-rate pyrolysis is assumed, where the rate coefficient is expressed by an Arrhenius rate expression

Char oxidation
Char is assumed to oxidize to carbon monoxide in a single-step irreversible heterogeneous surface reaction. The rate coefficient is calculated considering a combined rate limiting effects of kinetic and diffusion processes [34,35].
The kinetic rate is described by an Arrhenius rate expression. The diffusion rate is calculated as function of particle size, boundary layer temperature and diffusion coefficient.
It is considered that the released heat by combustion is partially absorbed by the particle itself. In the current study, heat of reaction absorbed by the particle is assumed to be 30% in ratio.

Gas phase reactions
The combustion in the gas phase is assumed to occur via a global reaction scheme comprising two irreversible reactions [23]. In the first reaction, the volatile matter is assumed to react to CO and H2O. The second reaction is the oxidation of CO to CO2. The effect of turbulence is considered by a rather simple approach. The resultant time-averaged volumetric species conversion rate is assumed to be limited by the smaller one of the kinetic and mixing rates. The latter, i.e. the rate of mass transfer to smallest scales via dissipation of turbulence eddies, is modelled by the Eddy Dissipation Model of Magnussen and Hjertager [36] using the original model constants.

Test Case and Fuel Specific Definitions
As mentioned above, the considered flames were obtained and measured at a test rig of RWTH Aachen University [37]. Sketches of the test rig and the swirl burner are presented in Figure 1.
The torrefied biomass flame has a thermal power of 40 kW. The operating conditions [37] are summarized in Table 1. It shall additionally be noted that the oxygen  mole fraction of the primary stream is 18.8% [37] (instead of 21% that is usual for air as injected as the secondary and staging stream). The pulverized fuel is injected through the primary nozzle. A swirling flame is achieved by swirling the secondary air. The geometric swirl number of the secondary air is given as 0.958 [37]. The ultimate and proximate analysis [37] of the fuel is given in Table 2. Please note that the sulphur content is neglected in the present calculations.
The measured particle size distribution [37] is provided in Table 3. The provided particle size distribution [37] does not necessarily have a high resolution. In the calculations, a Rosin-Rammler distribution is fitted to the measured particle size distribution, where the minimum and maximum particle diameters are assumed to be the half and the double of the minimum and maximum values given in Table 3, respectively. Totally 10 particle size classes are considered. For taking the turbulent particle dispersion by random-walks [31], the latter is taken into account by 10 trials.
A critical parameter that has a rather strong effect on the prediction of pulverized fuel combustion is the pyrolysis rate, since it can strongly depend on the fuel type and its determination is afflicted with rather high uncertainties. Within the framework of single-rate pyrolysis assumption that is approximated by an Arrhenius rate expression [23], the used rate constants are presented in Table 4.  For biomass, the experience and the data base is not that broad compared to coal, as there are many different types of biomass and the pyrolysis behavior strongly depends on the biomass type. In the present case of torrefied biomass, the constants by Tolvanen et al. [38] are used ( Table 4) that are suggested for a corresponding category.
For char oxidation, due to the lack of more detailed information, the same rate constants that are usually used for coal are used for biomass, and are presented in Table 5. The rate constants for the kinetics of the gas phase reactions (1st reaction: volatiles oxidation to CO and H2O, 2nd reaction: CO burn-out) are provided in Table 6.
These are the default values suggested by the used software. They represent a modified set of the empirical constants suggested by Dryer and Glassmann [39] and Westbrook and Dryer [40] for two-step oxidation of hydrocarbons.
It should be admitted that the rate constants for the first reaction (Table 6) are rather arbitrary, since the assumed, hypothetical volatile molecule structure as CxHyOz is not necessarily corresponding to the hydrocarbons underlying the empirical constants. Nevertheless, it is currently assumed to be a reasonable assumption to take the chemical kinetics effects on the reaction at least approximately into account, in the absence of more accurate information. The same rate constants are used for the oxidation of the volatile matter from coal and biomass.

Solution domain, boundary conditions, grid
As the geometry and boundary conditions are essentially axi-symmetric, the time-averaged solution of the problem needs to be axi-symmetric. Thus, as a first approach within a RANS formulation, a twodimensional, axi-symmetric analysis is applied.
On the other hand, it is also quite well known that the turbulence models within a RANS formulation may fall short in capturing the effect of coherent structures on the time-averaged solution, so that an Unsteady RANS or LES approach is needed for better accuracy, which is, however, not yet attempted at the present stage.
As seen in Figure 2, the length of the solution domain is defined to be shorter than the total length of the device (Fig. 1), since the processes in the remaining length do not have a substantial influence on the upstream region, where the measurements were performed. This shortening of the domain length is verified by preliminary calculations, and also in agreement with the approach of Sadiki et al. [20] in the analysis of a different flame in the combustor. It can also be seen in Figure 2 that the domain is closed by a wall, defining a purely radial outlet. In our previous studies, this approach was found to be an optimal choice for strongly swirling flows that may exhibit a sub-critical state at the outlet, and its suitability has also been verified in the present study by preliminary calculations.
At the inlets (primary, secondary, staging) block profiles for all variables are applied as boundary conditions, prescribing the velocity components, temperature and composition according to the operating conditions provided in Table 1. The boundary conditions of the turbulence quantities are derived based on an assumed turbulence intensity of 4% and the local hydraulic diameter of the considered inlet. The particle velocities and temperatures are assumed to be equal to those of air at the inlet. At the outlet, a constant pressure is applied along with zero normal gradient conditions, whereas, symmetry boundary conditions obviously hold at the symmetry axis. At walls, no-slip conditions hold, amended by the standard wall functions [27].
For the energy equation, a constant temperature of 800°C is applied as boundary condition as indicated in Zabrodiec et al. [37]. The emissivity of furnace walls is assumed to be 0.7, whereas an emissivity of particles is assumed to be 0.9. Zero normal gradient boundary conditions apply for the species transport equations at the walls.
The grid is generated by a block-structured strategy using quadrilateral finite volumes. The grid resolution is determined by a grid-independence study. A grid is sought that is sufficiently fine to assure a satisfactory grid independence, but not necessarily finer. The resulting grid that is consisting of about 20,000 cells is displayed in Figure 3 (outlet section is not displayed).

Results
The predicted temperature fields by different turbulence models in the burner near-field are depicted in Figure 4. One can see that the cold mixture enveloped by the reaction zone (flame brush) extends into the furnace with a simultaneous gradual radial displacement due to the effect of the centrifugal force generated by swirl, resulting in a diagonally extending "tongue" like shape of the flame brush. The temperature rise along furnace centerline is due to the internal recirculation zone created by the swirl and vortex breakdown, which recirculates the hot exhaust gases back to the burner and offers a reactor zone for the pyrolysis of recirculating particles of the smaller size classes.
These features are predicted qualitatively by all turbulence models. However, there are also differences in detail. Compared to S-KE, the flame brush predicted by R-KE more slender and longer in shape and has a smaller radial extension. A slightly thicker "tongue" and a delayed temperature rise along the axis are predicted by SST, compared to S-KE. In the prediction by SST, the strong interaction of the flame brush with the staging air is noticeable. The wall jet created by the staging air separates quite substantially, due to the injector effect of the flame jet, and interacts with it. Such a separation of the staging air wall jet is not observed in S-KE prediction. In R-KE prediction, a separation is indicated, but at a much smaller extent compared to SST.
Predicted fields of volatile matter mass fraction by different turbulence models in the burner nearfield are presented in Figure 5. One can see that the main volatiles release is in a region beyond the outer shear layer of the flame brush, predicted, however, in differing detail  patterns corresponding to predicted temperature fields. It is interesting to note that there is a local volatile matter release zone in the burner quarl according to S-KE and R-KE predictions. This is due to the recirculation flow mentioned before and appears to be more intense for R-KE due to the stronger recirculation of the latter. In the SST exhibiting a weaker recirculation, an isolated volatilization region in the burner quarl is not observed.
The predicted fields of O2 mass fraction by different turbulence models in the burner near-field are presented in Figure 6. Here, the interaction of the swirling burner jet and the staging air wall jet can more clearly be observed. One can see that the wall jet massively separates in the SST prediction, but remains attached to wall in S-KE. In the R-KE predictions, the separation of the wall jet is weaker compared to SST.
The predicted fields of CO mass fraction by different turbulence models in the burner near-field are presented in Figure 7. It is interesting to see that the CO fields are just enveloping the regions of volatile matter, due to the CO generation by the first reaction of two-step combustion reaction scheme assumed for volatiles combustion in the gas phase. CO is also produced by the heterogeneous char oxidation reactions, which seems to follow similar patterns since they are in a way following the devolatilization.
The predicted distributions of CO2 mass fraction by different turbulence models in the burner near-field are presented in Figure 8. CO2 is rapidly depleted by the second reaction of the combustion scheme assumed for volatiles combustion, to lead to practically fully burnt values in the downstream. The role of the staging air to achieve burnout is also observed.
Radial profiles of axial velocity at three axial stations are shown in Figure 9, where measurements are compared with the predictions. In all three axial stations, the internal recirculation zone with substantial negative velocities can be observed in the central region, which is neighbored by the forward flow with high positive velocities, the radial position of which is gradually increasing in the downstream direction. It can be seen  that the predictions agree rather well with the experiments, as far as the radial extension of the internal recirculation zone is concerned.
However, the backflow velocities on the centerline at x/d=0.3, and the peak velocities of the forward flow at all axial three stations are underpredicted in magnitude by all of the turbulence models. Especially, the underprediction of the peak velocity of the forward flow is observed to be quite significant. In this respect, the performances of the different turbulence models do not differ from one another very much. Still, one can see that R-KE performs slightly better than S-KE model and SST is performing slightly better than the R-KE.
Radial profiles of the swirl velocity at three axial stations, comparing the experimental results with the predictions are shown in Figure 10. One can see that the size of the vortex core is overpredicted and the peak swirl velocities are underpredicted substantially by all turbulence models. In this respect, the qualitatively best agreement is again provided by SST, where R-KE shows a performance between SST and S-KE, whereas the difference to S-KE is moderate.

Conclusions
Pulverized biomass combustion in a laboratory swirl burner is computationally investigated. For turbulence, three RANS models (S-KR, R-KE and SST) are applied. Results are compared with measurements. The SST model is observed to perform qualitatively better than the other models, where a slightly better performance of R-KE compared to S-KE is also observed. Still, the performances of the models leave much to be desired, which calls for higher order models such as LES.