Simulation of Biomass Combustion with Modiﬁed Flue Gas Tract

: The combustion of biomass is accompanied by the formation of particulate matter, the presence of which in the atmosphere harms human health. It is important to show the issues of reducing these pollutants and their impact on human health. This article focuses on the process of biomass combustion. The used model consists of two parts: the combustion chamber and the ﬂue gas tract. The article shows four types of modiﬁcation of the ﬂue gas tract designed to reduce the amount of particulate matter in the atmosphere. Bafﬂes are located in the ﬂue gas tract, which is designed to capture the particulate matter. The ﬁnal model is simulated by turbulent–viscosity models, k- ε realizable model, and then k- ω shear stress transport model. The interaction between turbulence and chemical reactions is expressed by using the Eddy Dissipation Concept model. The results then show different proﬁles of temperature, velocity, and particle distribution. Based on the evaluated data from two different calculations, it can be concluded that the bafﬂes have a signiﬁcant effect on the reduction of particulate matter in the atmosphere. The used bafﬂes are able to capture mainly particles with a diameter greater than 100 µ m. A signiﬁcant number of particles with a diameter lower than 100 µ m ﬂows from the ﬂue gas tract to the surrounding environment.


Introduction
Air quality is the most significant factor that affects the function of the ecosystem and the health of people. Air pollution is linked to a broad spectrum of acute and chronic illnesses. Based on the literature review, it is obvious that humans are exposed to various complex mixtures of particulate matter (PM). These particles are able to infiltrate into the lungs, some may even get into the bloodstream and contribute to rise in cardiovascular, respiratory problems, lung cancer, and heart disease [1]. Therefore, it is necessary to reduce the concentrations of PM getting into the atmosphere in every area of their formation in order to decrease environmental pollution [2,3]. The PM mainly comes from road transport, power plants, domestic heating (mainly wood burning), and industry [4,5]. Long-term exposure to high levels of these particulates has been associated with a rise in morbidity and mortality. In many European countries, problems with particulate matter concentration in ambient air, which also contribute to the deterioration of air quality were recognized. It is very important to trap the solid particles in a device before their release into the air.
In order to minimize the production of particulate matter (PM), different filters and separators are added to the boilers, but the installation and operation are costly for the owners and are often difficult to maintain.
Therefore, solutions are sought that will meet the requirement to reduce the concentration of PM, while at the same time being financially least demanding. Such a solution is also a modification of a fireplace construction, which has an effect on the discharged amount of PM and will increase the efficiency, and the lifetime of filters, electrostatic precipitators, etc. The construction can affect both the amount of produced particles (in the combustion chamber) and the concentration of trapped PM (in the flue gas tract-gas path). The results of Sulovcova et al. [6] indicate that particles were reduced in the optimized flue gas tract, where the change of the flow direction appeared. The deflector placed in the optimized flue gas tract influenced the particle stream and PM larger than 1.7 µm remained in the settling area. In the work of Menghini et al., the optimization of heat exchanger geometry in the fireplace was investigated. The two blades were replaced with one pipe in the construction of the fireplace, but no reduction of emissions was noticed [7]. The simulation of particle trajectories and deposit formation were analyzed in the work of Borello et al., where their results showed that the entire furnace was exposed to particle deposition, and the most deposits formed close to the outlet of the flue gases and inlet of fuel in the furnace [8].
Prokhorov et al. investigated the optimization of boiler construction and their work dealt with the effect of high aerodynamic drag of the gas path on boiler power. Numerical simulation showed that the proposed structures would provide an aerodynamic drag of less than 640 kPa and increase the boiler efficiency. Aerodynamic calculation of the boiler was performed in order to minimize aerodynamic losses, but there was no association to PM reduction [9].
There is a lack of knowledge about the motion of particle flow in the construction of heat sources, and especially in the flue gas tract. In terms of behavior, the velocity distribution of particles in the boilers is particularly important, which are entrapped by the flue gases. Computational fluid dynamics (CFD) are very helpful tools for designing and optimizing processes. CFD provides many ways to predict combustion phenomena with some accuracy [10], and it is a relatively simple way to optimize solutions for various physical processes. Turbulence is a very demanding process and sensitivity of the turbulence model can influence CFD predictions. Based on Stephen B. Pope, it has been shown that in shear layers with strong pressure gradients, the k-ε model is quite unreliable, while the k-ω model is much more accurate when calculating multiple flows [11]. Kumar and Ghoniem illustrated that the k-ω shear stress transport (SST) model achieved significantly better results in swirling flows compared to the k-ε model [12,13].
On the basis of these and other knowledge, it is possible to modify a boiler to decrease the particle discharge into the chimney and to the surrounding environment. The project results will determine the impact of geometrical optimization of the flue gas tract on the capture of PM. A CFD model of particle flow will be created and the results will show the effect of the nature of the flue gas flow on the particle motion. Modeling of flue gas tract modification in CFD can help to better understand the behavior of particles in the boiler and contribute to the reduction of PM. This could be achieved by the proper position of baffles in the flue gas tract. In this work, the k-ε model, the k-ω turbulence models were used and analyzed.

The Description of Model
The model used in the simulation is shown in Figure 1. This model is a modified wood stove, composed of two parts. The first part shows a combustion chamber with a retort burner and holes for the air supply. This part represents a prism with a length of 400 mm, a width of 420 mm, and a thickness of 400 mm. The second part shows a flue gas tract with a deflector and three baffles. It has a length of 500 mm, a width of 420 mm, and a thickness of 400 mm. These components are designed to capture solid pollutants such as particulate matter. Solid pollutants are particles consisting of carbon, ammonium, metals, organic materials, nitrates, and sulfates. as particulate matter. Solid pollutants are particles consisting of carbon, ammonium, metals, organic materials, nitrates, and sulfates. The position of the baffles placed in the flue gas tract was also observed and the best variant was used in the complete model. The proposed variants are shown in Figure 2. The baffles were placed perpendicular to the vertical plane in variant 1. In variant 2, baffles were also perpendicular, but the last third of the baffle was turned down by an angle of 30°. Variant 3 had perpendicular baffles, and the last third of the baffle pointed vertically downward and then turned by an angle of 65° from the vertical plane. In variant 4, they were parallel to each other and rotated at an angle of ±20° from the horizontal plane. The used material in the complete model was glass on the front wall of the model and then fireclay panels on the other walls. The air supply was realized by using three types of inlets. The supply of primary air was realized through the retort burner for primary combustion. The second type of air supply was provided through holes on the rear wall of the model for secondary combustion. The supply of the last tertiary air was used to blow the glass. The position of the baffles placed in the flue gas tract was also observed and the best variant was used in the complete model. The proposed variants are shown in Figure 2. The baffles were placed perpendicular to the vertical plane in variant 1. In variant 2, baffles were also perpendicular, but the last third of the baffle was turned down by an angle of 30 • . Variant 3 had perpendicular baffles, and the last third of the baffle pointed vertically downward and then turned by an angle of 65 • from the vertical plane. In variant 4, they were parallel to each other and rotated at an angle of ±20 • from the horizontal plane.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 3 of 11 as particulate matter. Solid pollutants are particles consisting of carbon, ammonium, metals, organic materials, nitrates, and sulfates. The position of the baffles placed in the flue gas tract was also observed and the best variant was used in the complete model. The proposed variants are shown in Figure 2. The baffles were placed perpendicular to the vertical plane in variant 1. In variant 2, baffles were also perpendicular, but the last third of the baffle was turned down by an angle of 30°. Variant 3 had perpendicular baffles, and the last third of the baffle pointed vertically downward and then turned by an angle of 65° from the vertical plane. In variant 4, they were parallel to each other and rotated at an angle of ±20° from the horizontal plane. The used material in the complete model was glass on the front wall of the model and then fireclay panels on the other walls. The air supply was realized by using three types of inlets. The supply of primary air was realized through the retort burner for primary combustion. The second type of air supply was provided through holes on the rear wall of the model for secondary combustion. The supply of the last tertiary air was used to blow the glass. The used material in the complete model was glass on the front wall of the model and then fireclay panels on the other walls. The air supply was realized by using three types of inlets. The supply of primary air was realized through the retort burner for primary combustion. The second type of air supply was provided through holes on the rear wall of the model for secondary combustion. The supply of the last tertiary air was used to blow the glass.
The created mesh of the complete model had nearly 575,500 number of elements, whose shapes were mostly tetragonal, then hexagonal. The quality of the mesh was assessed based on the orthogonal quality with the average value of 0.85, the skewness with the average value of 0.16, and the aspect ratio with the average value of 1.55. Based on the evaluation criteria, the quality of the mesh was recognized as sufficient. The mesh of the complete model is shown in Figure 3. Wall boundary conditions were set up to the baffles, deflector, and retort burner. The burner zone was also used. The exit from the chimney was called the outlet. Inlets were realized by using a retort burner, holes for the secondary inlet as well as a hole for the tertiary inlet.
The created mesh of the complete model had nearly 575,500 number of elements, whose shapes were mostly tetragonal, then hexagonal. The quality of the mesh was assessed based on the orthogonal quality with the average value of 0.85, the skewness with the average value of 0.16, and the aspect ratio with the average value of 1.55. Based on the evaluation criteria, the quality of the mesh was recognized as sufficient. The mesh of the complete model is shown in Figure 3. Wall boundary conditions were set up to the baffles, deflector, and retort burner. The burner zone was also used. The exit from the chimney was called the outlet. Inlets were realized by using a retort burner, holes for the secondary inlet as well as a hole for the tertiary inlet. The grid sensitivity test was realized only on the model of flue gas tract with baffles (variant 4). The denser mesh was created, where the number of elements was 6.7 times higher compared to the base model. The evaluated results showed differences of 2.4% based on the particle number and 2.3% based on the mass flow in the particle outlet stream. As differences in the results were insignificant, less dense mesh was used in all analyzed cases. The final calculation mesh was constructed based on similar element sizing.
The CFD model of the wood stove was developed in Ansys software, where the CFD solver was Ansys Fluent.

The Model Settings
Boundary conditions were set up to mass flow on the inlet and pressure outlet on the stove outlet. Mass flow was calculated by using stoichiometric calculations for maximum heat power of the stove, which was 15 kW. The excess air coefficient was included with a value of 1.6. Mass flow in the primary inlet included products of devolatilization. The ratio of the air in individual inlets was different with the highest mass flow in the primary inlet. The pressure at the outlet was set at −12 Pa and a suitable value for the chimney draft was from 10 to 14 Pa according to STN EN 13,240. For the side and back walls, the convective thermal boundary condition was used in which heat flux was calculated based on the heat transfer coefficient and the temperature on the outer side of the stove. The heat transfer coefficient was calculated based on material properties. For the wall representing glass in the front of the stove, a mixed thermal boundary condition was used, which was similar to the convective boundary condition but included the effect of thermal radiation. Boundary conditions are listed in Table 1. The grid sensitivity test was realized only on the model of flue gas tract with baffles (variant 4). The denser mesh was created, where the number of elements was 6.7 times higher compared to the base model. The evaluated results showed differences of 2.4% based on the particle number and 2.3% based on the mass flow in the particle outlet stream. As differences in the results were insignificant, less dense mesh was used in all analyzed cases. The final calculation mesh was constructed based on similar element sizing.
The CFD model of the wood stove was developed in Ansys software, where the CFD solver was Ansys Fluent.

The Model Settings
Boundary conditions were set up to mass flow on the inlet and pressure outlet on the stove outlet. Mass flow was calculated by using stoichiometric calculations for maximum heat power of the stove, which was 15 kW. The excess air coefficient was included with a value of 1.6. Mass flow in the primary inlet included products of devolatilization. The ratio of the air in individual inlets was different with the highest mass flow in the primary inlet. The pressure at the outlet was set at −12 Pa and a suitable value for the chimney draft was from 10 to 14 Pa according to STN EN 13,240. For the side and back walls, the convective thermal boundary condition was used in which heat flux was calculated based on the heat transfer coefficient and the temperature on the outer side of the stove. The heat transfer coefficient was calculated based on material properties. For the wall representing glass in the front of the stove, a mixed thermal boundary condition was used, which was similar to the convective boundary condition but included the effect of thermal radiation. Boundary conditions are listed in Table 1.
The fuel used in the simulated stove is wood pellets with the following characteristics: (I) carbon content with the value 49.8%, (II) hydrogen content with the value 6.4%, (III) nitrogen content with the value 0.3%, (IV) moisture content with the value 10%, and (V) ash content with the value 0.8%. The calorific value of wood pellets was 15.4 MJ·kg −1 [14].
In these equations were used the following quantities: the fluid density ρ; the velocity vector u; the species index k; mass fraction of k-species Y k ; the diffusion flux of k-species J k ; the pressure p; the stress tensor τ; the gravitational acceleration vector g; the enthalpy h; the heat flux vector q; the number of species in the system m; the source due to exchange of mass (between the continuous phase and the dispersed phase) and additional mass sources S mas ; the source due to exchange of mass of species S k ; the source due to exchange of momentum S mom ; the source due to exchange of energy and additional energy sources S en ; the source due to homogeneous and heterogeneous reactions in which k-species are produced or consumed R k ; the source as the contribution to the energy equation due to the radiation S rad and the source as the amount of energy released in chemical reactions S re [15].
The two-phase model was used in order to simulate the motion of ash particles in the continuous gas phase. Navier-Stokes equations were used for airflow. The motion of particulate matter was modeled by using Lagrangian particle tracking. The trajectories of particles were defined by integrating the force balance on the particle. The force balance equates the particle inertia with the forces acting on the particle, as an Equation (5) where F x is an additional acceleration; F D is the drag force; u is the fluid velocity; ρ is the fluid density; u p is the particle velocity; and ρ p is the particle density [16]. The mass distribution function can be assumed to be normally distributed and calculated by the Rosin-Rammler function [16]  where Y d is the volume fraction of particles with a diameter larger than d. d means the mean parameter and n means the spread parameter. To solve for these parameters, it is necessary to fit particle size data from made size distribution to the Rosin-Rammler exponential equation. The models of four variants of flue gas tract were simulated by a model k-ε realizable with a scalable wall function without the inclusion of combustion processes.
The complete model was simulated first by the model k-ε realizable with a scalable wall function and then also by using the model k-ω SST. The energy control equation was activated and radiative heat transfer was included. The radiative transfer equation for modeling radiation intensity was integrated over solid angles using the discrete ordinates method.
Due to the wide range of particle diameters, it is necessary to use a size distribution. The size distribution results shown in Table 2 were experimentally determined by using a vibrating sieve shaker. The used sieves had a mesh size from 500 µm to 20 µm.

Combustion Process
The combustion of biomass can be divided into three stages: moisture evaporation, devolatilization, and char combustion. Moisture evaporation was not included in the model.
In real stoves, the devolatilization process occurs in the fixed bed in the retort burner. Products of devolatilization process were included in the model as part of mass flow in the burner inlet boundary condition. The amount of released gases was calculated based on fuel analysis and assumed power of the stove.
Char combustion process was not considered in the model. To include heat and mass sources resulting from char combustion, the amount of main reaction products was calculated based on stoichiometry. Then, the combustion products were added into the burner as a volumetric mass source, and consumed oxygen was included as a negative mass source. Additionally, the heat released in the heterogeneous reaction was calculated and included in the model as a volumetric heat source in the burner. Such a simplification could be justified because the detailed analysis of heterogeneous reaction products was not a subject of this paper and because char combustion occurs in a fixed bed, the mass and heat source location is known and is independent of time for constant power of the stove. Values of volumetric sources are listed in Table 3. The combustion of volatiles is realized by two homogenous reactions described in Equations (7) and (8). As the exact composition as well as combustion reaction kinetics of volatiles is not known, it was assumed that volatiles were composed of methane and the amount of released gas was calculated to obtain an enthalpy source equal to the enthalpy of volatiles resulting from fuel analysis.
Appl. Sci. 2021, 11, 1278 Turbulence-chemistry interaction was realized by using the Eddy Dissipation Concept (EDC) model. This model is an extension of the Eddy Dissipation model to involve all the necessary chemical relationships and bonds in the turbulent flow. In EDC, it is assumed that the reaction occurs in so-called fine structures and the reaction rate depends on mass transfer between the fine structures and surrounding gas. Reactions in Ansys Fluent compiled according to the Arrhenius Equation (9) were mutually interconnected using the numerical algorithm ISAT.
where A is the pre-exponential factor; β is the temperature exponent; E is the activation energy for the reaction; and R is the universal gas constant [16]. The kinetic parameters and temperature exponent for the homogenous reaction (7) are as follows: (I) the value of pre-exponential factor A is 1.398762 × 10 10 ; (II) the value of activation energy E is 1.16712 × 10 8 J. kmol −1 ; and (III) the value of temperature exponent β is −0.062. The kinetic parameters and temperature exponent for the homogenous reaction (8) are as follows: (I) the value of pre-exponential factor A is 7.381123 × 10 11 ; (II) the value of activation energy E is 7.65969 × 10 7 J. kmol −1 ; and (III) the value of temperature exponent β is 0.215 [17].

The Flue Gas Tract
First, four variants of baffles geometry were analyzed. The impact of different baffles placing is shown in Table 4. Based on the results, the most suitable model was variant 4 with parallel baffles and rotated at an angle of ±20 • from the horizontal plane. This model was also used in the complete model of the wood stove. The results confirm that the smallest amount of escaping particles was by using variant 4. Despite this, a lot of particles escaped in the flue gas through the chimney.

The Complete Model
Based on the simulation results of different flue gas tract concepts, the geometry of the complete stove was developed. The amount of captured particles is shown in Table 5 and the results indicate that the baffles can mainly trap larger particles. The smaller particles get into the atmosphere transported by the flue gas due to its lower weight in comparison to larger particles, for which the removal effectiveness was at a good level. Particle distribution in the flue gas tract is shown in Figure 4. It can be observed that particles with a diameter close to 150 µm remained in the combustion chamber and did not even get in the flue gas tract. The particles with a diameter greater than 100 µm were captured by baffles located in the flue gas tract and they did not get into the chimney of the model. The particles with a diameter lower than 100 µm could be seen in the chimney area, which means that they escaped into the environment. It can be summarized that the baffles mainly captured particles with a diameter greater than 100 µm, and smaller particles moved out of the chimney to the environment carried by the flue gas.

Used Model Captured Particles [%]
-Depending on the Particle Number Depending on the Mass Flow k-ε [19] 50.95 22.58 k-ω 60.91 31.00 Particle distribution in the flue gas tract is shown in Figure 4. It can be observed that particles with a diameter close to 150 µm remained in the combustion chamber and did not even get in the flue gas tract. The particles with a diameter greater than 100 µm were captured by baffles located in the flue gas tract and they did not get into the chimney of the model. The particles with a diameter lower than 100 µm could be seen in the chimney area, which means that they escaped into the environment. It can be summarized that the baffles mainly captured particles with a diameter greater than 100 µm, and smaller particles moved out of the chimney to the environment carried by the flue gas.  Figure 5 shows a cumulative diameter distribution in the outlet compared with a diameter distribution in the inlet. Based on the data analysis, it can be confirmed that particles with a diameter greater than 100 µm were completely removed in the flue gas tract and the removal efficiency decreased for particles with lower diameter.
The other observed parameter was temperature. The temperature profile is shown in Figure 6 with the size range from 0 to 1500 K. In each case, the highest temperature was between the retort burner and the deflector. The k-ω model achieved only about 1% higher temperatures than in the k-ε model. A different temperature distribution could be observed mainly in the flame region. A higher temperature level resulted from the different combustion rates of volatiles as the turbulence-chemistry interaction rate of the homogenous reaction strongly depends on the calculated turbulence intensity.   The other observed parameter was temperature. The temperature profile is shown in Figure 6 with the size range from 0 to 1500 K. In each case, the highest temperature was between the retort burner and the deflector. The k-ω model achieved only about 1% higher temperatures than in the k-ε model. A different temperature distribution could be observed mainly in the flame region. A higher temperature level resulted from the different combustion rates of volatiles as the turbulence-chemistry interaction rate of the homogenous reaction strongly depends on the calculated turbulence intensity.  The last observed parameter was velocity. The velocity profile in Figure 7 can be seen with the size range from 0 to 3 m•s −1 . In each calculation, the higher speed could be observed in the space above the burner, around the lower part of the baffles, and finally in the flue. Higher values of velocity were achieved in the k-ω model. The maximum achieved velocity in the k-ω model was approximately 1.5 times higher than in the k-ε model. A higher velocity value above the burner and in the flue gas tract was the consequence of a higher reaction rate and higher temperature of the flue gas in the second case. The last observed parameter was velocity. The velocity profile in Figure 7 can be seen with the size range from 0 to 3 m·s −1 . In each calculation, the higher speed could be observed in the space above the burner, around the lower part of the baffles, and finally in the flue. Higher values of velocity were achieved in the k-ω model. The maximum achieved velocity in the k-ω model was approximately 1.5 times higher than in the k-ε model. A higher velocity value above the burner and in the flue gas tract was the consequence of a higher reaction rate and higher temperature of the flue gas in the second case.

Conclusions
Biomass combustion is associated with the formation of emissions such as particulate matter in the atmosphere. Therefore, the present work dealt with this issue. The aim was to optimize the flue gas tract and propose the modification of construction. The best variant was used for the complete model, which was simulated by turbulent-viscosity models. Results showed that application of baffles had a significant effect on the reduction of particles with a diameter greater than 100 µm. Fine particles (diameter less than 100 µm) flowed further through the chimney into the atmosphere. It can be concluded that these baffles could be useful as a pre-separator for larger particles and, then an electrostatic precipitator located in the chimney may be used to separate smaller particles.
In the numerical calculations, the k-ε realizable model with a scalable wall function and k-ω SST model were applied. The standard k-ε model was suitable for a gross estimation of the flow field and models with burning, multiphase flow, and chemical reactions of flow. Improved and modified of this model is a realizable model k-ε. The model

Conclusions
Biomass combustion is associated with the formation of emissions such as particulate matter in the atmosphere. Therefore, the present work dealt with this issue. The aim was to optimize the flue gas tract and propose the modification of construction. The best variant was used for the complete model, which was simulated by turbulent-viscosity models. Results showed that application of baffles had a significant effect on the reduction of particles with a diameter greater than 100 µm. Fine particles (diameter less than 100 µm) flowed further through the chimney into the atmosphere. It can be concluded that these baffles could be useful as a pre-separator for larger particles and, then an electrostatic precipitator located in the chimney may be used to separate smaller particles.
In the numerical calculations, the k-ε realizable model with a scalable wall function and k-ω SST model were applied. The standard k-ε model was suitable for a gross estimation of the flow field and models with burning, multiphase flow, and chemical reactions of flow. Improved and modified of this model is a realizable model k-ε. The model k-ω is more suitable for calculation in places of the boundary layer, where is an adverse pressure gradient The model k-ω SST combines the advantages of both previous models [20,21].
The used turbulent-viscosity models provide the following results: (I) temperature differences were only about 1%; (II) the velocities in some places were up to 1.5 times higher in the k-ω models than in the k-ε models; (III) the differences in the captured particles based on their number were 10% and based on the mass flow were 8.5%. Despite these differences, the models had a similar flow character.
Author Contributions: N.Č.K., investigation, conceptualization, writing-original draft; S.S., data curation; J.J. and A.Č., formal analysis; R.N., writing-review & editing, data curation. All authors have read and agreed to the published version of the manuscript.
Funding: This publication was produced with the support of the Integrated Infrastructure Operational Program for the project: Creation of a Digital Biobank to support the systemic public research infrastructure, ITMS: 313011AFG4, co-financed by the European Regional Development Fund and VEGA No. 1/0479/19: Influence of combustion conditions on production of solid pollutants in small heat sources.