Grid-independent Eulerian-Lagrangian approaches for simulations of solid fuel particle combustion

In this study, a computational fluid dynamics (CFD) model with three coarse graining algorithms is developed with the implementation of a layer based thermally thick particle model. Three additional coupling methods, cube averaging method (CAM), two-grid method (TGM) and diffusion-based method (DBM), are implemented. These coupling methods are validated and compared with the widely used particle centroid method (PCM) for combustion of a biomass particle in a single particle combustor. It is shown that the PCM has a strong dependence on the grid size, whereas the CAM and TGM are not only grid independent but also improve the predictability of the simulations. Meanwhile, a new parameter, the coupling length, is introduced. This parameter affects the sampling of the gas phase properties required for the particle model and the distribution of the solid phase properties. A method to estimate the coupling length by using empirical correlations is given. In general, it is found that a too small coupling length underestimates the heating-up rate and devolatilization rate, while a too large coupling length overestimates the O2 concentration at the particle surface. The coupling length also has an influence on the distribution of the gas phase products.


Introduction
Direct combustion of solid fuels, such as coal and biomass, is one of the main routes to generate heat and electricity [1,2]. An improved design of the combustor can increase the combustion efficiency and reduce emissions. With the rapid development of the computer hardware and numerical methods, computational fluid dynamics (CFD) with increasingly detailed sub-models is widely adopted by the industry as a powerful analysis tool to reveal the details of chemical and physical processes involved [3,4].
To model the multi-phase combustion system in CFD simulations, the gas phase is usually described with a continuum approach in the Eulerian framework. The solid phase is treated either as a continuum in the Eulerian framework, or as a dispersed phase by discrete methods in the Lagrangian framework, in which the particles are tracked individually. Irrespective of the framework, single particle conversion models are normally required as sub-models to describe the thermal decomposition of the solid phase. These include the sub-processes of heating, drying, devolatilization and char burnout. Such single particle models use local operating conditions from the gas phase to predict heat and mass release from the particles as boundary conditions or source terms for the gas phase.
The particle models and gas phase models are developed under different frameworks. When linking the two models together, bridging between them is critical. As shown in Fig. 1, the particle model requires information about its local gas field quantities as its boundary conditions, and the effects of the presence of the particle should be transferred back to the governing equations of the gas phase through a source term. For Eulerian-Eulerian approaches, although the massive computation of the particle-particle interactions is avoided, it is difficult to account for the distribution of different particles. In fixed bed simulations, Ström et al. [5,6] registered the particles into individual computational grids. They assumed that particles within one grid cell have the same degree of conversion. Such a method converts the particle sub-models into fully Eulerian models, but also makes it impossible to predict the deformation of the fuel bed from first principles. Some efforts have been made to account for this in similar studies, where extra grid transformation models are employed to include shrinkage of the fuel particles [7]. In contrast to Eulerian-Eulerian approaches, Eulerian-Lagrangian approaches on the other hand do not need extra averaging models for the solid phase properties, and can provide more detailed information for individual particles. This approach is more widely used when coupling with thermally thick particle models [8][9][10][11][12], revealing that the intra-particle temperature gradient has significant influence on the conversion process [13]. However, in the conventional discrete particle model (DPM) or discrete element method (DEM), the discretization of the governing equations of the continuous phase and the Lagrangian tracking of the dispersed phase employ the same grid system, and the coupling (heat, mass and momentum exchange) between the particle and the gas phase only happens inside the cell in which the particle's centroid is located. Such a coupling is called the particle centroid method (PCM). PCM requires that the particle length scale is much smaller than the grid size and that interactions between the particle and gas phase should not be significant [14].
A more direct way to couple the single particle to the gas phase in simulations is to resolve the boundary layers around the particle with a body-fitted grid at the gas-particle interface [15,16]. Considering the complexity and scales in typical industrial applications with many particles, this approach is not practically feasible. Instead, one-dimensional particle models formulated with uniform gas phase conditions are usually used. The heat and mass transfer from the gas phase to the particle surface are estimated by correlations between dimensionless numbers, such as Nusselt number (Nu) and Sherwood number (Sh). When coupled to the gas phase, the particles are treated as Lagrangian point-particles, which means that the particles' boundaries are unresolved and its geometry is neglected by the gas phase. This situation presents opposing requirements on the grid resolution from the two combined frameworks: the grid used for the gas phase must be fine enough for the solution of the governing equations to be grid-independent, while the particle conversion model requires that the grid should be large enough to allow for proper averaging to obtain the interphase properties, such as the porosity and the variables of local gas phase solution and physicochemical properties. When the size of the reacting particle is relatively large compared to the grid size, the opposing requirements on the grid cause conflict between the models' linkage. Besides, the reacting particle has strong interactions with its local gas phase, which are presented in source terms. Large positive source terms will strongly exacerbate solver robustness, giving rise to code errors and unreliability [17]. In direct numerical simulations (DNS) of conversion of pulverized coal particles, Krüger et al. [18] used Laplacian diffusion to diffuse the source terms produced from the particle model before addition to the gas phase governing equations in order to increase the numerical stability. Similarly, Farazi et al. [19] used a Gaussian kernel function to redistribute the source terms in a simulation of ignition and combustion of coal particles. In their work, the grid size is equal to coal particle's diameter. In order to get the particle's boundary condition, the gas phase properties were averaged from a cube at the particle's location. However, it is still not clear to what extent, if any, the parameters introduced by the coupling method influence the final results.

Nomenclature
In fact, the conflicting requirements on the grid is a common issue in several research communities when simulating relatively large particles. Sun et al. [20] made a brief review over these so-called "coarsegraining" methods, mapping from the particle-scale quantities to macroscopic fluid field quantities, and summarized the four main methods: PCM, the divided particle volume method (DPVM), the statistical kernel method and the two-grid formulation. As mention above, the conventional DPM uses PCM. In DPVM, the particle's volume is divided among all the cells it overlaps, so that each cell receives the actual volume inside it. Neither PCM nor DPVM completely resolves the underlying conflicting theoretical requirements for the grid resolution posed by the gas and particles phases. The statistical kernel method uses kernel functions, for example the Gaussian distribution function, to redistribute the solid phase properties to the computational domain. The two-grid formulation resolves the solid phase and gas phase under different grid systems with proper field mapping methods between different grids. These two methods are able to address the deficiency of the PCM in the cases of small cell size to particle diameter ratios.
Furthermore, Link et al. [21] proposed a porous cube representation method for a simulation of a spout-fluid bed. Every particle was represented by a porous media cube proportional to its own size when coupled to the CFD simulation. However, most of these studies focused exclusively on the hydrodynamics. To the authors' knowledge, there are no systematic studies on the influence of the coupling scheme on the predictions of solid fuel particle combustion. It is to be expected that the method employed to couple the Eulerian and Lagrangian frameworks will play an even more critical role in the presence of significant heat and mass transfer and chemical reactions, particularly due to the strong non-linearity of the latter.
The objective of this work is to study the coupling effects when simulating the combustion of solid fuel particles using coarse-graining methods. The porous cube representation method, the two-grid formulation, and a diffusion-based method, which is theoretically equivalent to statistical kernel method, are extended with reacting particles in this work. The mass and heat transfers are included under the same principle as momentum transfer in the original works. This study is focusing on combustion of thermally thick biomass particles modeled by a computationally efficient particle model [22], as biomass particles are normally larger in size compared to, for example, pulverized coal particles. The implemented coupling methods are compared and further discussed together with the conventional PCM through the simulation of a single particle combustor. Meanwhile, a method to estimate the additional coupling parameters, based on the physical non-dimensional Sh and Nu numbers, is proposed.

Mathematical modeling
The Eulerian-Lagrangian solver developed in this work is based on OpenFOAM. The gas phase is solved by using the Reynolds Averaged Navier-Stokes (RANS) equations, and a standard kmodel is used to account for turbulence [23]. The governing equations are summarized in Table 1. Here, is the volume fraction of the gas phase or the porosity; S S S , , m h U and S i are the source terms calculated from the single particle model, and Q gas and S Ri are the reaction heat and species source terms respectively due to the homogeneous gas phase reactions. Biomass particles are modeled by Lagrangian tracking scheme. Thermochemical degradation and conversion of the particles are calculated by a thermally thick single particle model, with the boundary conditions obtained from the solutions of the gas phase equations as prescribed by the coupling method in question.
The thermophysical properties of the gas mixture such as conductivity, thermal diffusivity and viscosity, are calculated by massweighted mixing laws. The ideal gas law is used to calculate the density of the gas phase. The effective dynamic thermal diffusivity eff and mass diffusion coefficient D eff for the species are calculated through the turbulent Prandtl number (Pr t ) and the turbulent Schmidt number (Sc t ), respectively [23]. The P-1 model is used as the radiation model. It is the simplest case of the more general P-N model and is formulated by the partial differential equation in incident radiation [24].

Layer based particle sub-model
The layer-based single particle model proposed by Ström et al.

Table 1
Gas phase governing equations.

Energy equation
Species transport equation  [5,6], which is based on Thunman's approach [22], is selected to describe the thermally thick particle. As shown in Fig. 2 the spherical particle is divided into four distinct layers: wet wood, dry wood, char and ash. The thermal conversion is assumed to occur at the infinitely thin fronts between the layers. For other typical biomass particle shapes, the model also applies if the surface area can be expressed as a function of the distance to the center. For non-spherical shapes, the temperature gradient along the radius direction could be corrected by the suggestions given by Ström et al. [6]. Each inner boundary between the layers is assigned a temperature, which is the reaction or conversion temperature for particle sub-models. The outermost boundary is the particle surface, and its temperature is determined by a balance between the gas phase and the intra-particle heat transfer process. The model can be simplified into a 1D discrete model along its radial direction. The intra-particle temperature gradient is predicted by resolving the heat conduction inside the particle. The layer mass is updated according to the reaction rate. The full details of the heat and mass transfer model is provided in Appendix A.

Devolatilization model
A two-stage wood devolatilization model is used in this study, shown in Fig. 3 [25]. Dry wood is converted into light gases, tar and char through three competing parallel reactions. Parts of the tar is further converted into light gases and char in the second step of the reaction, which is considered to occur inside the particle. By using this model, the char yield is determined by the temperature history. The light gases have a presumed composition which is listed in Table 2. In the simulations, the light hydrocarbons in the gas are represented by methane. The tar consists of heavy hydrocarbons which are lumped into a representative molecule C 6 H 6.2 O 0.2 , and its properties are given by those of benzene [26].
The kinetic constants are calculated by the Arrhenius expressions shown in Eq. (5). The kinetic data are listed in Table 3.
The heat balance of devolatilization includes the exothermicity of char formation and the endothermicity of generation of volatiles [27]. In this study, the devolatilization is considered as a heat neutral process, which means that it is neither exothermic nor endothermic.
The layer model assumes that the devolatilization also occurs in an infinitely thin front. However, the reactions could take place in a rather wide temperature range. The volumetric reactions are used instead of    Table 3 Kinetic data of reactions in particle sub-models.
Ref. the surface reactions to correct for this, by assuming that the temperature is piecewise linear between T b0 and T b1 [22].

Char conversion model
Char conversion reactions are heterogeneous reactions. The reactant gases reach the char surface by diffusion and convection. Thunman's model for char conversion [22] is selected in this study. The four main reactions and their rate equations are listed in Table 3. The char reaction process is a diffusion-controlled process. Hence the effective char conversion rate (R char i , ) also considers the mass transfer effects in Eq (6).
Here, C i is the species concentration in the particle's surrounding gas phase, which is calculated according to the coupling schemes. k r i , is the kinetic rate, which also follows the form of Eq. (5). Since the ash layer is also considered in this study, the diffusion rate k d i , has contributions from both the diffusion of gases to the particle surface and the diffusion through the ash layer. The mathematical framework describing this process is given by Eqs. (33)-(37) in a previous work [28], except Eq. (35) which is replaced by the Ranz-Marshall correlation [29] in this study.

Homogeneous gas phase reactions
The combustible gases released to the gas phase from devolatilization and char gasification participate in homogeneous reactions. A global reaction scheme is used in this work [7], and the kinetics are listed in Table 4. To account for the effect of turbulence, the partially stirred reactor (PaSR) combustion model is employed. The species mixing time scale is calculated from the turbulent properties. The reaction rate is adjusted according the reaction time scale and mixing time scale according to Mohseni et al. [30].

Coupling between gas and particle
The coupling scheme between the gas phase and the particle submodels should provide two things: the properties of the gas as seen by the particle, and the effect of the particle as seen by the gas. Therefore, a coupling scheme consists of two components. The first is used to obtain the local gas phase properties, which provide the boundary condition for the particle sub-models. The second is to distribute the source terms from the sub-models and update interphase information such as the phase volume fraction. In the conventional PCM, the coupling occurs only between a grid cell and the particles whose centroids fall within the cell. The gas properties are represented by the average values of the particles' host cell, or interpolation values. The gas phase fraction and source terms for the host cell are calculated as: where j is particle index. Such an arrangement is neither valid nor stable when the computational cell, in which the particle resides, approaches or even becomes smaller than the particle itself. In order to study the coupling in more detail, three other coupling strategies are used and described in the following sections.

Cube averaging method (CAM)
In order to overcome the grid dependency in the hydrodynamic simulation of a fluidized bed, Link et al. [21] proposed a porous cube model. As shown in Fig. 4, instead of directly coupling the particle to its owner grid, a cubic region is created as an interaction media between the particle and gas phase. By doing this, the Lagrangian point particle is transferred into an Eulerian porous media. The original work only focuses on the momentum transfer, and the calculation of solid volume fraction. In this work, all the mass and heat transfer terms due to thermal conversion of the particle are also coupled to the gas phase through the porous media cube.
The side length of the cube d s depends on the particle diameter d p and a constant factor a: where a is a free parameter. In the hydrodynamic simulations, Link et al. [21] used a value of 5. The gas phase property is mapped to the cube based on the volume-average: where f j cube is the volume fraction of grid cell j occupied by the cube.
The cube property cube provides the boundary condition for the particle. All the source terms calculated according to the particle model return their values to the cube. And then the source terms calculated by Eq (11) are mapped back to the grid cells from the cube.
The source terms in the governing Eqs. (5)-(8) are replaced by the results from the above equations. The solid volume fraction is calculated at the cube level, and is also mapped onto the grid in the same way as the source terms. V cube and f j cube need to be updated as long as a new d p or the particle's location is calculated. To calculate f j cube , there are many scenarios that need to be accounted for when dividing fluid cells that intersect with the cubes. To apply CAM on unstructured grid cells, certain interpolation schemes are required, for example, a conservative interpolation scheme developed by Su et al. [31] in which the cells are decomposed into tetrahedrons to calculate the intersection volume. Such an implementation may be tedious work but is fully feasible. To simplify the implementation, a Cartesian grid is used in this study. The cubic shape of the porous media region is also chosen on the basis of this concern.
As for multi-particles systems, the cubes can overlap with each other, due to that the calculation for each cube is independent, and every grid cell is restricted with a maximum solid phase fraction of 0.99. In the near-wall region, the part of a cube exceeding the calculation domain will be discarded before the mapping calculation.

Two-grid method (TGM)
As described above, the theoretical requirements of the grid size from the fluid phase and the particle can be in conflict. The fluid requires fine grids to resolve the flow, however, the particle sub-models require that the grid size should be large enough to represent the local volume fraction of the particulate phase properly. Deb et al. [32] and Farzaneh et al. [33] proposed to use different grids for fluid and particles separately. As shown in Fig. 5, a virtual coarse grid is thereby created based on the fine grid of the gas phase. The coarse grid information is calculated by averaging the gas phase properties, including velocity, temperature, species concentrations and thermal properties. The particle model is resolved on the coarse grid and the source terms are mapped to the fine grid. The phase fraction is calculated on the coarse grid, and the overlapped coarse and fine grid cells share the same phase fraction value. The averaging of the gas phase properties and the reallocation of the source term are weighted by the volume of the fine grid to the coarse grid.
The mapping uses the same equations as in CAM (Eqs. (14) and (15)), but the cube is replaced by the coarse grid cell that the particle locates. There are two main differences between the two methods. The first one is that the particle is always in the averaging region center for CAM, while for TGM, the treatment of the particle at coarse grid is in the same way as PCM, which means that the particle moves in the fixed particle grid. The other difference is that, as the particle shrink during the conversion, the porous cube also becomes smaller, keeping the porosity of the cube constant, until the cube gets smaller than the fluid grid. For the multi-grid method, however, the size of the coarse grid is fixed.

Diffusion-based method (DBM)
Another way to transfer a Lagrangian point particle into an Eulerian field is to use the statistical kernel functions. The particle will be distributed onto a domain according to a weight function called a kernel function h x ( ), as shown in Fig. 6. For example, the solid volume at location x consists of the distributed volumes from all particles, which can be expressed as: For Gaussian distribution, the Gaussian kernel function as shown in Eq. (13) is applied with a free parameter, bandwidth b.
The kernel function method is difficult to implement into a CFD solver, especially when the calculation domain has non-orthogonal boundaries. Capecelatro et al. [34] proposed a method to resolve a diffusion equation of the distributed properties to represent the results using the statistical kernel functions. One main advantage of such a method is that no special treatment of physical boundaries is required. Before updating the solid phase volume fraction and the source terms to the gas phase governing equations, these terms are dispersed by a passive scalar diffusion equation: The diffusion equation is solved from =0 to the time scale =T. Sun et al. [20] proved that the diffusion-based method and the Gaussian kernel based averaging method Eq. (13) are mathematically equivalent when the bandwidth b of the Gaussian kernel function and the diffusion time scale T satisfy = b T 4 . As for the implementation, the redistribution can only be applied to the solid phase properties, and the diffusion-based method is directly inherited from PCM. This means that the gas properties required by the particle sub-models are sampled from the particles' host grid, and the source terms are first calculated under PCM. Then Eq. (14) is solved with OpenFOAM's standard Laplace operator for all the source terms and phase volume fraction. This method smooths the particle's influence on the gas phase, and the changes of the gas phase properties in the region near the particle become more moderate than when using PCM.

Numerical simulation
To understand the effects of heat and mass transfer caused by the thermal conversion of reacting particles on the coupling methods, CFD simulations employing the different coupling approaches are configured according to the experiments from a single particle combustion reactor [35,36]. The geometry of the reactor as well as the calculation domain are shown in Fig. 7. At the bottom, the inlet gas is provided by a flat flame burner. The biomass particle is suspended in the center of the reactor with a distance of 300 mm to the burner. The conversion process of the particle was recorded by a camera. The domain is generated as a cuboid to apply the Cartesian grids. As mentioned above, the Cartesian grids greatly simplify the implementation of CAM, and also help to avoid unintended errors when different coupling methods are employed. Since the reactor is a single particle reactor and this study is focusing on the coupling effects, as long as the near-particle region is well represented, the transformation of the walls into rectangular shapes is of minor significance. The size of the computational domain is 50 mm × 50 mm × 150 mm. Four sets of Cartesian fluid grids are generated with 7 × 7 × 21 (Coarser), 11 × 11 × 33 (Coarse), 17 × 17 × 51 (Fine) and 33 × 33 × 99 (Finer) cells across the entire domain and the side length of the grids are 2.38d p , 1.52d p , 0.98d p and 0.51d p respectively.
The thermophysical properties of the gas phase, as well as the reaction heat of the gas phase reactions, are evaluated by the standard NASA polynomials [37]. The particle's properties are summarized in Table 5. Uniform inlet boundary conditions are used. In order to match the gas phase temperature at the particle's location as given by the experimental measurement, the inlet temperature is set as 1473 K with a uniform wall temperature of 1250 K. The average gas phase velocity is set to 1.38 m s −1 [35] with an oxygen concentration of 20%. Test cases using both the rectangular domain and the domain with the actual cylindrical geometry showed no significant differences with regard to the temperature and incident radiation at the particle location.
The standard PISO (Pressure-Implicit with Splitting of Operators) algorithm is used to calculate the coupling between the velocity and pressure fields. It should be noted that the single particle model and the gas phase solver use different time steps. For the thermally thick particle model, the time step is × 5 10 7 s. Adjustable time steps with Courant number of 0.1 (time steps are in the range between × 2 10 5 s and × 8 10 5 s) is used for the gas phase calculation. A brief numerical scheme is given as follows: • Step 1. Resolve the gas phase governing equations together with homogenous gas phase reactions and update the fluid fields.
• Step 2. Calculate the average gas phase properties according to the chosen coarse-graining method.
• Step 3. Resolve particle sub-models in one particle time step. Update particle information, and restore all the mass, momentum and heat transfer source terms between particle and gas phase. Forward particle time step and repeat step 3, until one fluid time step has elapsed.
• Step 4. Update the gas phase volume fraction field. Redistribute accumulated source terms according to the chosen coarse-graining method. Advance time to the next fluid time step, go back to step 1 and repeat until finished.
The coupling of the source terms is done using a semi-implicit scheme, except in DBM simulations, which use an explicit scheme.
The different coupling methods, except PCM, introduce an additional parameter, which is the side length of the cube d s in CAM, the side length of coarse particle grid x coarse in TGM and the bandwidth b in DBM respectively. These three parameters have a similar physical meaning. They are the length scales in which that the particles can still be treated as point particles. The mass, momentum and heat transfers between the particle and the gas phase can be directly coupled at such scales without resolving the transfer process inside the coupling region. Here, the three length parameters are defined as the coupling length x c . The default values are taken from the recommendations in the reference papers which are 5d p , 5d p and 6d p for CAM, TGM and DBM respectively [21,32,20]. Parameter studies with varying x c have also been documented in a later section.

Grid independence of different coupling methods
Grid independency studies have been conducted with different coupling methods. In Fig. 8, the particle's surface temperature and residual mass ratio, which indicates the extent of the conversion with different grid resolutions, are presented. The particle surface temperature reflects different conversion stages. When the conversion starts, the particle is heated by the gas phase and the drying process begins. After a short period, the devolatilization starts and causes a rapid mass loss. The released gases from the devolatilization also undergo homogeneous gas phase reactions. The clear inflection point on the residual mass ratio curve implies the end of devolatilization. The residual mass then mainly consists of char. The particle surface temperature rapidly increases to its peak due to the char oxidization. Afterwards, the rate of char conversion gradually decreases owing to the shrinkage of the reacting surface at the end of char burnout stage. The residual mass ratio as well   Ash layer porosity 0.65 - [28] Particle emissivity 0.8 - Fig. 8. Surface temperature and conversion ratio of the particle with different coupling methods and mesh resolutions. Solid lines are particle surface temperatures, dashed lines are residual mass ratios. as mass loss rate versus the particle temperature are presented in Appendix B.1, to show the changing of the conversion stages more clearly.
The results from the PCM shows a strong dependence on the grid size. The predicted particle's surface temperature becomes lower with decreasing cell size, and the conversion processes slow down consequently. This happens since the local effect of the source terms from the gas-particle coupling increases with decreasing cell size. It leads to increasingly poorer predictions of the state of the far-field gas phase properties when sampling inside the cell to obtain boundary conditions for the particle conversion models. These observations indicate the pronounced dependence on the coupling scheme, which may significantly influence the dynamics of the conversion process.
The simulations are transient. In order to quantify the deviations of the results between different grids with the same coupling method, the estimated time needed to achieve the same degree of conversion is compared. The deviations are evaluated based on the results using Finer grid and averaged over the entire process. The predictions of the CAM and TGM are almost identical. These two methods are considered to be grid independent, and the deviations between different grid resolutions are less than 3%. However, the deviations with DBM are 8.1%, 13.5% and 15.0% for Fine, Coarse and Coarser grids respectively. Eq. (14) is solved with an independent time step. In this simulation, each time before resolving the gas phase, Eq. (14) is solved for the source terms and the particle properties from 0 to time T within six time steps, as recommended by Sun et al. [20]. The numerical diffusion of the solution of Eq. (14) is dependent on both the time steps and grid size.
Since the different grids as well as the different coupling methods predict different conversion rates, the results obtained at the particle residual mass ratio of 50% are compared. The field of the gas phase volume fraction, , with the fine grid are shown in Fig. 9 (different color scales are used to show the figures more clearly). The PCM predicts a very sharp change of the fields at the particle's location, while for the other coupling methods, the fields are almost unity. The fields also show how the source terms are distributed in space when coupled to the gas phase equations. The temperature profiles along the radial direction of the combustor at the particle's location are presented in Fig. 10 at the same residual mass ratio. In PCM, all the source terms are coupled with one grid, resulting in a large temperature gradient in the gas phase. The gradients get larger as the grid size becomes smaller. This is due to that the source term in Eq. (3), which contains the enthalpy of the released gases, is distributed over different sizes in space by using different grids. The released gases, which are in heat balance with the particle, have a lower temperature than the gas phase, thus cooling down the gas phase cell. The smaller grid size means the source terms are returned to a narrower region, resulting in a sharpened gradient. The source terms should be coupled to the particle surface region, which is independent from the grid size. Meanwhile, a relatively large source term will reduce the robustness of the solver. How to control the distribution of source terms, must be considered empirically.
For CAM and TGM, the temperature gradients under different grid systems are quite similar to each other. The main reason why these two methods show better grid-independence performance is that the coupling regions are determined by x c instead of the fluid grid size. Meanwhile, the cooling effects on the particle nearby gas is smoothed, and so is the heat release from the homogeneous reactions. This makes it hard to evaluate whether the averaged gas temperature obtained from CAM and TGM is overestimated or underestimated. There is an asymmetry in the profile of the temperature for the TGM method with the Coarse grid. This is because the particle grid is placed asymmetrically with respect to the particle due to the limitation of the fluid grid. For the DBM, the influence of the particle to the gas is even more smoothed. The temperature profiles with the DBM are slightly higher. The reason could be that the conversion rate predicted by DBM method is slightly higher. The prediction of the oxygen concentration has similar issues. For the PCM, the large mass source terms of the released devolatilization gases, which consume O 2 rapidly, result in a large gradient in the O 2 concentration. As shown in Fig. 11, the high concentration of combustible gases in the central cell leads to an underestimation of the oxygen concentration, resulting of an underestimation of char conversion rate during devolatilization. It explains why the PCM predicts a high residual mass ratio at the end of devolatilization. For the DBM, the O 2 concentration is almost uniform and the char oxidization is not limited by the devolatilization, which leads to a low residual mass ratio when the devolatilization ends.
In general, when the particle size is close to the grid size, the grid size has significant influence on the PCM, while the CAM and TGM show good independence of the grid size. DBM reduces the grid dependence, but as for the method itself, the numerical diffusion of the solution of Eq. (14) is grid dependent.

Effects of the coupling parameter
Although the CAM, TGM and DBM improve grid independence, these methods introduce a new parameter, the coupling length x c which needs to be determined. In fact, this parameter has impact on the results in a similar way as the grid size in the PCM simulations. Sensitivity studies of this parameter have been conducted using the Fine grid. It is worth noting that only the initial value of x d / c p can be set for the TGM and DBM. As the particle shrinks, x d / c p becomes larger, since for these two methods the coupling length x c is a constant value. For CAM on the other hand, the ratio of x d / c p is kept constant. x c has a minimum value of one grid cell, because when x c is smaller than one grid cell, CAM and PCM become equivalent.
The predicted devolatilization time and total burnout time with the different coupling methods are compared against experimental data in Fig. 12. Devolatilization is primarily a heat transfer controlled process. All three coupling methods predict similar devolatilization times. As x c increases, the results converge. However, the numerical convergence does not necessarily imply that the results are physically correct. Due to the high temperature at the particle surface, char oxidization in the simulated case is a diffusion-controlled process. This implies that the O 2 concentration is the dominating factor in determining the char conversion rate. In fact, the char oxidization competes for O 2 with the homogeneous gas phase reactions. When the coupling region is too large, the O 2 consuming region will also be enlarged and the O 2 concentration gradient will be smoothed. As the x d / c p increases, the O 2 concentration becomes closer to the far-field Fig. 10. Gas phase temperature profile at particle location and particle surface temperature with residual mass ratio of 50%. The solid lines are gas phase temperature, and the circle markers are the particle surface temperature.
J. Zhang, et al. Chemical Engineering Journal 387 (2020) 123964 condition (Appendix B.2, Fig. 17). The larger x c overpredicts the local O 2 concentration, and correspondingly predicts a higher char oxidization rate during devolatilization. The overpredicted char oxidization results in a higher particle temperature, which also causes overprediction of the devolatilization rate. We therefore argue that the coupling length should be interpreted as an additional model parameter for a coupled reactive Eulerian-Lagrangian framework, and has a non-trivial impact on the results obtained. The gas phase temperature history at the particle location is shown in Fig. 13. The first peak at the beginning of the conversion is due to the dry wood accumulated from the drying process. The released gases from devolatilization react with oxygen raising the gas phase temperature rapidly. Afterwards, the devolatilization rate is limited by the drying rate. When x d / c p is close to 1, CAM gives the same results as Fig. 11. Gas phase O 2 profile at particle location with residual mass ratio of 50%. The different coupling lengths also have an impact on the gas phase reactions. As discussed above, for the CAM a larger coupling length results in a higher gas phase temperature during devolatilization and lower gas phase temperature during char burnout. The CO/CO 2 ratio at the outlet is shown in Fig. 15. For the char burnout, according to the kinetics used in the char oxidization, the CO/CO 2 ratio increases with the increase of temperature. The effects of the coupling length on different coupling methods are summarized in Table 6.
It is worth noting that if the coupling methods are applied to multiparticle systems with high levels of turbulence, for example fluidized bed combustors, the above conclusions may not be valid. The Nu and Sh correlations need to consider the turbulence and bed voidage [41], and the estimated from Eq. (17) should be thinner. Such industrial-scale systems form meso-scale structures, such as particle clusters, which is between the particle scale and the system scale. The x c should be smaller than the meso-scale, otherwise the coupling methods will oversmooth the gas-solid interactions. The coarse-graining methods are likely to increase the numerical stability [18], but could also break the meso-scale structures. x c need to be carefully studied, and the validation of x c will be rather empirical.

Computational efficiency
The computational cost of the different coupling methods is not only influenced by the model parameter settings, namely the choice of x c , but also related to the grids resolution. The computational cost for the particle part using PCM is taken as the reference, which means that the   15. Gas phase products distribution at outlet in CAM. A CO/CO 2 ratio of 1 means the char is totally burned out.

Table 6
Effects of coupling methods on the particle conversion.
computational time for the gas phase governing equations and homogeneous reactions are excluded. This is because adjustable time steps are used in these two parts of the calculations, and usually the time consumed by resolving homogeneous reactions is dominating in the whole simulation. The increase in computational cost associated with the different coupling methods are shown in Fig. 16.
For all the three coupling methods, it is shown that the computational time increases with the increasing of the grid number. For CAM and DBM, the computational time increasing is much faster than that of TGM. CAM requires grid searching in every fluid time step, while for TGM similar searching only required once at beginning, then the mapping relation between the particle grid and fluid grid is registered as a constant table. It is worth noticing that Fig. 16 only shows results of the single particle simulations. If more particles are added into the system, it is expected that the computational time increases for TGM and DBM are the same as the single particle simulation, because these two methods already applied to the whole calculation domain and are independent from the particle number. However, for CAM, the additional computational time required is multiplied by the particle number. In addition, even though DBM requires a relatively large computational cost for the fine grid simulation, it is still of the same order as the cost for resolving the governing equations of the gas phase, which means this method is feasible. Considering the computational efficiency is a practical issue for multi-particle simulation, TGM is the most efficient coupling method, while CAM could be the most expensive one.

Conclusion
CFD simulations of combustion of a solid fuel particle are conducted using an Eulerian-Lagrangian approach, employing different coupling methods between the reacting gas and solid phases. The three coupling methods, CAM, TGM and DBM have all been extended for reacting particles, and are able to improve the grid independence of the CFD solver. When linking the single particle model and the Eulerian gas phase model, the interaction between the particle and the gas phase is shown to occur within a certain coupling length scale. When the particle size is comparable to the grid size, the coupling length becomes a critical model parameter. The coupling length affects the boundary condition of the particle sub-model, which is sampled from the resolved gas phase model, and also influences the distribution of the solid phase properties and interaction source terms. The results show that for all the three methods, a small coupling length underestimates the heating and devolatilization processes, while a large coupling length overestimates O 2 concentration and weakens the influence of the gas phase reactions. The coupling length can be evaluated by the estimation of the boundary layer thickness using given correlations. In this study, the coupling length of 3d p is shown to be a reasonable estimation. All the coupling methods introduce additional diffusion of the source terms, which further influences the gas phase profiles and products distribution. The source terms should be limited to the region near the particles' surfaces, but this aspect of the implementation is subject to the robustness of the solver in question. The computational efficiency of the three methods is also evaluated. TGM is believed to be the most efficient method for potential multi-particle simulations.

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.
T g is the local gas temperature, which must be obtained through the coupling scheme. Another boundary condition is the heat flux to the particle centroid, which is represented by the wet wood layer temperature. The heat flux is modeled by a thermal drying model, where the drying rate is determined by the evaporation heat transferred from the dry wood layer to the drying front. The mass transfer process of vapor diffusing out of the particle is considered by a correlation function [5]. The maximum temperature of the drying front is limited by the water boiling point at the given gas phase pressure. The mass balance for each layer is calculated by the reaction rate on the boundary. The ash is assumed to be an inert component in every layer and will be transferred to the outer layer according to the mass loss of each layer. The particle ash content needs to be pre-defined. Similar to Thunman's method [22], a shrinkage model using empirical shrinkage factors, i , is employed to update the particle volume according to the mass changes. The volume change of the ith layer is calculated by: where m pi is the mass consumed on the ith boundary.