Large-Eddy Simulation of two secondary air supply strategies in kraft recovery boilers

Kraft recovery boilers are large scale combustion applications operating on black liquor, a common side-product of the pulp industry. Here, a simplified boiler model utilizing Computational Fluid Dynamics (CFD) with Large-Eddy Simulation (LES) and Lagrangian Particle Tracking (LPT) is explored to better understand the secondary air supply system and the dispersion of sprayed droplets. The unsteady nature of the air jets and droplet dispersion in such context advocates the usage of scale-resolving simulations such as LES. In the present exploratory study, the usage of LES in recovery boiler air jet simulations is piloted for the first time. The air supply system is modeled as high-momentum-flux jets injected to a uniform cross-flow. First, the set-up was verified by performing a mesh sensitivity study. The main observed global flow features included the mixing zones, wall jet and jet impact regions, jet bending in cross-flow, and reverse flow downstream of the jets. Two different engineering relevant air supply systems, a staggered and an in-line configuration, were studied and compared in terms of mixing and droplet dispersion. The main findings of the present study are as follows. First, out of the two studied configurations, the staggered one was observed to provide a more uniform downstream temperature field. Second, the in-line configuration was noted to outperform the staggered configuration in droplet dispersion by having 22% less spray-wall impingement and 15% more droplets landing at the bottom of the furnace. Third, different modes for droplet trajectories were identified based on the global flow structures.


Introduction
The growing concern of climate change and the scarcity of natural resources increase the need for renewable energy sources.In 2017, on a global scale, the share of bioenergy in final energy consumption was 12.4%, providing half of the renewable energy supply [1].In Finland, the share of wood-based fuels was 28% of the final energy consumption in 2019, from which the share of black liquor (BL), a side-product of the kraft pulping process, is the most significant [2].The share of black liquor in the electricity production in Finland was approximately 10%, which is 22% of the electricity produced by renewable sources.In addition, 53% of industrial heat was produced by black liquor [3].
The combustion air for the droplets in a kraft recovery boiler is mostly provided by the primary and secondary air supply, which are designed to provide a desired temperature profile and thorough mixing to prevent fouling and to reduce emissions [4].It would be a significant technological advantage to be able to accurately model the multiscale physics and chemistry processes inside such boilers for higher efficiency and emission control.In this paper, the secondary air supply system and black liquor droplet dispersion in a kraft recovery boiler is studied using high-resolution computational fluid dynamics (CFD) simulations.
Fundamentally, the secondary air supply of a kraft recovery boiler relates to the canonical flow case of a Jet in a Cross-Flow (JICF).The supply consists of air jets injected in the horizontal direction while a vertical cross-flow is induced by the primary air supply.In addition, the flue gases and combustible gases from the char bed at the bottom of the boiler feed the cross-flow.Holdeman et al. [5,6] studied single, multiple and opposed low-momentum-flux jets in cross-flow experimentally and numerically.They proposed that the key parameters for JICF (single or multiple) are the velocity ratio of the jet to the cross-flow,  =   ∕ ∞ , the momentum-flux ratio  = (   2   )∕( ∞  2 ∞ ) and the penetration constant  = (∕) √  where  is the jet spacing and  is the distance to the confining opposite wall.A thorough review on the fluid dynamical aspects of JICF was given by Mahesh [7].More recently, Liu et al. [8] studied fuel-air mixing in JICF with  = 1.5 and  = 4 using LES.They reported more thorough mixing and deeper jet penetration with higher velocity ratio.Shan and Dimotakis [9] studied JICF experimentally with  = 10 and  = 32 within the Reynolds number range of 1000 ≤  ≤ 20000.They reported better mixing with higher Reynolds number.The velocity ratio in the present cases is between 31 ≤  ≤ 38.
https://doi.org/10.1016/j.applthermaleng.2022  A schematic of a kraft recovery boiler and a boiler CFD modeling approach.(a) A schematic of a kraft recovery boiler [27]. 1 (b) a modeled flow field with RANS showing velocity iso-surfaces of the secondary air jets [13].
Bain et al. [10] studied numerically the effect of the momentumflux ratio and jet spacing.They noted that with larger jet spacing the jet penetration was deeper, forcing more of the flow between the jets.Choi et al. [11] studied experimentally the mixing of high momentumflux jets injected to an industrial furnace.They noted that the previous findings by Holdeman n [5] for the optimal design parameters for low-momentum-flux jets were not applicable for higher values of  .Holdeman [12] suggested that with high momentum-flux ratio, opposing staggered jets provide better mixing than in-line jets or one-sided injection.Maakala et al. [13] studied the effect of the momentum-flux ratio, jet spacing and the penetration coefficient on the mixing of highmomentum-flux jets in confined cross-flow in a recovery boiler context.They noted that the mixing and the variation of the vertical velocity component were increased with higher  .They also reported that the previously defined optimum penetration coefficients were not directly applicable to high-momentum-flux jets.
The dispersion of black liquor droplets is studied in this paper using the Lagrangian Particle Tracking (LPT) method.CFD with LPT has been previously used in the industrial boiler context, e.g. to study ash deposition and fouling [14][15][16][17][18]. Gallen et al. [19] introduced a semi-deterministic LPT method for combustion chamber soot prediction using LES with LPT.They reported encouraging results for the usability of their model when compared to experimental data.Ghasemi et al. [20] studied multi-plume sprays injected into cross-flow jets experimentally and numerically using LES and an Eulerian/Lagrangian multiphase approach.Good agreement between the simulations and the experiments was observed.Modlinski and Hardy [21] studied the applicability of a validated CFD model introduced by Modlinski et al. [22] for O 2 and CO measurements in a pulverized coal boiler.They noted good agreement between the numerical and experimental results.Other related studies include, among others, combustion in grate boilers [23,24], and heat transfer prediction in the superheater region [25,26].Previously, some of the co-authors utilized RANS modeling in the recovery boiler heat transfer optimization context [13,25,27].In contrast, in the present study, LES is used for the first time in such large scale, highmomentum application.Fig. 1 shows a schematic of the boiler and the modeling approach. 1Reprinted from [27] with permission from Elsevier.
For unsteady flows, LES is an attractive numerical method.Wegner and Sadiki [28] studied the mixing of a JICF with LES.They reported good accuracy of LES capturing the unsteady flow and mixing features.Ivanova et al. [29] studied the turbulent Schmidt number in a JICF using LES and RANS.They reported superior accuracy of LES compared to RANS in capturing the turbulent mixing process.Tu et al. [30] conducted an experimental and numerical study using LES of opposing planar jets.They stated that a scale-resolving method, such as LES, is needed to capture the unsteady behavior of such a flow.Thysen et al. [31] compared LES and RANS for opposing plane wall jets in a mixing ventilation related application.They reported that RANS fails to predict correctly the mean flow characteristics in the jet-jet impingement zone.Van Hoof et al. [32] assessed the accuracy of different RANS models and an LES model in a cross-ventilated room with strong air jets.They reported underpredicted turbulent kinetic energy with the RANS models, while the LES produced the mean flow field along with the turbulent kinetic energy profiles with good accuracy.The LES was noted to perform better due to the reproduction of the transient unsteady flow features.
In this paper, the secondary air supply and droplet dispersion in a kraft recovery boiler are studied.Two configurations of the secondary air with engineering relevance are investigated.Based on the literature survey, we note that LES has not been previously used in an industrial boiler high-momentum jet context.Based on the studies cited above, LES is particularly suitable for capturing the flow types associated with such boilers, including JICF and opposing jets.Hence, by recognizing this research gap, LES is selected as the simulation method for this paper.In the present industrial context with relatively high Reynolds number, the numerical quality of LES requires careful consideration: in the present work, both sub-grid metrics and a grid sensitivity analysis are assessed to verify the numerical results.The main goals of the paper are to (i) pilot the usage of LES for the first time in such large scale application, (ii) compare the performance of two different air supply configurations in terms of mixing and droplet dispersion, and (iii) link the large scale flow patterns and droplet dispersion with one another.

Governing equations
The fluid in the present study is considered compressible due to the temperature-induced density variation.The conservation of mass, momentum and energy are governed by the equations where  is the density,  is the filtered velocity field,  is the pressure,  is the gravitational constant,    is the sum of molecular and subgridscale (SGS) eddy viscosity, () = 1 2 (∇+(∇)  ) is the strain rate tensor,  is the sensible internal energy,  = || 2 ∕2 is the kinetic energy per unit mass and    is the sum of laminar and SGS thermal diffusivities.
The momentum of the droplet is determined by where is the drag coefficient when   < 1000.For   ≥ 1000 the drag coefficient is   = 0.424.Droplet relaxation time is defined as where   is the droplet diameter and the subscripts  and  stand for particle and gas, respectively.LES is a turbulence modeling method where the large scale turbulence is resolved and the sub-grid-scale turbulence is modeled.Herein, the -equation model is used [33], which is a one-equation model where the SGS turbulent kinetic energy is given by The eddy viscosity is defined as   =    √  where the filter parameter  =  1∕3  depends on the cell volume   .Furthermore, the constant parameters are   = 0.094 and   = 1.048.

Numerical approach
The governing equations are discretized using a second order accurate finite volume method (FVM) with the open source software Open-FOAM [34].The LES and RANS simulations are both performed with the standard buoyantPimpleFoam [35] solver which is a transient solver for solving the momentum and the internal energy [36].The results of the RANS simulations are summarized in Appendix.The Lagrangian particle tracking library in OpenFOAM is added to the solver by the authors to simulate the particle phase.The pressure-velocity coupling relies on the PIMPLE algorithm, which is a combination of the Pressure Implicit with Splitting of Operations (PISO) [37] and the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) methods [38].For the convection term of the momentum equation, linear central difference scheme was used, while flux limiters are used for the convection terms of the temperature and turbulent kinetic energy equations.For all the governing equations, the diffusion terms were discretized using linear central differencing.A second order implicit scheme is used for temporal discretization.A flow chart for the numerical process is presented in Fig. 2. The flow field is calculated one time step at a time in the PIMPLE loop until a desired end time is reached.The solution of the flow field is used to calculate the droplet trajectories for each time step, but the droplets do not affect the flow field.The flow field is first calculated until it reaches a converged solution before the droplets are added to the domain.More details on the numerical process is presented in [36].

Computational domain
A kraft recovery boiler air supply system includes primary, secondary and tertiary air.The black liquor settling on the bottom of the furnace forms a char bed, which affects the primary air velocity distribution and provides flue gases and combustible gases to the system.Here, the char bed, primary air and tertiary air are not modeled.However, the mass flow rate from the primary air and the gases from the char bed are indirectly considered via the uniform cross-flow.The in-line configuration includes 10 air jets, five on each walls (front and rear walls) as seen in Fig. 3(c).The locations are similar to the staggered configuration, however, in addition to the main jets, smaller brake jets are introduced in-line to the main jets.The purpose of the brake jets is to prevent the main jets from impinging at the opposite wall [13].Hence, the jet spacing here is halved to  =  ∕6.The jet openings in both cases are rectangular and the width of the openings is .The total mass flow rate introduced by the jets is the same in each configuration, and hence, the jet openings in the in-line configuration are smaller in height.
The mesh contains 22M hexahedral cells.The mesh is finer in the jet region and coarser in the downstream region.The coarser mesh in the downstream region is justified since the focus of this study is in the jet region and the downstream region has relatively small influence on the results in the jet region.The number of cells per jet width in the jet openings are 24.In the in-line configuration the cells per jet height in the opening varies between 18-24 and in staggered configuration it is 42.The mesh for the in-line configuration can be seen in Fig. 4. The mesh for the staggered configuration is otherwise identical except for the different number of jets.A schematic of the inlet flows is shown in Fig. 4(c) along with the mesh.

Boundary conditions
The cross-flow from the bottom of the domain is identical in both configurations.The mass flow rate is ṁ = 43.705kg∕s at a temperature of  = 1100 • C. As mentioned, the shape of the velocity profile due to the primary air and the char bed is neglected and the flow is modeled with a uniform velocity profile.The outlet condition for the velocity is set such that it prevents reverse flow.Zero gradient is applied when the velocity vector is pointing out of the domain, while a fixed value ( = 0) is utilized for inward pointing flow.A respective condition is applied for the temperature and  as well, where the fixed value for the reverse flow temperature is set to the downstream mean value.At the walls, a no-slip condition is used for the velocity and zero-gradient (adiabatic) for the temperature.The droplets are modeled as solid particles with one way coupling and with no swelling, evaporation or chemistry.While this is a strong assumption, the main focus here is to link the large scale flow patterns with the droplet trajectories.Hence, droplets with constant density and size are used.The droplets are injected in a 45 • cone angle at a velocity of 9 m/s and a −42 • tilt angle from the horizontal level.The droplet diameter varies between 0.236 mm and 10.238 mm as a continuous uniform distribution and the droplets have a constant density of 346.9 kg/m 3 .The density is based on an estimation of the density of the droplets after swelling in the furnace right after the liquor burner.The droplets are modeled to stick on the walls upon impact.The particle injectors are located at  = ℎ above the secondary air level, all on the same horizontal plane.There is one injector at the front and rear walls respectively, and three injectors at the left and right walls, adding up to a total of eight injectors.The injectors at the left and right walls are located so that the middle injector is at the center and the spacing is .The orientation of the injectors is shown in Fig. 3(d).

Results and discussion
In the following, a mesh sensitivity analysis is first shown for the more complex in-line configuration, followed by an analysis of the flow features in both air supply configurations.The configurations are then compared regarding mixing and droplet dispersion.

Grid sensitivity analysis
The grid sensitivity analysis is shown here for the in-line configuration, since it is considered to be fluid dynamically a more complex configuration.The complexity arises from the impingement of the opposing jets as well as the smaller jet openings, resulting in a lower number of cells per jet diameter.Three different meshes are compared with similar refinement regions.In the finest region near the jets, the cell size is  = ∕16, where  is the size of the largest cell.The cells in the second, third and fourth refinement regions are  = ∕8,  = ∕4 and  = ∕2 respectively.The three different meshes are generated by varying the largest cell size in the background grid.The resulting total cell counts are 6M (coarse), 11M (medium), and 22M (fine) cells.Fig. 5 shows the temperature field, the total (resolved + modeled) turbulent kinetic energy, and the ratio of the SGS turbulent kinetic energy to the total for the three different meshes.The figure shows the 6M, 11M and 22M cell cases from left to right.Figs.5(a) and 5(b) represent cutplanes at the level of jet center while Fig. 5(c) is considered one jet opening height above the jet centers, with relevance to the temperature mixing process above the secondary air jets.In the temperature field, the jet shear layer Kelvin-Helmholtz instability is the clearest at the finest grid.On the coarsest grid, the Kelvin-Helmholtz vortices are not observable before the impingement zone.The jet-level turbulence magnitude between the 11M and 22M cases is similar, however, the 6M case differs significantly from the two other cases.Above the jet level, the percentage of the modeled  to the total  =   ∕(  +   ) is less than 10% in the 22M case, however, even with the coarsest grid the percentage is less than 20%, which has been considered sufficient [39].The results between the 11M and the 22M cell cases do not differ significantly, indicating that either one of the meshes would sufficiently well resolve the turbulent flow field.However, the 22M cell mesh is chosen as the main mesh and from here on, all the results are shown for that mesh only.Based on the grid sensitivity study, we conclude that the 22M mesh can be used for realistic and physical results.
A comparison between LES and RANS is also carried out as provided in Appendix.The mean flow features are noted to be qualitatively similar with certain differences noted in jet bending, jet core length, shear layer growth rate, and jet-jet impingement locations.However, the most prominent differences can be noted in the unsteady flow features captured by LES, such as the shear layer unsteadiness and jetjet impact.The results of the comparison advocate the usage of LES for the main simulations herein.

Flow topology
In this section, the global coherent flow features are identified.The analysis is performed for both air supply configurations separately.
Staggered configuration.Fig. 6(a) illustrates the iso-surface of the velocity magnitude (|| = 20 m∕s) of the five jets in the staggered configuration along with a vertical cut-plane visualization of the vertical velocity component.2D streamlines from the time-averaged flow field are drawn on the vertical cut-plane to further illustrate the coherent flow structures.The jets (i) penetrate readily and impinge on the opposite wall.The jets are noted to bend slightly due to the crossflow.The impingement of the jets combined with the upward bending creates wall jets (ii).The wall jets create coherent upward currents that induce smaller vortices (iii) and even reverse flow (iv) in the center of Fig. 7 shows the instantaneous and time-averaged velocity fields on a horizontal cut-plane at the jet level.2D streamlines are drawn on the time-averaged field.In the instantaneous field, it is noted that the jet tip almost reaches the opposite wall (i).It can be observed from the streamlines, that due to the staggered alignment of the jets, the antiparallel shear layers of the jets meet (ii), creating coherent swirling motion between the jets.The corner vortices (iii) noted in Fig. 6 can be observed here as well.
In-line configuration.Fig. 8(a) illustrates the iso-surface of the velocity magnitude (|| = 20 m∕s) of the opposing jets in the in-line configuration along with a vertical cut-plane visualization of the vertical velocity component.2D streamlines are drawn on the cut-plane to illustrate the flow structures.The jets (i) penetrate until the impingement zone.The bending of the jets is scarcely visible due to the shorter penetration distance.However, the minor bending of the jets can be observed to induce a small region of strong upward current (ii) in the impact region.Here, the momentum-flux ratio is between 3000 <  < 4500, and hence, the jets are not expected to experience significant bending.Qualitatively, the velocity field is more uniform in the downstream region, however, more coherent upward flow (iii) is observed in the center of the domain.An illustration of the 3D streamlines can be seen in Fig. 8(b).Vortical structures (iv) induced by the opposing jet impingement are observed in the mixing zone.The more coherent upward current in the center of the domain (v) is observed while the upward flow near the walls (vi) is noted to be less significant.The jets in the in-line configuration are noted to bend less with the jet core remaining intact for longer compared to the staggered configuration, consistent with previous studies [10].
The instantaneous and time-averaged velocity fields are shown in Fig. 9.The 2D streamlines are drawn on the time-averaged field.In the instantaneous field, stabilized regions (i) are noted between the jet openings.The stabilization is primarily explained by the parallel motion of the jets which reduces the relative shear and hence avoids merging of the opposite shear layers.The self-stabilization of the jets is noted to increase the transition length.The rapid mixing motion induced by the jets is noted to extend to only a small distance from the impingement zone (ii).The center point of the impact is noted to shift slightly from the average location (iii).The previously noted coherent vortical structures in the impingement zone (iv) are clearly observable from the streamlines here as well.

Mixing of cool and hot gases
In kraft recovery boilers the secondary air acts as an oxidizer for combustion.A common objective is to prepare a well-mixed thermal environment for combustion.In the following, the mixing characteristics in the two configurations are compared.The instantaneous temperature fields are marked in Fig. 10 as a vertical cut-plane.The mean temperature is defined as the downstream bulk temperature, and is calculated as   = ∫    ∫   ≃ 1030 K.It is the same for both configurations due to the similar inlet mass flow rates and temperatures.The hot gases flow  upwards from the bottom of the boiler (i) and mix with the cool gases supplied from the secondary air jets (ii).The most intense mixing is noted in the mixing zone (iii), induced by the shear layers of the jets.It is noted that the mixing zone in the staggered configuration is larger compared to the in-line configuration.In the staggered configuration, the wall jets create regions with lower temperatures (iv), however, the temperature variation is qualitatively not large.In the in-line configuration, large pockets of non-mixed hot gases (v) are noted downstream of the jets.The hot gases flow upward from the self-stabilized region between the jets outside of the mixing zone.The upward flow due to the jet impact creates a large pocket of non-mixed cool gas (vi).
The probability density functions (PDFs) for the temperature at the jet level and at the droplet injection level are illustrated in Fig. 11.The calculated downstream bulk temperature (  = 1030 K) is illustrated as a dashed line in the figures.As seen in Fig. 11(a), in the staggered configuration, at the secondary air level, the temperature variance is higher with certain level of mixing.The cool jets are observed as a temperature peak at 440 K.At the injection level, the temperature variation is relatively small, indicating thorough downstream mixing.It is also noted that the temperature peak in the PDF is positioned at the average downstream bulk temperature value.
As noted in Fig. 11(b), in the in-line configuration, the temperature variance at the secondary air level as well as at injection level is considerably higher compared to the staggered configuration.At the secondary air level, a peak at the cross-flow temperature (1373 K) can be noted.This indicates that the hot cross-flow passes by the jets with minimal level of interaction due to the smaller area of the mixing zone, as previously noted in Fig. 10(b).At the droplet injection level, the temperature PDF is observed to be bi-modal, indicating the presence of pockets with hot and cool gases.This is further addressed by noting the bulk temperature between the two modes with a local minimum value.Based on the present numerical observations, the level of mixing is higher with the staggered compared to the in-line configuration.A more thorough mixing with the staggered configuration  compared to the in-line jets is also consistent with the findings of Holdeman [12], who studied the effect of the penetration constant  in high-momentum-flux jet systems.

Droplet dispersion
Staggered configuration.The residence time of the droplets as a function of the total distance traveled is shown in Fig. 12.For clarity, the droplets are only shown for a single injector (right injector 2 in Fig. 3(d)).The droplet injection locations are illustrated in the figure by the vertical cylinders at the walls (i).The residence time is defined as the time the droplet has spent in the domain from the injection time until the time when it attaches to a boundary.The droplets are modeled here to attach the wall upon impact.In addition, the trajectories of several droplet modes are marked in the figure.The different modes are identified based on the global coherent flow structures as well as the droplet diameter.
Mode I.These are the largest droplets ( > 8 mm) that experience little to no jet influence, depositing exclusively on the bottom.Some interaction with the jets can be noted (ii), however the droplet trajectories seem to retain the initial conical injection pattern.These heavy droplets show an almost linear correlation between the traveled distance and the residence time, indicating that the velocity and the path of the droplets are not significantly affected by the jets.
Mode II.This mode includes the medium-sized droplets (4 mm <  < 8 mm) that interact more strongly with the jets.The droplets are noted to be pushed back and forth by the jets (iii) at the jet level, however, most of the droplets are still observed to deposit at the bottom of the furnace.

Mode III.
Here we show the medium-small droplets (2 mm <  < 4 mm) that are pushed by the jets to the opposite walls above the jet level (iv).These droplets are carried by the shear layer of the jets  until they deposit at the walls.The residence time compared to the traveled distance is relatively scattered, indicating that the droplet-jet interaction is much less predictable compared to the larger droplets.
Mode IV.This mode is identified based on the large number of the droplets shown in the residence time scatter plot in a small area.These droplets are all of smaller size ( < 2 mm) and they are noted to be carried away by the coherent upward current induced by the wall jets (v) shown in Fig. 6.The larger deposit of these droplets in the residence time scatter plot indicates a predictable flow pattern.
In-line configuration.The droplets residence time as a function of the total distance traveled from a single injector can be seen in Fig. 13 Mode I. Here, similar to the staggered configuration, the heavier droplets ( > 8 mm) pass the jets with little or no interaction and retain the conical injection shape.Some minor interaction with the jets is noticed (i), however, the almost linear pattern of the residence time compared to the traveled distance reveals predictable flow patterns.Mode II.Similar back and forth motion of the medium (4 mm <  < 8 mm) droplets (ii) can be observed compared to the staggered configuration.The main difference here is that due to the opposing jets, the droplets keep going back and forth (iii) and are less likely to reach the walls.Nevertheless, some wall deposition is noted (iv).
Mode III.The medium-small (2 mm <  < 4 mm) droplets are noted to be affected by the upward motion of the jets due to the jet impact (observed in Fig. 8).These droplets bounce upwards from the jet impact (v), however, due to the weaker upward currents near the walls, the droplets travel back downwards due to gravity (vi).Although possible, it is unlikely that such a bouncing droplet would encounter a strong upward current multiple times at the jet level and hence those droplets will fall eventually to the bottom of the boiler.
Mode IV.This group of the smallest droplets ( < 2 mm) shows the most unpredictable behavior since they are not observed to follow any coherent flow structure (see Fig. 8).Instead, the droplets seem to drift slowly upwards due to the upward bulk flow, seemingly randomly (vii).

Droplet deposition
In real kraft recovery boilers, ideally, the droplets land on the bottom of the furnace.Fouling of the furnace walls is not desired since the black liquor layer building up at the walls might fall down to the char bed as large chunks, causing possible problems in the boiler combustion process.Additionally, the droplets carried away by the flow field, if not incinerated in the air, end up fouling the superheaters at the top part of the boiler.Fouling of the superheaters lowers their heat transfer performance and should also be avoided.
Staggered configuration.The particle deposition in the staggered configuration is presented in Fig. 14.The deposition is visualized in Fig. 14(a).The figure shows that while most of the droplets land in the bottom of the domain (i), a significant number is pushed by the jets to the front and rear walls in the jet-wall impingement zone (ii).The deposition pattern of the droplets at the front and rear walls correlates with the configuration of the jets which transport the spray particles.Droplets can also be observed at the top of the domain on the rear side (iii).The droplet deposition is noted to correlate with the large scale flow structures near the boundaries shown in Fig. 6.
The relative share of droplets on each wall from each injector is illustrated in Fig. 14(b).Here, droplets depositing above  = ∕4 are considered to have been carried away by the flow field.The numbering of the left and right injectors is shown in Fig. 3(d).It is noted that the most significant amounts of droplets are found at the bottom, while a considerable amount deposits at the front and rear walls.The number of droplets depositing at the left and right walls is noted to be insignificant, which is related to the lack of strong horizontal flows normal to the left or right walls.
Fig. 14(c) shows the z-coordinate of the droplets' location as a function of the particle diameter.It can be noted that only smaller droplets get carried away by the flow, while the largest droplets are found below the jet level exclusively.The largest concentration of wall deposited droplets can be observed between the particle injection and the secondary air levels, at the jet-wall impingement zone (see Fig. 6).
In-line configuration.Fig. 15(a) visualizes the droplet deposition in the in-line configuration.Compared to the staggered configuration, a much less significant number of droplets can be observed at the front and rear walls.A more significant amount of deposition is found at the bottom (i).However, droplet wall deposition peaks at the jet-wall impingement zones (ii).A small amount of droplets is noted to deposit at the top boundary (iii), correlating with the observations of the flow characteristics made in Fig. 8.
The relative amount of the droplets on each wall is shown in Fig. 15(b).A larger portion of the particles settle at the bottom compared to the staggered configuration, and a significantly smaller number of droplets are deposited on the rear and front walls.Similar to the staggered configuration, the relative amount of droplets depositing at the left and right walls is small.The total relative amount of droplets depositing at the walls for both configurations is quantified in Table 1.Compared to the staggered configuration, 15% more droplets landed on the bottom of the boiler and a 22% decrease in wall deposition is observed.A small decrease is noted in the particles carried away by the flow field.
The vertical location of droplet deposition as a function of the particle diameter is depicted in Fig. 15(c).It can be noted that in the in-line configuration the diameter at which droplets are carried away is considerably smaller than in the staggered configuration.This is expected since the vertical velocity component near the injectors is weaker reducing droplet transport upwards.Larger droplets are not found above the injection level.It is noted that the addition of the brake jets reduces the velocity magnitude of the flow interacting with the walls.This reduction is directly linked to the extent of droplet deposition at the rear and front walls.Hence, in terms of the droplet dispersion, the in-line configuration turns out to outperform the staggered one.

Future work
The present work focused on jet-jet interaction and mixing investigations using LES in a simplified numerical setup.As a limitation of the present work, we note that additional case configurations with different number of jets would be needed to draw more general conclusions between staggered and in-line configurations.The effects of the primary air and the shape of the char bed at the bottom of the boiler could still be further considered.Those effects could have a major effect on the formation of large scale flow structures.In the present configuration the jets may bend downwards while in reality the char bed would act as a physical barrier preventing such motion.The effect of the jet spacing to the shear layer and transition length could be studied further.The swelling, evaporation and combustion of the black liquor droplets could be considered as well.The heat release from the combustion is expected to have a further effect on the buoyancy-driven flow in the boiler.In  addition, the geometrical shape of the boiler and the heat exchangers could be taken into account in future studies.Last, modeling the wall roughness due to fouling is presently not accounted for and wall-models for such porous zones could be further developed.

Conclusions
The secondary air supply and droplet dispersion in a kraft recovery boiler with two different air supply configurations (staggered   cross-flow, and the black liquor droplets were simulated as simplified Lagrangian particles with one way coupling.In terms of coherent fluid mechanical phenomena, a mixing zone, jet impingement, wall jets, jetjet impact, jet bending and reverse flow were observed.The two air supply configurations were compared in terms of mixing and droplet dispersion.The core physical difference, the staggered configuration poses merging opposite shear layers, while the in-line configuration involves jet-jet impingement.The main differences in the simulation results are explained by these physical differences.From the two studied configurations, the level of mixing was noted to be stronger with the staggered configuration from the viewpoint of temperature uniformity.However, the in-line configuration was observed to perform slightly better than the staggered configuration based on the droplet dispersion.For the in-line configuration, a total of 15% more droplets landed in the bottom of the boiler as desired while 22% less wall deposition was noted.A relation between the global flow characteristics and the droplet dispersion was observed from droplet trajectories.The larger droplets were not essentially influenced by the flow field while the dispersion of smallest droplets was noted to follow the observed flow structures closely.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix. Comparison to RANS
Here, the results obtained with LES are compared to RANS.For the RANS simulations,  −  SST model [40] is chosen. −  SST model was used previously by the authors [13,25,27] in recovery boiler context.Both secondary air jet configurations are compared separately.We remind the reader, that the LES results are verified with a grid convergence study in addition to sub-grid metrics.
Staggered configuration.Fig. A. 16 shows a vertical cut-plane of the (time-averaged) velocity fields in addition to the streamlines of LES and RANS in the staggered configuration.The velocity fields are qualitatively similar, however, differences can be seen in the shear layer and the jet core.The jet core is noted to be slightly longer in the RANS.
The total turbulent kinetic energy field in a vertical cut-plane for the staggered configuration can be seen in Fig. A.17.The character and magnitude in the considered metrics are mostly similar, however, the RANS model is observed to predict a lower level of turbulence in the shear layer.Clear differences are noted between the two methods in the turbulence transition region.This can be explained by the growth of the Kelvin-Helmholtz instability which is captured by the LES.
The velocity fields and the total turbulent kinetic energy fields in a horizontal cut-plane at the jet level can be seen in   and RANS fields are noted to be mostly similar with minor differences.Differences are noted in the jet core length, jet shear layers and the level of turbulence.
In-line configuration.impingement point of the jet-jet impact between the two approaches.Such an opposing jet flow is inherently unsteady, which is expected to be challenging to capture using RANS.The results between the LES and the RANS reveal differences mostly in the unsteady features of the flow, such as the shear layer growth and the jet-jet impact.These features are of importance in terms of overall flow pattern formation, mixing, and droplet dispersion.The comparison advocates the usage of LES for the simulations presented in this paper.

Fig. 1 .
Fig. 1.A schematic of a kraft recovery boiler and a boiler CFD modeling approach.(a) A schematic of a kraft recovery boiler [27]. 1 (b) a modeled flow field with RANS showing velocity iso-surfaces of the secondary air jets [13].

Fig. 2 .
Fig. 2. A process flow chart for the numerical solution.The solution of the flow field affects the solution of the droplets, but the solution of the droplets does not affect the flow field.

Fig. 3
shows a schematic representation of the boiler.The areas that are modeled are indicated in the figure.In Fig. 3(b) the horizontal dimensions are shown.The width of the domain is  while the domain depth is .The width and depth are roughly ten meters, and the two configurations chosen are relevant to real industrial boilers.The total height of the domain is  with the bottom of the domain located at  = −∕12.5and the secondary air level at  = 0 m.The height is chosen arbitrarily and is considered sufficiently large to avoid the outlet boundary influencing the results in the upstream regions.The domain is simplified since the main focus is on the flow structures of the secondary air supply and their effect on the dispersion of the droplets.The two air supply configurations can be seen in Figs.3(b) and 3(c) shown as a cut-plane at the secondary air level at  = 0 m.The black liquor injector configuration is shown in Fig. 3(d) at the injection level at ℎ = 0.15 H. Two different secondary air supply configurations are studied, a staggered and an in-line configuration.The staggered configuration includes five air jets, three on the front wall and two on the rear wall.The locations of these jets are indicated in Fig. 3(b).The jet spacing is  =  ∕3, i.e. the jets are distributed evenly along the walls.

Fig. 3 .
Fig. 3.The main dimensions of the boiler used in the present study.(a) A schematic of the boiler indicating the regions that are modeled.(b) Horizontal dimensions of the domain indicated from a cut-plane at the secondary air jet level at  = 0 m (staggered configuration).(c) The jet arrangement with the in-line configuration.The shorter jets indicate smaller mass flow rate.(d) The configuration of the black liquor injectors at the black liquor injection level at  = ℎ.The numbering of the left and right injectors is depicted at the right wall.
The total mass flow rate from the secondary air supply in both configurations is the same ṁ = 28.64 kg∕s at a temperature of  = 170 • C. In the staggered configuration, each of the five openings are similar, and the opening area, jet velocity, Reynolds number, momentum-flux ratio and the penetration constant are  = 0.1027 m 2 ,  = 70 m∕s,  = 697800,  = 3800 and  = 18, respectively.In the in-line configuration, the parameters for the main jets are  = 0.05625 m 2 ,  = 76.8 m∕s,  = 590400,  = 4570 and  = 9.9.For the brake jets at the front wall the parameters are  = 0.04125 m 2 ,  = 62.7 m∕s,  = 407100,  = 3050 and  = 8.1, and for the brake jets at the rear wall  = 0.04125 m 2 ,  = 74.2m∕s,  = 481800,  = 4270 and  = 9.6.

Fig. 4 .
Fig. 4. The structure of the computational mesh used in this study.(a): Vertical cut-plane showing the main refinement regions.(b): Horizontal cut-plane at the jet level showing the smallest refinement regions at the jet openings.(c): a zoomed horizontal cut-plane showing the refinement regions along with a schematic of the velocity inlets.

Fig. 5 .
Fig. 5. (a): Temperature field at the jet level.(b): sum of the diagonal components of the resolved and SGS Reynolds stresses at the jet level.(c): The ratio of the SGS turbulent kinetic energy to the total turbulent kinetic energy at one jet opening height above the jet level.From left to right: 6M, 11M, 22M cell cases.

Fig. 6 .
Fig. 6.Global flow structures in the staggered configuration.(a): The jets are visualized as isosurfaces of the velocity magnitude at  = 20 m∕s while the vertical cut-plane depict the vertical velocity component.The streamlines are drawn on the 2D cutplane to illustrate coherent flow structures.(i) Jet, (ii) wall jet, (iii) vortices induced by the wall jets, (iv) reverse flow.(b): 3D streamlines illustrating coherent flow structures.(v) corner vortices, (vi) coherent upward currents, (vii) less coherent reverse flow.

Fig. 7 .
Fig. 7.The instantaneous (left figure) and the time-averaged (right figure) flow fields at the jet level in the staggered configuration.The streamlines are drawn on the mean field to illustrate coherent features.(i) Jet tip, (ii) merging of the anti-parallel shear layers, (iii) corner vortices.

Fig. 8 .
Fig. 8. Global flow structures in the in-line configuration.(a): The jets are visualized as isosurfaces of the velocity magnitude at  = 20 m∕s while the vertical cut-plane depicts the vertical velocity component.The streamlines are drawn on the 2D cutplane to illustrate coherent flow structures.(i) Jet, (ii) upward momentum due to the jet impact, (iii) coherent upward flow.(b): 3D streamlines illustrating coherent flow structures.(iv) vortical structures in the mixing zone, (v) coherent upward currents, (vi) less coherent currents near the rear wall.

Fig. 9 .
Fig. 9.The instantaneous (left figure) and the time-averaged (right figure) flow fields at the jet level in the in-line configuration.The streamlines are drawn on the mean field to illustrate coherent features.(i) Stabilized region between the jets, (ii) edge of the mixing zone, (iii) center point of the jet impact, (iv) vortical structure in the impingement zone.

Fig. 10 .
Fig. 10.Vertical cut-planes of the instantaneous temperature fields in the staggered (a) and in-line (b) configurations.(i) Hot gases induced as the cross-flow, (ii) cool gases from the secondary air jets, (iii) rapid mixing in the mixing zone, (iv) cooler gases in the wall jet region, (v) pockets of non-mixed hot gases, (vi) non-mixed cool gases in the upward flow region.

Fig. 11 .
Fig. 11.PDFs of temperature at two horizontal planes at the secondary air level and at the droplet injection level in staggered (a) and in-line (b) configurations.
. The orientation of the figure is similar to Fig. 12.The different modes of the droplet trajectories are identified in the figure.

Fig. 12 .
Fig. 12.The residence time of the droplets in the staggered configuration as a function of the total distance along with the trajectories of several main groups of droplets from a single injector (right injector 2).The main groups are indicated in the scatter plot.The figure shows the main droplet patterns as particle trajectories related to the global flow characteristics.Mode I: larger droplets, mode II: medium-large droplet, mode III: medium-small droplets, mode IV: small droplets.(i) illustration of the droplet injectors, (ii) minor droplet-jet interaction, (iii) strong droplet-jet interaction, (iv) wall deposition, (v) droplets getting carried away.

Fig. 13 .
Fig.13.The residence time of the droplets from all injectors in the in-line configuration as a function of the total distance along with the trajectories of several main groups of droplets from a single injector.The main groups are indicated in the scatter plot.The figure identifies the main droplet patterns.Mode I: heavier droplets, mode II: medium-heavy droplet, mode III: medium-small droplets, mode IV: small droplets.(i) minor droplet-jet interaction, (ii) strong droplet-jet interaction, (iii) droplet bouncing in the jet-impingement zone, (iv) wall deposition, (v) droplets bouncing upward from the jet impact, (vi) bounced droplets traveling back downwards due to gravity, (vii) seemingly random droplet trajectories.

Fig. 14 .
Fig. 14.Droplet deposition with the staggered configuration.(a) A visualization of the droplets with iso-surfaces of the velocity magnitude.(i) Droplets landing at the bottom, (ii) wall deposition, (iii) droplets getting carried away by the flow field.(b) The relative amount of droplet deposition on each boundary from each injector.(c) Vertical location of the deposited droplets as a function of the droplet diameter.

Fig. 15 .
Fig. 15.Droplet dispersion with the in-line configuration.(a) A visualization of the droplets with iso-surfaces of the velocity magnitude.(i) droplets landing at the bottom (ii) wall deposition, (iii) droplets getting carried away by the flow field.(b) The relative amount of droplet deposition on each boundary from each injector.(c) Vertical location of the droplets on the walls as a function of the droplet diameter.The figure shows the droplets final location and the relation to the jets.

Fig. A. 16 .
Fig. A.16.The time-averaged LES velocity field and streamlines compared to RANS in a vertical cut-plane for the staggered configuration.

Fig. A. 17 .
Fig. A.17.The total turbulent kinetic energy field in a vertical cut-plane for the staggered configuration.

Fig. A. 18 .
Fig. A.18.(a) Time-averaged velocity field and (b) total turbulent kinetic energy field in a horizontal cut-plane at the jet level for the staggered configuration.
Fig. A.18.The LES

Fig. A. 19 .
Fig. A.19.The time-averaged LES velocity field and streamlines compared to RANS in a vertical cut-plane for the in-line configuration.

Fig. A. 20 .
Fig. A.20.The total turbulent kinetic energy field in a vertical cut-plane for the in-line configuration.
Fig. A.19 shows a vertical cut-plane of the (timeaveraged) velocity fields and streamlines captured by the LES and RANS methods in the in-line configuration.Some differences are noted in the

Fig. A. 21 .
Fig. A.21.(a) Time-averaged velocity field and (b) total turbulent kinetic energy field in a horizontal cut-plane at the jet level for the in-line configuration.
Fig. A.20 shows comparison of the total turbulent kinetic energy fields in a vertical cut-plane.The magnitude of the turbulent kinetic energy is noted to be similar in LES and RANS.The LES shows a clear transition to turbulence.The stagnation point characteristics, including jet bending and size of the high-turbulence zone, are noted to be slightly different in LES compared to RANS.The level of turbulence in the shear layer in the RANS is noted to be more uniform compared to the LES.The time-averaged velocity fields and the total turbulent kinetic energy fields in a horizontal cut-plane at the jet level are compared in Fig. A.21.The figure shows clear differences in the jet-jet impact zone between LES and RANS.In particular, on the considered cutplane, the RANS field shows clear jet-jet impact only for the central jets.Horizontal bending of the outer-most jets is noted in the RANS simulations.

Table 1
Relative amount of droplets depositing on each wall in both configurations.