The present study is devoted to a two-phase, gas–solid flow in an annular pipe (hollow cylinder) at an elevated pressure of 15 bars and moderate Reynolds number of circa 6000. The influence of the particle loading, the interaction between the phases, and turbulence dispersion on the flow dynamics is systematically studied by means of computational fluid dynamics simulations, employing the Ansys FLUENT commercial package. The cases with a particle volumetric fraction of 1.2% are referred to as “high particle loading,” and those with 0.13% are denoted as “low particle loading.” The following cases are investigated: (1) pure gas flow; (2) low particle loading two-phase flow with one-way coupling and with turbulence dispersion; (3) low particle loading two-phase flow with two-way coupling but without turbulence dispersion; (4) low particle loading two-phase flow with two-way coupling and with turbulence dispersion; (5) high particle loading two-phase flow with one-way coupling and with turbulence dispersion; (6) high particle loading two-phase flow with two-way coupling but without turbulence dispersion; and (7) high particle loading two-phase flow with two-way coupling and with turbulence dispersion. The boundary layer is found to grow without fluctuations of the turbulent kinetic energy (TKE) for cases 1, 2, and 5. For case 4, the TKE fluctuations have been identified, although they appear to be less substantial than those in cases 6 and 7. The authors attribute the semi-chaotic nature of the TKE fluctuations to the particle loading and two-way coupling. In addition, the onset and development of the flow instability have been observed at a random axial distance in cases 4, 6, and 7. Such instability is also attributed to the two-way coupling with turbulence dispersion in the flow. It is concluded that the particle loading, one-way, or two-way coupling between the phases, and the turbulence dispersion models significantly influence the development of the flow dynamics with the same inlet and boundary conditions. Consequently, it is not a trivial question, which result a user should trust. The present computational results inspire to perform verification as well as experimental validation of the simulations, so the simulation results can subsequently be used with confidence for design analysis.

## I. INTRODUCTION

Ubiquitous in nature and engineering applications,^{1} multiphase flows, especially two-phase flows, are categorized as (i) *dispersed* and (ii) *separated*.^{2,3} Two-phase flows are further classified into three distinct groups: (i) *gas–solid*, (ii) *gas–liquid*, and (iii) *liquid–solid* flows.^{4} Since modeling of turbulent two-phase, gas–solid flows still faces challenges, comprehensive and fundamental studies of such flows are important as they can help to enhance the effectiveness of the industrial processes, for example, the overall efficiency of the coal combustion systems.

While turbulence is a complicated phenomenon, one of the basic traits of a turbulent flow is the presence of random and chaotic fluctuations, with a wide spectrum of length and time scales.^{5} Enhanced diffusion that prompts mass, momentum, and energy transfer makes turbulence a useful phenomenon in most engineering applications;^{6} however, it may need to be avoided in some cases, for example, in order to reduce the particle deposition on the reactor walls.

Giacinto *et al.*^{7} demonstrated that particle loading (volumetric fraction of particles) is one of the primary parameters determining the selection between the approaches of one-way or two-way coupling for a computational simulation of a two-phase flow. In the parametric study of a gas-particle flow in a vertical duct, Liu and Glasser^{8} identified a segregation of the particles carried toward the wall from the center, with the increase in the solid phase loading. In Ref. 8, the particles were considered as perfectly elastic. However, when inelastic, the particles were found to be segregated toward the center of the channel. Sinclair and Jackson^{9} investigated a gas-particle flow in a vertical pipe accounting for the particle–particle interaction and observed a segregation of both phases in the radial direction. Likewise, Kartushinsky *et al.*^{10} performed a similar numerical study of a three-dimensional (3D) gas–solid flow but in a horizontal pipe. Here, an asymmetric flow field was observed, with the asymmetry being attributed to gravity. In a semiempirical model for an upward vertical two-phase flow system, Rhodes *et al.*^{11} implied that the radial concentration profile of the solids in the flow is governed by the cross section solid concentration. Ding *et al.*^{12} observed a laminar-to-turbulent transition at a Reynolds number of $Re=620$, which is far below the conventional $Re$ of circa 2300.

As the flow configuration without a swirl, employed in the present work, is not commonly used in industrial applications, our literature overview resulted in a limited number of relevant studies. In particular, most of the literature on two-phase flows in annular pipes is devoted to a gas–liquid flow in the so-called “core-annular” flow regime. Other fundamental studies are mostly devoted to a flow between concentric cylinders, with one or both cylinders being rotating. Unlike geometry of the present work, a relevant study by Hajidavalloo *et al.*^{13} was devoted to a gas–solid flow in a vertical annulus, with the flow direction opposing gravity. Moreover, Ref. 13 employed the inner cylinder off-center, with a rotating drill bit at the entrance of the annular region, although the authors of Ref. 13 mentioned that the effect of rotation was negligible. The main finding of Ref. 13 was that the fluid (air) velocity profiles exhibit two peaks near the inner and outer walls in the presence of a significant particle loading. In the present study, this feature is observed when a two-way coupling is used (cases 4 and 7). However, Hajidavallo *et al.*^{13} did not show any experimental evidence for this phenomenon. Since Ref. 13 and the present work employ the same model (and the same, Ansys FLUENT software), our study also questions the validity of the model. Earlier studies such as those by Supan and Adewumi^{14} and by Sharma and Chowdhry^{15} dealt with a vertical gas–solid flow opposing gravity, and they focus primarily on a pressure drop rather than the details of the turbulent flow dynamics. In this light, the configuration employed in our work constitutes a rather unique application and deserves a careful study.

Although numerous studies have shown the influence of multiple model parameters in the gas–solid phase flow simulations, systematic studies considering one parameter at a time are lacking. The present work therefore encompasses a systematic study of the influence of the three key factors—the particle volumetric fraction (i.e., particle loading); phase coupling; and turbulence dispersion—on the gas–solid, two-phase annular flow with a Reynolds number of about 6000 inside the inlet section of a burner of a pressurized coal combustor. The solid phase is therefore the coal particles.

## II. METHODOLOGY

### A. Turbulence modeling

While various computational tools or methods are available to simulate a turbulent flow, here we employ the Reynolds-averaged Navier–Stokes (RANS) simulations with a commonly used version of the $k$-$\omega $ model.

#### 1. RANS simulation

Solving most of the fluid flow problems by means of RANS starts with the flow analysis, focusing on the time-averaged velocity field. Substituting the Reynolds decomposition $u=u\xaf+u\u2032$ into the continuity and momentum equations, followed by the ensemble averaging, yields the following equations:

where $u\u0303i,\u2009u\u0303j,\u2009and\u2002u\u0303k$ are the three Favre-averaged velocity components (mass-weighted averaged values for compressible flows), $ui\u2032\u2009and\u2009uj\u2032$ are the fluctuating velocities, and $\mu $ is the dynamic viscosity of the fluid, with $\rho \xaf\u2009and\u2009p\xaf$ being the time-averaged density and pressure, respectively,^{16} coupled to the temperature $T$ via the ideal gas law, $p=\rho RT$, where $R$ is the specific gas constant. The applicability of the ideal gas model for a pressure of 15 bars, employed in this work, is certified by the classical thermodynamic criterion that the compressibility factor $z\u2261p/\rho RT$ should be within the range $0.95\u2264z\u22641$. In our case, according to the Nelson–Obert compressibility chart, the compressibility factor of the carrier phase is $z\u2009\u2248\u20090.98$, which justifies the ideal gas approximation. The term $\xaf\rho ui\u2032uj\u2032\xaf$ in Eq. (3) is known as the Reynolds stress, proper modeling of which tackles the issue of the closure problem—the cornerstone behind the RANS simulations.^{17} The closure problem is numerically solved by employing the Boussinesq hypothesis, which is postulated based on the observation that the mixing caused by the large energetic turbulence eddies dominates over the momentum transfer in a turbulence flow.^{18} Analogous to a laminar flow, the dependence of the turbulence shear stress on the mean strain rate is assumed to be linear and obeying the formula

where the factor $\mu t$ is the eddy viscosity, and $k$ is the specific turbulent kinetic energy (TKE),

where $u\u2032,\u2009v\u2032$, and $w\u2032$ are the fluctuating velocity components in the *x*-, *y*-, and *z*-directions, respectively. Based on the additional transport equations, there exist various RANS models, among which we chose the $k$-$\omega $ shear stress transport (SST) model for the present work. For the $k$-$\omega $ model—an extension of a much simpler $k$-ϵ model—the eddy viscosity is defined as

where $\omega =\u03f5/k$ is the specific dissipation rate, while the constants *s*, F_{2}, and $\alpha *$ are the strain rate, the blending function, and the damping coefficient of the turbulent viscosity, respectively. The quantities $k$ and $\omega $ are obtained by solving the following equations:

Here, $G\u0303k$ and $G\u0303\omega $ represent the generation of the TKE and the specific dissipation, respectively. Likewise, $\Gamma k$ and $\Gamma \omega $ show the effective diffusivity, while $Yk$ and $Y\omega $ stand for the dissipation of the TKE and the specific dissipation, whereas $D\omega $ represents the cross diffusion, while $\sigma k$ and $\sigma \omega $ are the turbulence Prandtl numbers for *k* and $\omega $, respectively.^{16}

### B. Numerical modeling of a gas–solid flow

Mainly, there are two methods of modeling gas–solid flows, namely: (i) Eulerian–Eulerian (EE) and (ii) Eulerian–Lagrangian (EL) approaches.^{19} The latter, EL approach is employed in this study as it provides more latitude in the analysis of a flow with a wide range of particle sizes, types, and shapes as compared to the EE methodology.^{19}

#### 1. The Eulerian–Lagrangian (EL) approach

In the EL approach, the fluid/gas phase is treated as a continuum medium, while the particles are tracked using the Lagrangian approach, according to which they are considered as a discrete phase. The momentum equation governing the gas phase is

where $F\u2192$ stands for an external force emerged due to interactions with the discrete phase. Likewise, the momentum equation for the particles is based on the Newton's 2nd law of motion in the form

where $Fd\u2192,Fb\u2192,F\u2192VM\u2009,F\u2192fp,and\u2009F\u2192PP$ represent the drag force, buoyancy force, virtual mass force, fluid–particle interaction, and particle–particle interaction. In the present work, the *discrete-phase model* (DPM) is adopted. The volumetric fraction of the particles plays a pivotal role in the modeling of a particle-laden flow,^{20} and for a gas-particle flow with a volumetric fraction of particles less than 10%, the DPM usually provides a good prediction of a flow field.^{21} For such a low volume fraction, the effects of particle pressure and viscous stresses due to particles, i.e., particle–particle interactions, are not considered.^{22} Consequently, the dominating item in the Lagrangian frame is the drag force defined as

where $u\u2192g\u2009and\u2009uP\u2192$ stand for the gas and particle velocities, with the factor $FD$ given by

with the Reynolds number

and the drag coefficient $CD$ given by the formula

where

and the shape factor $\xd8\u2261s/S$ is defined as the ratio of the surface area of a sphere having the same volume as the particle, $s$, to the actual surface area of the particle, *S*. In this study, the spherical drag law is employed such that the shape factor is waved.

### C. Coupling between the phases

In a gas–solid flow, drag and turbulence of the gas or the carrier phase always influence the solid or particle phase. However, only at relatively high particle loading, the carrier phase is influenced by the particles through the modulations in momentum and turbulence. The DPM available in the Ansys FLUENT software package provides both the one-way and two-way coupling options.^{22} A flow, in which one phase influences the other phases, is defined as *one-way coupling*, while the presence of mutual influence of phases on each other is regarded as *two-way coupling*. The coupling of the phases can be performed via mass, momentum, and energy transfer.^{22} Within the DPM approach, the mass, momentum, and energy gained or lost by the particles are calculated while computing the trajectory of the particles. The change of these quantities is incorporated into the continuous phase calculations for two-way coupling,^{23} and Eqs. (12) and (13) are solved for gas and particle phases, respectively, until the solutions for both phases converge.^{24}

### D. Turbulence dispersion

The dispersion of the particles due to the turbulent component of the fluid phase velocity is known as *turbulence* dispersion,^{25} which can be modeled for particles in a gas–solid flow by either

the particle cloud model or

the stochastic model/the discrete random walk (DRW) model.

In the stochastic model, which is the choice for the present simulations, the instantaneous equation for the particle motion is solved.^{26} When turbulence dispersion is not included, Eq. (13) is solved using only the mean velocity component $u\xaf$. However, to account for turbulent dispersion, the fluctuating velocity component $u\u2032$ in all three directions is also included; hence, the particle motion equation is solved using the instantaneous gas velocity. The particles are assumed to interact with the turbulent eddies in continuous phase, which are characterized by the normally distributed velocity fluctuations. Specifically, the random fluctuating velocity is given by

where $\zeta $ is a normally distributed random variable. For a particular turbulence model, the TKE ($k)$ is calculated and, hence, the root mean square (RMS) velocity is found as

In most cases, isotropic turbulence is assumed due to the lack of detailed information about local fluid turbulence.

### E. Geometry and domain discretization

The present work is devoted to a fundamental behavior of a gas–solid, two-phase flow in a hollow cylinder with an annular geometry and a uniform cross section area. This study is conducted to understand the effect of the key model assumptions and parameters on the flow dynamics inside the inlet section of a burner of a coal combustor. Figure 1 represents the wireframe and the solid structure of the configuration considered, with the inlet section, the outer wall, the inner wall, and the outlet section depicted by the blue, green, red, and yellow regions, respectively. The details of various geometrical dimensions are listed in Table I.

### F. Solution methods and boundary conditions

Transient RANS simulations were performed for a pure gaseous flow as well as for a gas–solid phase flow. The $k$-$\omega $ SST model has been selected to better capture the flow phenomena near the wall and away from it. This model has been validated by Tian *et al.*,^{27,28} who conducted a thorough numerical investigation of coal combustion in a utility furnace and gas-particle rectangular jets in a crossflow. They have compared their numerical results to the experimental data, with the $k$-$\omega $ SST model providing best agreement with the experimental data.

The present study targets an understanding of the flow conditions to which the burner of a coal combustor is exposed. Let us briefly describe the setup: the burner is the inlet section of the combustor, where the flow is cold and non-reacting. Specifically, the section considered in the simulations constitutes the inlet section of the actual combustor, i.e., this section is used to deliver/convey the reactants to the reaction chamber (reactor), and the reactor starts at the end of the section being simulated here. Consequently, all properties of the solid particles are chosen to be those of coal. At the same time, we believe that our numerical results can be applicable for other types of particles as long as the Stokes number for such particles is on the same order of magnitude as the Stokes number for the coal particles considered here. A more rigorous limitation is that the EL approach will not work if the particles volume fraction exceeds 10% (see above), but in the present work we are very far from such a limit (see below). Also, various studies^{29–32} have shown that the coal plants operating at higher pressures have higher overall efficiencies. Therefore, an operating pressure of the domain is set to 15 bars as that for a whole combustor, which is being experimentally investigated at Washington University in St. Louis (WUSTL).^{29,30} The geometry employed in our simulations resembles that in the WUSTL experimental facility, but in the ongoing WUSTL experiments, the diameter of the inner cylinder grows gradually, starting from setting an axial point to accelerate the flow in the annular region, before the reactor is reached. Unfortunately, measurements in such a pressurized system face difficulties such that we do not have the experimental data to validate our work. In this light, the present simulations are considered as predictions, the results of which are being used as a guidance to improve the experimental design.

The simulation has been started with a pure gaseous flow, followed by injection of the particles after sufficient simulation time passed, namely, circa five gas-flow through times (FTT). Here, the FTT is defined as the ratio of the volume of the computational domain to the volumetric flow rate of the gaseous phase. The gas–solid, two-phase flow dynamics is modeled using the EL approach, with the coupling between the phases devoted to the source term in the governing equation. Particularly, the DPM was selected for this purpose. The particles of a uniform size of 65 *μ*m, which mimics a typical average particle size in the ongoing experimental work, were injected uniformly through the inlet cross section. The DRW model was selected for the simulation runs accounting for turbulence dispersion of the particles. The spatial discretization for the scaler parameters was implemented using both the 2nd-order upwind and 3rd-order monotone upstream-centered schemes for the conservation laws (MUSCL). It is, however, noted that no notable difference was observed between these two schemes. The transient analysis was performed using the bounded 2nd-order implicit method. The vector parameters have been spatially discretized using the 2nd-order upwind method. As for pressure-velocity coupling, the semi-implicit method for the pressure-linked equation (SIMPLE) was used. Table II lists the details of the boundary conditions.

Parameters . | Values . |
---|---|

Inlet mass flow rate of air | 2.45 g/s, uniform gas velocity u, v = 0, w = 0 |

Turbulent intensity at the inlet | 5% |

TKE | 2.365 × 10^{−4} m^{2}/s^{2} |

Specific dissipation rate (ω) | 88.487 66 1/s |

Hydraulic diameter | 0.026 806 m |

Inlet mass flow rates of particles | 3.98 and 0.39 g/s |

Size of injected particles | 65 μm, uniformly distribution over the inlet cross section |

Wall boundary condition | No-slip for gas, perfect reflection for particles |

Parameters . | Values . |
---|---|

Inlet mass flow rate of air | 2.45 g/s, uniform gas velocity u, v = 0, w = 0 |

Turbulent intensity at the inlet | 5% |

TKE | 2.365 × 10^{−4} m^{2}/s^{2} |

Specific dissipation rate (ω) | 88.487 66 1/s |

Hydraulic diameter | 0.026 806 m |

Inlet mass flow rates of particles | 3.98 and 0.39 g/s |

Size of injected particles | 65 μm, uniformly distribution over the inlet cross section |

Wall boundary condition | No-slip for gas, perfect reflection for particles |

## III. RESULTS AND DISCUSSIONS

Based on the volumetric fraction of particles in the flow, two different particle loading conditions have been created and are denoted as ** high** and

**particle loading—with the volumetric fractions of 1.2% and 0.13%, respectively. It should be noted, in this respect, that while there is no rigorous mathematical threshold between the high and low particle loading, by an order of magnitude, we generally consider the particle loading as high if the particles volume fraction substantially exceeds 0.1%, while it is considered as low otherwise. Nevertheless, 0.1% is not a cutoff value but rather a guide in the field of two-phase, gas–solid flows for a transition at which the two-way coupling—i.e., the influence of the particles (the discrete phase) on the gas (the continuous phase)—cannot be neglected anymore. To be specific, it is accepted in the literature that beyond such a transitional point of 0.1%, one-way coupling is not fully valid, i.e., the impact of the particles on the continuous phase starts becoming important. Regarding the particle loading quantities employed in the present work, 1.2% is the actual loading targeted in the WUSTL experiments. We therefore treat it as high particle loading. In addition to this quantity, we employed that reduced by an order of magnitude (0.13%) in order to come close to the limit of 0.1%. These two particle loadings are compared aiming to identify whether the high particle volume fraction might cause a significant gradient of the DPM concentration observed in the boundary layer. It is also noted that both quantities of 1.2% and 0.13% are well below the limit of 10%, thereby making the EL approach viable. Overall, in the present work, seven cases have been distinguished and investigated by the combinations of phase coupling and turbulence dispersion. Table III presents details of the various cases considered.**

*low*Cases . | Phases . | Particle loading . | Coupling . | Turbulence dispersion . |
---|---|---|---|---|

1 | Pure gas phase | No particle | No coupling | Off |

2 | Two-phase | Low | One-way | On |

3 | Two-phase | Low | Two-way | Off |

4 | Two-phase | Low | Two-way | On |

5 | Two-phase | High | One-way | On |

6 | Two-phase | High | Two-way | Off |

7 | Two-phase | High | Two-way | On |

Cases . | Phases . | Particle loading . | Coupling . | Turbulence dispersion . |
---|---|---|---|---|

1 | Pure gas phase | No particle | No coupling | Off |

2 | Two-phase | Low | One-way | On |

3 | Two-phase | Low | Two-way | Off |

4 | Two-phase | Low | Two-way | On |

5 | Two-phase | High | One-way | On |

6 | Two-phase | High | Two-way | Off |

7 | Two-phase | High | Two-way | On |

### A. Case 1: Pure gas flow

First, the simulation was initiated with a pure gaseous (air) flow, constituting a purely single-phase problem. The axial velocity at four different locations was plotted, namely, *x* = −0.2 m (the inlet), *x* = −0.1 m, *x* = −0.05 m, and *x *=* *0 (the outlet). The locations for *x* = −0.1 m and *x* = −0.05 m are shown by the two vertical dark lines in Fig. 2.

A horizontal line, going from the inlet to the outlet, inside the boundary layer and in the vicinity of the wall, at *r *=* *0.014 m, was created in order to capture the fluctuations inside the boundary layer. The TKE and the axial velocity along the horizontal line were monitored to observe the flow field fluctuations. The plots for the TKE and the axial velocity along the horizontal line are shown in Figs. 3 and 4, respectively.

Initially, the gaseous phase is initiated with a uniform axial velocity at the inlet. The results show the expected developing flow in the annular region of concentric cylinders. First, the TKE grows. Subsequently, after a certain distance, it diminishes and saturates as expected. No fluctuations of the velocity and the TKE profiles are observed for this case.

The axial velocity, plotted in Fig. 5 vs the radial distance at various axial locations, shows the developed flow field. However, again, no sign of fluctuations is seen. The axial velocity of the gas phase at the center grows, indicating that the thickness of the boundary layer along both walls is growing, thus showing an axially developing flow. There is no indication of any flow instability.

### B. Verification of mass conservation

After the particles were injected into the gaseous flow, the contours of the DPM concentration were inspected. Beyond expectations, significant concentration gradients and clustering are seen in Fig. 6. A series of simulations was initiated to understand such peculiar flow dynamics. It is noted that the results shown in Fig. 6 are from case 7 (see Table III). First, the conservation of the particles mass and particles number flow rates were checked before proceeding with other cases. The entire flow domain was treated as a control volume. Obviously, the particles entered through the inlet and exited through the outlet. Since the mass flow rate of the particles injected at the inlet was known, the mass of the particles and, hence, the number of particles entering and exiting the control volume were calculated as

Here, $M\u0307$ is a mass flow rate of the particles, $M$ is a mass of the total number of particles entering the control volume at each injection time step, and $Np$ is the number of particles entering through the inlet during each injection. Each particle has a volume $V\u2212$ and a density $\rho p$. The mass flow rate of the particles at the inlet was found to be constant; however, it was fluctuating at the outlet. The difference between the number of particles entering and exiting the domain is shown in Fig. 7. Since such a difference (random fluctuations according to the DRW model) is relatively small, and its time average is zero, this is taken as a proof of the conservation of the particles number and mass. It should be noted that the particles used in the simulations are of constant size. The DPM concentration gradient was also checked by injecting the particles with a Rosin–Rammler size distribution in the range from 5 to 200 *μ*m. However, those DPM concentration gradients did not notably differ from the cases presented in this paper. Hence, for the purpose of simplicity and brevity of the paper, those cases are not included in the present work. Such an analyse of the role of various particle sizes will be presented elsewhere.

### C. Impact of gravity in a flow

The axial velocities, with which the gas and the solid particles were injected at the inlet boundary, were the same (about 0.25 m/s). However, immediately after the inlet, a difference in the particle velocity and the gas velocity was observed, circa 0.02 m/s, being close to the terminal velocity of a 65 *μ*m particle, as expected. The particle and gas velocities at the outlet are shown in Fig. 8. It is seen that the change in the gas velocity as compared to the inlet velocity of 0.25 m/s is not significant. However, there is a change in the particle velocity due to gravity. To address this point, an additional simulation run was performed to see if a clustering of particles would persist without gravity; the DPM-concentration contours of Fig. 9 certify this, albeit with more coherent flow structures in the near-wall regions.

### D. Case 2: Low particle loading with turbulence dispersion and one-way coupling

In this case, the same process as in case 1 was repeated. The TKE and the axial velocity of the gas, plotted along the boundary layer, were inspected. Also, the axial velocities at the various axial locations were studied. The flows did not differ from the pure gaseous flow case, indicating that the low particle loading with one-way coupling did not influence the gaseous flow field, as expected.

### E. Case 3: Low particle loading with no turbulence dispersion and two-way coupling

For this case, the axial velocity and TKE were also monitored at various locations and times. While the gas-phase velocity profiles at different cross sections were found to be the same as in the previous cases, the respective particle concentration contours are depicted in Fig. 10. The results from this case did not differ much from that of the previous one (case 2). Despite the use of two-way coupling, no instability or fluctuations in the flow field of the gas phase were observed. Consequently, it was suspected that such phenomena appearing in the flow should be the result of turbulence dispersion and/or high particle loading. Extra cases were set and simulated to validate this hypothesis.

### F. Case 4: Low particle loading with turbulence dispersion and two-way coupling

This case was modeled using two-way coupling with the incorporation of the turbulence dispersion of the particles. The axial velocity and the TKE profiles differed from those obtained in the earlier cases. Specifically, the DPM-contours in Fig. 11 showed the clustering of the particles at different locations, appearing near the walls when two-way coupling is imposed. This clearly shows that two-way coupling (as implemented in the Ansys FLUENT) significantly modifies the flow field (despite a low particle loading with a volumetric fraction of circa 0.13%), when disturbed by the small turbulent fluctuations arising from the DRW model. The TKE, Fig. 12, and the axial velocity, Fig. 13, plots along the boundary layer show an instability in the flow field in the form of coherent oscillations. In this respect, the axial velocity profiles of the gas-phase at various locations exhibit relatively large distortion/oscillations with time; see Fig. 14. Hence, the results obtained from all three low particle loading cases show that turbulence dispersion of the particles, along with the two-way coupling model, triggers an instability in the flow field. Thus, turbulence dispersion of the solid particles can impact the gas-phase when the terms describing two-phase coupling are included in the gas-phase equations, despite a relatively low particle loading.

### G. Case 5: High particle loading with turbulence dispersion and one-way coupling

This case presents the results obtained for the high particle loading with the effect of turbulence dispersion of the particles, but without the influence of discrete particles on the continuous phase. For this case, the TKE profile and the axial velocity profiles resembled those for case 2. This again shows that two-way coupling and turbulent dispersion are the two primary factors that contribute to instability in the flow dynamics by triggering a transition-like behavior in the boundary layer. Here, we use the term “transition-like” because the flow is already treated as turbulent, with an inlet turbulence intensity being circa 10% of the inlet flow velocity.

### H. Case 6: High particle loading with no turbulence dispersion and two-way coupling

Unlike case 5, although case 6 underlines the two-way interaction between the phases, turbulence dispersion of the particles was not considered. All the parameters for this case are identical to case 3, except for the particles mass flow rate. In this case, two-way coupling causes the boundary layer to become unstable, with a transition-like phenomenon observed around *x* = −0.1 m (see Fig. 16). Also, Figs. 15–18 delineate the nature of the large-scale flow structures caused by the particle clustering, plausibly, in the low vorticity regions inside the boundary layer. The oscillations are excited more due to the high particle loading. Indeed, the mass flow rate of particles for this case exceeds that for the low particle loading by a factor of an order of 10.

### I. Case 7: High particle loading with turbulence dispersion and two-way coupling

Finally, turbulence dispersion of particles is considered along with two-way coupling in the case of high particle loading. The DPM concentration contour plot, Fig. 6, showed particle clustering at various locations near the wall, this time, intenser than that in the other cases. The corresponding TKE and the axial velocity plots (Figs. 19 and 20, respectively) show the large fluctuations along the boundary layer. Likewise, the velocity profile delineated in Fig. 21 show an asymmetry with respect to the centerline of the tube. Hence, activating turbulence dispersion of particles in cases of high particle loading enhances and sustains the flow field instability, thereby causing an earlier transition at circa *x* = −0.16 m, similar to case 4. Another phenomenon occurring at high particle loading is an appearance of two peaks in the gas velocity near the walls, as opposed to one peak near the duct center for low particle loading; compare Fig. 5 with Figs. 18 and 21.

### J. Influence of mesh and time step

In order to assess the snesitivity of the present results to the grid size, a mesh sensitivity study was performed for case 3. Since this is not a primary objective of the current sudy, we present the results of the grid-sensitivity analysis in the Appendix, which confirms that the grid resolution used in our analysis is adequate for the purpose of the present work.

## IV. CONCLUSIONS

With an objective to understand the fundamental physics governing a gas–solid, two-phase flow inside a vertical annular pipe, representing the inlet section for an envisioned low turbulence coal-combustion reactor without swirl, seven different cases have been investigated. The roles of phase coupling, particle loading, and turbulence dispersion were studied.

For the cases with high and low particle loading, one-way coupling, and turbulence dispersion (cases 2 and 5), the presence of the particles in the flow field did not have a notable influence on the gas-phase, as expected. It is found that one-way coupling between the phases, by intention, does not show any significant influence on the solid phase, even with a relatively large particle loading of 1.2% by volume. However, the two-way coupling model (even when used with a low particle loading of 0.13%) and the turbulent dispersion model did show a significant impact on the flow dynamics by triggering a transition-like flow instability in the boundary layer and clustering of the particles near the walls.

In addition, for a flow with a two-way coupling but no turbulence dispersion (cases 3 and 6), different flow behaviors were observed. For a low particle loading, case 3, the flow field seemed to have no flow instability, while case 6, with a high particle loading, showed pseudoperiodic coherent flow structures in the flow field. High particle loading with two-way coupling was found to create the particle clustering in the latter case and was attributed to this difference.

Finally, for cases with a two-way coupling flow and turbulence dispersion (cases 4 and 7), the particle clustering accompanied by pseudoperiodic flow structures was observed for both cases. Hence, it is concluded that the phase coupling, turbulence dispersion, and particle loading play significant roles in determining the two-phase flow dynamics inside the annular region of two circular-concentric pipes. Particularly, even in the case of low particle loading (0.13%), while one-way coupling exhibits neither flow instability nor particle clustering, the situation is completely different if two-way coupling is adopted, thereby questioning the validity of the two-way coupling model.

The results of this study call to question which of these models reflect the reality. In other words, is there particle clustering near the walls accompanied by coherent/chaotic flow structures or not for the two-phase flow investigated? If yes, then at what location does the transition occur along the boundary layers? These two different scenarios would lead to completely different combustion scenarios inside the reactor, unless the flame is stabilized by other mechanisms. Once again, the necessity of the validation process is demonstrated even when using a presumably “validated” code. It is expected that there will be some experimental data available in the near future to validate these numerical results, while we will be continuing the simulation work with more sophisticated but computationally extensive and expensive methods such as Large Eddy Simulations (LES).

## ACKNOWLEDGMENTS

This work was supported by the U.S.-China Clean Energy Research Center (CERC)—Advanced Coal Technology Consortium (ACTC) under DOE Contract No. DE-PI0000017.

### APPENDIX: A GRID SENSITIVITY STUDY

In order to investigate the sensitivity of the results to the mesh refinement, we took case 3 as an example and ran it employing two other time steps, namely, 0.002 and 0.0005 s, in addition to 0.001 s. It is recalled that case 3 is the case with a low particle loading, two-way interaction, without turbulence dispersion. Observation of the TKE, the mean gas velocity profiles, and the DPM concentration contours did not exhibit any significant difference. Based on these observations, we concluded that the time step of 0.001 s used in our simulations was sufficiently small to achieve suitable numerical accuracy. We also studied the effect of spatial mesh size by both coarsening and refining the mesh we used in most of the simulations. The meshes used are characterized by the minimum, *h*_{min}, and maximum, *h*_{max}, values (Table IV) calculated from

We chose the TKE (in m^{2}/s^{2}) and DPM concentration (in kg/m^{3}) profiles as the analysis variables, while the others did not exhibit much sensitivity to the mesh refinement. These profiles are depicted in Figs. 22–25. It is seen that there is significant grid dependence observed from the variation of both variables. The mesh we used for most of the analysis, marked as “medium,” in fact, is seen to be relatively coarse. Since the current study is not a validation study, i.e., our aim is not to assess modeling errors using measured data, we find the mesh used in our calculations sufficiently fine for the stated purpose. The radial profiles of Figs. 22 and 23 show that the boundary layer near the walls becomes thinner with the mesh refinement. The axial profiles show that a relaxation from a uniform plug-flow to a classical bounday-layer type flow occurs much faster as the mesh is refined.

Mesh . | h_{max} (m)
. | h_{min} (m)
. | h (m)
. _{avg} |
---|---|---|---|

Fine | 0.001 134 | 0.000 818 | 0.0010 |

Medium-fine (m-f) | 0.001 283 | 0.000 930 | 0.0011 |

Medium | 0.001 446 | 0.001 068 | 0.0013 |

Coarse | 0.002 181 | 0.001 616 | 0.0019 |

Mesh . | h_{max} (m)
. | h_{min} (m)
. | h (m)
. _{avg} |
---|---|---|---|

Fine | 0.001 134 | 0.000 818 | 0.0010 |

Medium-fine (m-f) | 0.001 283 | 0.000 930 | 0.0011 |

Medium | 0.001 446 | 0.001 068 | 0.0013 |

Coarse | 0.002 181 | 0.001 616 | 0.0019 |

An unusual behavior is seen in the DPM profile (Fig. 25), whereby some superficial oscillations (or wiggles) are seen in the DPM concentration profiles. The frequency of these oscillations do not change much as the mesh is refined. It is difficult to assert whether these are some instabilities in the calculations or they are some artifacts of numerics. Nevertheless, these oscillations are probably the precursors to the chaotic behavior which was observed in the DPM concentration contour plot (see, e.g., Fig. 6).

#### DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## References

*Proceedings of*

*, Lemont, IL, USA, 4–7 October 2004*(