Incorporation of flamelets generated manifold method in coarse-grained Euler-Lagrange simulations of pulverized coal combustion

(cid:1) Application of FGM and coarse-graining approaches in the framework of CFD-DEM. (cid:1) State-of-the-art pulverized coal combustion models. (cid:1) Concentrated fuel release induced by coarse-graining can be mitigated by smoothing. (cid:1) Simulations are validated with experimental data of a jet ﬂame of pulverized coal.


Introduction
The world is undergoing a transition from a fossil-fuel based to a renewable-energy based economy (IPCC, 2018).Before the complete transition to renewable energies, the dependency on solid fuels remains relevant for at least the next 20 years (IEA, 2021).It is important that the transition can be as clean as possible to limit the rise of global temperature within 1.5°C to pre-industrial levels, as manifested in the 2021 Glasgow agreement (UNFCCC2021,).The efficient use of solid form fossil fuel, i.e. coal, and the technological leap toward the minimization of pollutant emissions are priorities for industries.To achieve these goals, cost-effective numerical simulations have become a reliable and promising tool for understanding the complexities involved in solid fuel combustion, which can be infeasible for direct measurement due to harsh conditions, such as high temperature and high pressure in blast furnaces (Geerdes et al., 2015) or in pressurized coal combustors (Hong et al., 2009).In this work, we utilize the Computational Fluid Dynamics -Discrete Element Method (CFD-DEM) to numerically study pulverized coal combustion (PCC) and we aim to incorporate different costreduction techniques in order to build a basis for solid fuel combustion-related industrial applications.
In the CFD framework, Large Eddy Simulation (LES) has been widely applied to PCC (Hasse et al., 2021;Yamamoto et al., 2011;Stein et al., 2013;Franchetti et al., 2013;Vascellari et al., 2013;Watanabe and Yamamoto, 2015;Wen et al., 2017;Rieth et al., 2017;Akaotsu et al., 2020).Often LES of PCC employs an Eulerian-Lagrange approach coupled to a flamelet-based tabulated chemistry approach such as the Flamelet Generated Manifold (FGM) method (Vascellari et al., 2013;Watanabe and Yamamoto, 2015) to account for detailed kinetics with a reasonable computational cost for reaction modeling.Previous studies have developed a multi-stream FGM approach.Watanabe and Yamamoto (2015) first proposed a two-stream FGM method to account for two types of fuels coming from devolatilization and char oxidation processes.Wen et al. (2017) and Rieth et al. (2017) considered the effect of turbulent combustion by using a presumed probability density function while maintaining a relatively cheap cost.Akaotsu et al. (2020) extended the work to study the effect of pulverized coal supply on the flame structure.
Although many studies addressed the effectiveness of the FGM approach, not much research has investigated the reduction method on Lagrangian objects in the scope of PCC.Often EL simulations of particulate flow are bottlenecked by the discrete phase due to the enormous amount of tracked objects, especially in the applications of industrial processes.As a result, a collision-based Coarse-Graining (CG) method (Bierwisch, 2009) is introduced.The essence of the CG method is to ease the computational burden by replacing a number of small primary particles with a representative large parcel while keeping the energy lost per collision event equivalent.This CG approach has been widely applied for granular flow studies (Radl et al., 2011;Nasato et al., 2015;Queteschiner et al., 2018) and been extended to Eulerian-Lagrangian simulations of the hydrodynamics of particulate flows (Chu et al., 2015;Napolitano et al., 2022;Jiang et al., 2021;Tausendschön et al., 2020).However, its incorporation in EL simulations of reactive particulate flows is so far less reported in the literature.
In this work, we combine both the FGM and CG methods with the CFD-DEM model to reduce the computational burden on both the continuum and the discrete phases for simulations of reactive particulate flows.The model (referred to as CG-EL-FGM in the rest of the study) is validated with a laboratory-scale pulverized coal jet flame experiment performed at the Japanese Central Research Institute of Electric Power Industry (CRIEPI) (Hwang et al., 2007;Hwang et al., 2007).The goal is to demonstrate the feasibility of using much fewer computational resources to achieve a consistent level of performance.The CG-EL-FGM simulation results of particle trajectories, velocities, and species concentration are compared with that from the experiment.Additionally, the reduction in computation time is given to show the strength of the current reduction models.Finally, the effects of particle coarse-graining factor on the hydrodynamics as well as thermochemical behavior are investigated, from which an optimal scaling factor is suggested for future CRIEPI burner studies.

Numerical methods
In this section, the numerical methods of the discrete phase including the coarse-graining adaption, a description of the continuous phase, and the numerical setup are outlined.

Particle motion
The motion of particles are calculated using the Discrete Element Mehotd (DEM) due to its ability to resolve realistic contact forces.The translational and rotational motion of particles is described using Newton's second law of motion.Particles are subjected to gravity, contact forces with other particles or walls, and drag force (Golshan et al., 2020): where m p is the mass of a particle, ṽ is the translational velocity of the particle, f D is the drag force between gas and particle, g the gravitation, f C the summation of all contact forces exerted on the particle including collision and friction, I p the moment of inertia, x the angular velocity of particle, TT the tangential torque, and Tr the rolling friction torque.In this study, Lagrangian objects are considered soft spheres.The linear-spring-dashpot model is employed for the estimation of pair-wise collision, as described in detail by Goniva et al. (2012).The drag force model proposed by Beetstra et al. (2007) is used for drag closure: where l eff the effective dynamic viscosity, e is the local void fraction, d p the particle diameter, ũ the fluid velocity, and Re p the particle Reynolds number.Note that the sub-grid particle turbulent dispersion is not considered as Olenik (2020) has shown that the dispersion of particles has a negligible effect on particle dynamics.

Particle mass balance
Pulverized coal undergoes a gasification process, in which solid fuel coal gasifies into chemically reactive gases (Ishii, 2000).Gasification involves several sequential sub-processes.First, the coal is heated and the moisture evaporates.Further heating induces devolatilization, in which volatile matter is released from the coal.The remaining char is then converted to the gas phase by the oxidizers.These sub-processes occur sequentially with some overlapping in time.In this study, the gasification process is simplified into two sequential processes, first devolatilization and then char conversion.The mass balance can be described as follows where the m VM;p and m char;p is the remaining mass of volatile matter and char in the particle, respectively.The devolatilization is modeled using the Single First Order Reaction approach by Badzioch and Hawksley (1970): The rate of devolatilization has an Arrhenius form, where A V is the pre-exponential factor, E V is the activation energy, R is the gas constant, and T p is the particle temperature.Eq. ( 6) shows that the initial mass of volatile matter is calculated by multiplying the volatile matter obtained from proximate analysis m VM0 by a factor of Q to compensate for the higher yield of volatile matter at a higher temperature.It is assumed that during the devolatilization process, the particle does not shrink or swell in size.Thus, during the devolatilization, the particle diameter is considered constant.When the remaining volatile matter mass drops below 0.1% of the initial volatile matter mass, the devolatilization process ends and the pulverized coal particles proceed to the char oxidation stage.The 0.1% is an arbitrary value to prevent numerical error and the results are not sensitive to it.The char oxidation is modeled using the Baum and Street approach (Baum and Street, 1971), with the expression given by where p O 2 is the partial pressure of oxygen, K kin is the coefficient of the chemical reaction rate, K diff is the coefficient of the oxygen diffusion rate, and r p is the radius of the particle.This model assumes that the overall reaction includes the effects of both chemical reaction and bulk diffusion of oxygen to the surface of the particle.K kin and K diff can be evaluated using where A char is the pre-exponential factor and E char is the activation energy.The diffusion process of oxygen to the surface of particle is modeled by Eq. ( 9), where T f is the fluid temperature interpolated from the Eulerian grid to the particle location.The value C 1 ¼ 4:998 Â 10 À12 follows the study of Akaotsu et al. (2020) and the unit of C 1 depends on the unit of oxygen partial pressure.Unlike the assumption of constant particle size during the devolatilization, the char oxidation consumes particle mass, and correspondingly the particle size shrinks (thus a constant density).

Particle energy balance
Pulverized coal particles can be heated up or cooled down through multi-phase heat transfer pathways and heterogeneous chemical reactions.The particle energy balance is described as follows where C p;p is the particle heat capacity, and the energy source/sink terms due to Q conv the convective heat transfer between gas and solid phase, Q devol devolatilization process, Q char char surface reaction, and Q rad the radiation.For the convective heat transfer, where k f is the thermal conductivity of the fluid, A p is the particle surface area, and Nu is the Nusselt number which is calculated using Gunn's model (Gunn, 1978).The radiative heat transfer is expressed as, The P1 approximation is used to estimate the incident radiation (Modest, 2013).In Eq. ( 12), e p is the emissivity of the particle, r is the Stefan-Boltzmann constant, and G is the incident radiation term.The incident radiation is transported according to the following equation, where the a g ; a p ; s p , and E p , are the absorption coefficient of fluid, the absorption coefficient of particle, the scattering coefficient of particle, and the emission from the particles to the fluid, respectively.a g is assumed to have a constant value of 0.075 (Akaotsu et al., 2020).a p ; s p , and E p are calculated on the Eulerian cell where the particles locate, where N p is the number of particles in cell, V cell the cell volume, and A pp is the projection area of a particle.Finally, following the previous study (Rieth et al., 2016), the reaction heat of devolatilization Q devol and char oxidation Q char in Eq. ( 10) is assumed to be zero and À9.21 MJ/kg, respectively.

Coarse-graining adaption
The essence of the coarse-graining method is to reproduce physics and chemistry by tracking fewer Lagrangian objects.In the perspective of solid-phase physics, the energy lost per coarsegrained parcel collision should be consistent with the energy loss of the represented primary particles' collisions.Various scaling methods have been proposed to ensure consistency in collision energy loss (Radl et al., 2011;Bierwisch et al., 2009;Thakur et al., 2016;Sakai et al., 2012).In this study, the scaling law proposed by Radl et al. (2011) is applied.The coarse-grained normal stiffness k parcel ¼ ak p , and the damping coefficient c parcel ¼ a 2 c p , in which a the coarse-graining factor is defined as where d parcel is the coarse-grained parcel size and d p is the primary particle size.From the perspective of gas-solid interaction and chemistry, the scaling law is relatively simple compared to the collision scaling.Take the devolatilization process as an example, considering a parcel representing a 3 particles with homogeneous properties (i.e.density, temperature, velocity, and composition), its devolatilization rate scales a 3 times to match the rate of a 3 primary particles.Thus, the coarse-graining adaption on Eq. ( 5) becomes This scaling law applies to parcel's momentum exchange with fluid as well as heat and mass transfer, as also applied by Wang and Shen (2022), Tausendschön et al. (2020).Table 1 summarizes the coarsegraining adaption of the equations.

Governing equations
The Favre-filtered conservation equations of mass and momentum are used to describe the LES (Dufresnes et al., 2016): where the '' $ " operator denotes the density-weighted Favre average terms, q f the average fluid density, e the cell void fraction, S p the mass production rate from particle chemical reactions, i.e. char gasification, l eff ¼ l þ l t the effective dynamic viscosity where l is the laminar viscosity and l t the subgrid-scale viscosity is closed using a static compressible Smagorinsky model (Smagorinsky, 1963), and f n D is the drag force of particle n exerted on the fluid.The summation over n p includes all particles residing in the cell.(See Table 2)

Flamelets model
For the gaseous reactions, there are two fuel streams originating from the pulverized coal gasification: one fuel stream consists of the volatile matter released from the devolatilization process, and the other fuel stream consists of gaseous products from char surface reactions.We follow the flamelet approach in previous studies (Watanabe and Yamamoto, 2015;Akaotsu et al., 2020;Wen et al., 2018), in which two mixture fractions are introduced to describe the mixing process between volatile matter, char-off gases, and oxidizers.The definition of mixture fractions for volatile matter and char-off gases are given as follows, where m devol is the mass of gas originating from the devolatilization, m char is the mass of gas originating from the char oxidation (char-off gas), and m oxid is the mass of oxidant gas.The transport equations of the two mixture fraction Eqs. ( 21) and ( 22), and progress variable Y PV are solved in the CFD solver to track the state of combustion.
@e q f e Z devol @t þ @e q f u i e where Sc is the laminar Schmidt number, Sc t is the turbulent Schmidt number, e S devol and e S char are the mass release from devolatilization and char oxidation, and e _ x PV is the source of progress variable.
In this study, both Sc and Sc t are set to a value of 0.5.There is no mass release term in the transport equation, Eq. ( 25), of progress variable, because it would stimulate the combustion process.As a result, the choice of progress variable should not contain components from the volatile matter.The definition of Y PV is given as (Wen et al., 2017).For the convenience of tabulation, we used the following transformations for the look-up variables, suggested by Hasse et al. (2005): To construct the FGM look-up table, a collection of one-dimensional non-premixed flamelets are generated using Chem1d (Somers, 1994) with the GRI3.0 mechanism (Smith et al., 1999), which contains 325 reactions and 53 species.Since the enthalpy lost is not considered in the one-dimensional flamelet calculation, enthalpy was introduced as an additional dimension in the FGM library to account for the heat exchange between gas and dispersed phase and between gas and environment.Thus, the transport equation of absolute enthalpy is also solved online, where a f is the thermal diffusivity, Pr t is the turbulent Prandtl number, h devol and h char are the enthalpy change due to devolatilization and char oxidation, respectively.

Convection
Eq. ( 11) Eq. ( 12) Composition of fuel streams volatile matter and char-off gas (Wen et al., 2017 In the one-dimensional flamelet calculation, the different enthalpy level was achieved by changing the boundary temperatures T b .Wen et al. (2018) suggested that the multiphase heat transfer can be covered by equalizing the fuel temperature T f and oxidizer temperature T o at the flamelet boundaries, namely For each value of mixture fraction ratio X and boundary temperature T b , a steady-state condition is applied and the strain rate is varied from 10 to the extinction limit.After the extinction limit has been reached, a transient simulation is performed to capture the quenching behaviour.The parameter X and T range from 0 to 0.9 and from 200 K to 1200 K, respectively.The composition of the fuel can be calculated as , where i denotes the index of chemical species in the volatile and char-off gas.In this study, the volatile matter is assumed to be composed of CH 4 , C 2 H 2 , and CO, whereas the char-off gas consists of CO, as proposed by Wen et al. (2017).
A 3D visualization of the FGM library is shown in Fig. 1a for X ¼ 0:0.The lowest energy level corresponds to the flamelets generated from a 200 K boundary temperature while the highest energy level from 1200 K boundary temperature.When storing the tabulation data, the normalized enthalpy is used for interpolation.Since enthalpy is not dependent on the progress variable, as can be observed in Fig. 1a, the enthalpy is normalized with respect to the minimum and maximum enthalpy with given Z and X, It can also be seen in Fig. 1a that the progress variable Y PV is dependent on Z; X, and enthalpy.Following a similar procedure, the progress variable Y PV can be normalized by its minimum and maximum values at a given Z; H norm , and X, The visualization of the To account for the effect of turbulence on combustion, a presumed top-hat PDF approach for Z is employed (Floyd et al., 2009;Olbricht et al., 2012).As a result, an additional pre-integrated table is used for retrieving the integral.

Simulation setup
The CRIEPI pulverized coal jet flame experiment by Hwang et al. (2007,) is chosen for the validation study thanks to its simplicity and sufficient amount of data for quantitative and qualitative benchmark.The simulations are performed using the setup and parameters from the CRIEPI experiment.The experiment operating conditions are summarized in Table 3.The corresponding boundary conditions are given in Table 4.A cylindrical domain is constructed for the simulations, as illustrated in Fig. 2a.The length is 210 mm and the radius is 60 mm.The injection lance consists of a laboratory-scale annular burner with the main port of 6 mm inner diameter and a pilot burner of 0.5 mm co-axially installed outside the main port.The inner wall and the outer wall are 0.5 mm and 1 mm thick, respectively.The pulverized coal carried by air is injected through the main port while the methane is supplied to the pilot burner for the stabilization of the flame.The main port has a length of 2 mm reserved for particle injection.
The mesh consists of 1.8 M hexahedral cells in which three levels of mesh refinement are employed, as shown in Fig. 2b.The center region has the finest level and the mesh becomes coarser toward the outer region.The finest part of the mesh is comparable to the fine mesh used in Wen et al. (2017).The number of grid points at the main port is 22 in the radial direction.Pope's criterion (Pope, 2004) is used to assess the mesh quality.Fig. 3a shows Pope's criterion on the mesh calculated at the pseudo-steadystate.More than 90% of the cells have a Pope's criterion larger than 0.8, indicating the kinetic energy is well resolved.A small portion of the cells has Pope's criterion lower than 0.8.Those cells locate closely to the coflow air boundary, which is outside of the region  of interest due to the absence of flame and particles.Thus, it is concluded that the quality of the mesh is good for the current study.
The carrier air follows a turbulence velocity profile provided by the CBC workshop (Hasse et al., 2019).The provided velocity profile does not correspond to the spatial and temporal discretization of this study.Thus, the velocity is linearly interpolated in space and time.The particle size distribution follows a Rosin-Rammler distribution with the scale parameter of 28 and the shape parameter of 3.07, giving a number-based average diameter of 25 lm and a mass-based average diameter of 33 lm (Hara et al., 2015), as shown in Fig. 3b.The injection process places particles in the 2 mm reserved region in the main port.Within this region, the particles' velocity is interpolated from the CFD cells that they are located in.This method assumes particle has reached terminal velocities and is dominated by the fluid.In the original experiment, methane is supplied from the pilot port.However, the two-mixture fraction model considers volatile matter and char-off gas as fuel.An additional mixture fraction would have to be introduced to cover the methane stream as well as an additional dimension in the tabulation (Wen et al., 2019).In this study, for simplicity, the methane stream is replaced by volatile matter, following the approach in the study of Rieth et al. (2016).
The simulations are realized by combining several open-source toolboxes: OpenFOAM Ò for CFD (Weller et al., 1998), LIGGGHTS Ò for DEM (Kloss et al., 2012), and CFDEM Ò for the multiphase coupling (Kloss et al., 2012).A compressible solver cfdemSolverFGMPimple, which is based on the solver cfdemSolverRhoPimple developed by Lichtenegger (Leitz et al., 2018), has been developed.The simulations were carried out on the Dutch national supercomputer SURFsara Snellius.Each simulation utilized 64 cores for a physical runtime of 0.5 s.

Results
The simulations are run from their initial conditions until reaching a pseudo-steady-state.A sampling window of 0.5 s is used for the statistical sampling (Stein et al., 2013).As the experimental  data was measured along the centerline parallel to the particle jet direction, the numerical sampling for the continuum phase is performed by time-averaged the LES-filtered data along the centerline, whereas for the discrete phase, a cylindrical sampling probe with a diameter of 1 mm is used.To ensure the statistic can capture the high-frequency features, a 1 ms sampling frequency is used for the discrete phase and 0.5 ms is used for the continuum phase.

Non-reacting flow
For the validation of the CG-EL-FGM model, three CG factors a = 1, 2, and 3 are considered.Note that a = 1 represents the non-CG case or the unscaled particles and acts as a reference case.It is prefixed with a for simplification.The simulation results using different a are expected to align with each other, indicating the validity of the CG method when the computational cost is reduced.In the non-reacting flow, although the flame is not ignited and the temperature remains at the ambient, the FGM look-up process is still considered to observe the mixing behavior of fuels and air.The mixture of fuels and air is more viscous than the air, causing the dynamic viscosity to increase and reach an asymptotic value along the centerline.
The snapshot of the non-reacting flow captured at pseudo steady state is shown in Fig. 4. As can be seen in the fluid velocity contour, the jet leaves the main port with high velocity until breakup, where the turbulence disrupts the integrity of the jet.The velocity after break-up quickly reduces and the turbulence disperses further downstream.The particle jet also behaves in a similar manner.Upstream particles remain mostly concentrated and after the jet break-up, a plume of particles forms.To quantitatively compare with the experimental results, the particle axial distributions of mean velocities v p and rms velocities v 0 p from 3 different a are shown in Fig. 5.In general, the profiles follow the trend from the experiment well and there are no significant differences with the use of the coarse-graining method.An offset upstream (x   < 40 mm) can be found in the mean velocity compared to the profile from the experiment.This offset can be attributed to the artificial inflow generation method taken from the CBC workshop (Hasse et al., 2019).The increase of v 0 p when x < 60 mm indicates the development of jet turbulence upstream.After that, v 0 p drops due to the break-up of the jet.
In Fig. 6, radial profiles of the mean and rms velocities are compared.The velocities from all a cases agree well with the results of the experiment.No clear distinction is found between different a, except for the rms velocities at the axial location x ¼ 60 mm and the radial location r > 12 mm.The reason for this slight deviation in rms velocities is because only a very small number of particles is sampled over the period (at r ¼ 13 mm, a ¼ 3 samples 27 particles; whereas at r ¼ 4 mm, a ¼ 3 samples 4299 particles).Before the break-up of the jet, particles are concentrated around the center.As a result, fewer particles can be sampled in the outer region.The lack of sampling data gives uncertainty to the prediction.It can be observed that when the particle jet travels further downstream, the particles spread toward the outer region.The statistics at radial location r > 12 mm then become more significant.
Concluding from the above results, the CG-EL-FGM model can describe the cold flow jet break-up behavior well.The application of the CG method reduces the computational time in 10 folds from a ¼ 1 to a ¼ 3 while still reproducing an accurate particle dynamics and flow field.The good performance in the cold environment concludes the validity of the CG-EL-FGM model in predicting hydrodynamics.
3.2.Reacting flow: reference case with primary particles Fig. 7 shows the snapshot of the particle temperature, gas temperature, mixture fraction Z devol , and mixture fraction Z char captured from the reference case (a ¼ 1) in which the primary particles are simulated.Figs.7a and b show the temperature distribution of the particles and gas.It can be seen that an annular high temperature region is formed due to the combustion of volatile matter and char-off gas, whereas the temperature in the center region is relatively lower.Further downstream the inward diffu-sion of high temperature gas as well as the turbulent fluctuation cause the increase of center temperature.Particle temperature, consequently, increases through multiphase convection and radiation.Figs.7c and d show the mixture fraction of volatile matter Z devol and char-off gas Z char , respectively.It can be observed that upstream Z devol is distributed in the annular region.Pulverized coal in the outer ring of the jet that gains heat from the flame exhibits a faster devolatilization rate.Further downstream, pulverized coal located in the jet center has a high enough temperature for devolatilization, whereas particles in the outer region have depleted volatile matter and, thus, undergo char surface reaction.As a result, Z devol is higher in the downstream center and Z char is higher in the outer region.
Fig. 8 shows the comparisons of the axial distribution of mean and rms velocities between the experiment and numerical results.Results from the literature that also incorporated FGM in LES are included for reference.It is found that the upstream mean velocity v p is well predicted but a discrepancy can be found in the downstream location.The decrease of v p can be explained by the more viscous fuel-air mixture present downstream, as observed in Fig. 7c.The rms velocity v 0 p in the very upstream (x < 20 mm) is under-predicted due to the particle injection method.The pulverized coal particles are placed at the tip of the supply port and are imposed with a velocity interpolated from the cell where they are located.Since the inlet condition provided by the CBC workshop (Hasse et al., 2019) was used and the prescribed fluctuation level is 0.32 m/s, the injected particles also possess the same fluctuation.The rms velocity v 0 p increases to a peak value at x = 50 mm corresponding to the same location where maximum v 0 p is observed in the cold flow.The sampled location is along the centerline where the upstream particles are cold and less reactive.As can be expected, the upstream particle jet resembles the behavior of the cold flow case.The more viscous fuel-air mixture downstream suppresses the development of turbulence and, thus, the v 0 p also decreases.
In Figs.9a and b, the comparisons of O 2 and CO 2 concentration along the centerline are shown.The numerical solutions from this study are compared to solutions in the literature (Wen et al., 2017;Akaotsu et al., 2020;Rieth et al., 2017).The prediction of O 2 from this study agrees with the measurement from the experiment.However, the prediction of CO 2 does not, which is in alignment with what other studies have found.It is considered that in the original experiment, the CO 2 concentration is overestimated due Fig. 8. Axial distribution of particle mean velocities vp and rms velocities v 0 p from the experiment, literature (Wen et al., 2017;Rieth et al., 2017;Akaotsu et al., 2020), and the primary particle simulation a ¼ 1.The solid blue line represents the mean velocity; the dashed blue line represents the rms velocity.to the large sampling resolution as well as a portion of CO is converted into CO 2 (Hara et al., 2015).It has been shown that the current model combining CFD-DEM with FGM succeeds in predicting the dynamics and chemistry of the pulverized coal combustion.Next, the coarse-graining method is applied and the results are discussed in the following section.

Reacting flow: coarse-grained parcels
The performance of the coarse-graining method is evaluated in this section.The coarse-graining factors a ¼ 2 $ 4 are applied to the original system.The performance is considered good if both the physical and chemical responses of the coarse-graining are  consistent with that from the original (a ¼ 1).In Fig. 10a the axial distribution of particle/parcel mean velocity v p and rms velocity v 0 p are compared for different coarse-graining factors.Overall, the mean particle velocity profiles from all coarse-grained simulations a ¼ 2 $ 4 agree very well with the un-coarsened simulation a ¼ 1.However, deviations to the a ¼ 1 case can be observed for the resulted v 0 p when the a increases, especially when x > 50 mm.This deviation can be explained when we take a close look at the viscosity for both laminar and turbulent components along the axial direction, shown in Fig. 10b.It can be seen that when x > 50 mm, the laminar viscosity (solid lines) increases with the increasing a.Higher viscosity not only causes the drop of mean velocity but also suppresses the development of turbulence, which can be implied from the decreasing turbulent viscosity (dashed line) as a increases.Consequently, the v 0 p downstream is lower than the original system, as shown in Fig. 10a.Note that the turbulent viscosity is one order magnitude smaller than the laminar viscosity due to the weak turbulence in this system (Reynolds number 2544).The over-prediction of the laminar viscosity from larger a is a result of concentrated fuel release problem when coarse-graining is activated.The next subsection further discusses this phenomenon in detail.

Coarse-graining-induced localized concentrated fuel release
Fig. 11 shows the radial distribution of particle/parcel temperature and the scaled number of particle/parcel sampled from three axial locations x = 40, 60, 120 mm.It can be observed that the upstream x = 40 mm particle/parcel temperature and scaled number are consistent among all four simulations.This suggests that the application of the coarse-graining method has little impact on the upstream region.However, in the outer region where the devolatilization rate is faster due to high temperature, parcels release volatile matter in a more concentrated manner.Fig. 12 illustrates the concentrated fuel release phenomenon.The essence of the coarse-graining method assumes that a parcel with a size of a times the primary particle diameter can represent a 3 number of primary particles.In other words, all of the a 3 particles are virtually located in the same spatial location.As a consequence, the volatile matter is released at the same cell where the parcel resides, which creates a chemical composition that is different than the case without using coarse-graining method.The higher the a is applied, the more prominent is the localized concentrated fuel release phenomena.Fig. 13 gives the distributions of the devolatilization rate from each a as well as the iso-contour of the mixture fraction Z devol at stoichiometric mixture fraction Z st ¼ 0:096.As can be seen, the devolatilization rate for higher a exhibits a spotty pattern and a higher magnitude, which results in a higher fraction of Z devol .When the volatile fuel is transported, this high gradient can cause a faster diffusion toward the center region.Eventually, more volatile fuel is present in the center region for a higher a value.This explains why in Fig. 10b the flow is more viscous downstream when a increases.Note that the viscosity in a ¼ 4 peaks at around x = 90 mm due to the depletion of volatile matter in pulverized coal.After that, the char conversion releases CO which is less viscous than the volatile matter components.
Another effect of the localized concentrated fuel release for higher a is that the gas in the downstream center region is thermally more conductive and diffusive.Fig. 14 shows a snapshot of thermal conductivity k f and thermal diffusivity a f for coarsegraining factor a ¼ 1 and 4 captured at pseudo steady-state.Note that both thermal conductivity and thermal diffusivity are quantities interpolated from the FGM database during runtime.These values correlate to the species composition.For a ¼ 4, in the downstream center region, because of the concentrated fuel release phenomenon, the gas is thermally more conductive and diffusive.The increase in thermal conductivity enhances the convection heat transfer, see Eq. ( 11).The increase in diffusivity accelerates the transport of enthalpy between the outer flame region and the center region.As a result, the increase in both conductivity and diffusivity leads to a higher parcel temperature in the center region r < 5 mm at x = 120 mm, as can be seen in Fig. 11a.
To mitigate this concentrated fuel release from the parcel concept, smoothing strategies can be applied.The application of smoothing has been a standard for the CFD-DEM method (Tausendschön et al., 2020;Capecelatro and Desjardins, 2013;Pirker et al., 2011;Radl et al., 2014;Nijssen et al., 2020) as the particle size becomes comparable to or larger than the fluid grid size.It improves the numerical stability as well as prevents local extremity.When it comes to the coarse-graining method, smoothing not only inherits the benefits of stability but also tries to reproduce the particle cloud information, as illustrated in Fig. 12c.The coarse-grained parcel represents a cloud of primary particles.As the goal is to replicate the dynamics of primary particles, the smoothing method can virtually spread the presence of particles from parcels.The smoothing operation is often realized by solving a diffusion equation for any quantity / where the diffusion coefficient is defined as and L smooth is the characteristic smoothing length.The advantages of this diffusive smoothing are its conservative nature and its wide applicability to different meshes.Conventionally, the smoothing length is chosen as L smooth ¼ 3d p for mono-disperse particulate flow (Radl et al., 2014).However, what the parcels' smoothing length should be is unclear.In this study, a smoothing length of 3 max d parcel À Á is arbitrarily used for a > 1 and the exchange fields, i.e. void fraction, drag force, convection heat transfer, and mass transfer, are smoothed.Fig. 15 presents the comparisons of the O 2 and CO 2 mole fractions along the centerline for simulations with and without smoothing.The solid lines display results from different coarsegraining factors without using smoothing, whereas the dashed lines display the results from simulations that enable the diffusive smoothing method, as described in Eq. ( 31).In the case of nosmoothing simulations, the response of applying larger a deviates from the reference case after x > 50 mm.When a ¼ 4, the middlestream (50 mm < x < 100 mm) O 2 concentration plummets owing to the fast release of volatile matter.Furthermore, because CO is one of the components of volatile matter, the fast reaction between CO and the rich O 2 in the center jet consumes O 2 as well, leading to both the decrease of O 2 concentration and the increase of CO 2 , as can be seen in Figs.15a and b.In Fig. 15b, the overall under-prediction of CO 2 concentration is attributed to the sampling resolution, as discussed in Fig. 9 in Section 3.2.The predictions from a ¼ 2 are in good agreement with the reference case and the experiment.Once the diffusive smoothing is enabled, the response of both O 2 and CO 2 converge slightly toward the reference case.The use of constant diffusivity implies that the area of the results from simulations that enable the diffusive smoothing method, as described in Eq. ( 31).The smoothing length of 3 max d parcel À Á is used for a > 1. influence is ubiquitous for all sizes of parcels.On one hand, the smoothing operation relaxes the concentrated fuel release for large parcels, spreading information as if there were a particle cloud; on the other hand, because of the universal smoothing length, it accelerates the transport of fuel for smaller parcels, leading to an overspreading.The acceleration of the fuel transport for smaller parcels might counteract the benefit that smoothing brings.The chemistry response, as a consequence, improves only slightly compared to the results without applying smoothing.Nevertheless, the prediction is improved compared to the simulations without using smoothing.Smoothing is, thus, recommended even when the parcel size is smaller than the grid size for the sake of proper thermochemical response (in the case of a ¼ 2, the largest parcel size is 140 lm and the smallest grid size is 250 lm.).However, the exact smoothing length is subject to further study.

Radial distribution of particle dynamic
Fig. 16 shows the radial distribution of particle mean velocities v p and rms velocities v 0 p sampled at three axial locations (60, 120, 180 mm).For the v p values, in general, the predictions from the CG-EL-FGM simulations follow the experimental data well, despite of a small underestimation downstream.Rieth et al. (2017) suggested that the differences in heat and mass release from devolatilization and volatile combustion can cause a change in v p .On the other hand, the results of v 0 p show a deviation from the experimental data in the center and outer region.Such deviation has been reported in other LES (Franchetti et al., 2013;Akaotsu et al., 2020;Rieth et al., 2017) and DNS (Hara et al., 2015) studies as well.In the current study, the overall trend of the radial distribution of v 0 p resembles the results of Akaotsu et al. (2020).Possible explanations for this deviation between simulation and experimental results are: 1) due to the insufficient mesh resolution of the shear layer for resolving turbulence, and 2) because of the exclusion of the particle turbulent dispersion model.The particle turbulent dispersion is deemed negligible because both the axial and radial distribution of particle velocity v p are insensitive to the dispersion coefficient C 0 , as reported by Olenik (2020).However, he also showed that the v 0 p in the outer region is sensitive to the C 0 .The larger the C 0 (< 10 4 ), the higher values of v 0 p in the outer region.This might suggest the necessity of a turbulent dispersion model on particle motion for a better prediction.One is referred to the work of Olenik for further detail (Olenik, 2020).

Computation time reduction
Fig. 17 shows the particles/parcels' spatial distribution colored with velocity.It visualizes the impression that the coarsegraining method can largely reduce the number of Lagrangian objects being simulated.Fig. 18a shows the wall time of each simulation, and Fig. 18b shows the computational time speedup in the portion of DEM and coupling.The simulation of using primary particle a ¼ 1, costs 76 h while a ¼ 2 takes merely 33 h.The significant reduction in computational time can be attributed to the speedup in DEM and coupling.Higher coarse-graining factor (a > 2) brings marginal benefit in time reduction since the simulation is bottlenecked by the continuum phase computation.When a > 2, both discrete phase dynamics (Fig. 10a), and chemistry responses (Fig. 15), are inconsistent with the primary particle case.It is thus concluded that in this CRIEPI pulverized coal burner configuration without the application of a smoothing strategy, a ¼ 2 is the optimal coarse-graining factor.It brings the most substantial cost reduction while preserving a good accuracy in particle dynamics and gaseous phase chemistry.

Conclusion
To achieve efficient Euler-Lagrange modeling of large-scale reactive particulate flows with complex (combustion) reactions, reduction methods for effective modeling of both reactions in the continuum phase and particle dynamics in the discrete phase are important.In this work, we incorporate FGM and CG methods simultaneously in Euler-Lagrange simulations for the above objective.The capability of the current CG-EL-FGM model is demonstrated using the CRIEPI burner case for pulverized coal combustion.The gasification process of the pulverized particle is modeled by two processes in serial, the devolatilization and char surface conversion.The multiphase heat transfer includes the convection and radiation between fluid and gas.The gas-phase combustion is modeled using the FGM method, where the properties of gases are interpolated from a look-up table.
The CG-EL-FGM model is first validated in the non-reacting condition.It is shown that the use of the CG method can yield a consistent result compared to the reference case and the results are in good agreement with the experiment.However, in the reacting flow, only the results from CG factor a ¼ 2 match the reference case.When a larger than 2 is applied, the localized concentrated fuel release causes a different spatial fuel composition, indicating the inconsistency in flame structure.To remedy this concentrated fuel release consequence, smoothing of the parcel influence is found to be necessary even when the parcel size is smaller than the grid size.When a ¼ 2 is used, the computational cost is reduced by a factor of 2.3 and it produces an accurate prediction of particle dynamics and oxygen concentration.It is thus concluded that for the future study on the CRIEPI burner validation, an optimal coarse-graining factor a ¼ 2 can be applied for the cost reduction.The developed model is applicable to large-scale simulations of solid-fuel combustion processes.For large-scale simulations, an optimal coarse-graining factor needs to be considered in terms of the representatives of statistics of particle dynamics.In a future study, we will use this model for simulations of the industrial-scale blast furnace as an example.

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.

Fig. 1 .
Fig. 1.Visualization of flamelets library when X ¼ 0. (a) flamelets temperature as a function of Z, enthalpy H, and PV.(b) flamelets temperature as a function of Z, normalized enthalpy Hnorm, and normalized progress variable YPV;norm.

Fig. 2 .
Fig. 2. (a) Schematic of computational configuration.The dimension is not up to scale.(b) Cross-section of the computational grids.

Fig. 3 .
Fig. 3. (a)  Histogram of the Pope's criterion(Pope, 2004).(b) The mass-based probability density function and cumulative density function of initial particle size distribution.The solid line follows a Rosin-Rammler distribution with the scale parameter of 28 and the shape parameter of 3.07.The circle indicates the sampling particle diameter and the interval of 2.5 lm is used.The dashed line indicates the sampled cumulative density function.

Fig. 5 .
Fig. 5. Axial distribution of particle mean velocities vp and rms velocities v 0 p .The solid lines represent the mean velocity; the dashed lines represent the rms velocity.

Fig. 6 .
Fig. 6.Non-reacting flow result: radial distribution of particle mean velocities vp and rms velocities v 0 p at three axial locations (60, 120, 180 mm) from the experiment and LES simulations.Top row: mean velocities; Bottom row: rms velocities.

Fig. 7 .
Fig. 7. Snapshot of transient quantities captured from reference case a ¼ 1: (a) Particle temperature.Particles are enlarged 4 times their size for visualization.(b) Gas temperature.(c) Mixture fraction Z devol .(d) Mixture fraction Z char .

Fig. 10 .
Fig. 10.(a) Axial distribution of particle mean velocities vp and rms velocities v 0 p from the experiment and coarse-graining simulations a ¼ 2 $ 4. The solid lines represent the mean velocity; the dashed lines represent the rms velocity.(b) Axial distribution of gas phase dynamic viscosity.Solid lines represent the laminar dynamic viscosity; dashed lines represent the turbulent dynamic viscosity.

Fig. 12 .
Fig. 12. Illustration of the concentrated fuel release phenomenon.The light blue color indicates lower magnitude of fuel release, whereas the deep blue color indicates higher magnitude of fuel release.(a) Primary particles with particle size dp.(b) Coarse-grained parcel with d parcel ¼ 3dp.One parcel is a representative of 27 particles.(c) Coarsegrained parcel with smoothing on the fuel release term.

Fig. 13 .
Fig. 13.Snapshot of devolatilization rate for each a captured at pseudo steady-state.The solid line indicates the iso-contour of the mixture fraction Z devol at stoichiometric mixture fraction Zst ¼ 0:096.

Fig. 14 .
Fig. 14.Snapshot of (a) thermal conductivity k f and (b) thermal diffusivity a f for a ¼ 1 and 4 captured at pseudo steady-state.Thermal conductivity and thermal diffusivity are both interpolated from the FGM database during runtime.

Fig. 15 .
Fig. 15.Axial distribution of O 2 and CO 2 from each coarse-graining factor a. The solid lines display the results from different coarse-graining factors.The dashed lines display

Fig. 16 .
Fig. 16.Reacting flow result: radial distribution of particle mean velocities vp and rms velocities v 0 p sampled at three axial locations (60, 120, 180 mm) from the experiment and LES simulations.Top row: mean velocities; bottom row: rms velocities.

Fig. 17 .
Fig. 17.Snapshot of particles/parcels velocities cross-section from different coarse-graining factors.From left to right, a = 1, 2, 3, and 4, respectively.The values indicate the average number of particles/parcels at pseudo-steady-state. Particles and parcels are enlarged 8x for visualization.

Table 1
Coarse-graining adaption on mass and heat transfer and the drag force closure. ).

Table 3
Experimental operating conditions(Hwang et al.,  2007).The flow rates are used for simulations.

Table 4
Simulation boundary conditions.