Investigation of CO2 Photoreduction in an Annular Fluidized Bed Photoreactor by MP-PIC Simulation.

Carbon dioxide (CO2) photoreduction is a promising process for both mitigating CO2 emissions and providing chemicals and fuels. A gas-solid two-phase annular fluidized bed photoreactor (FBPR) would be preferred for this process due to its high mass-transfer rate and easy operation. However, CO2 photoreduction using the FBPR has not been widely researched to date. The Lagrangian multiphase particle-in-cell (MP-PIC) simulation with computational fluid dynamic models is a new and robust approach to explore the multiphase reaction system in the gas-solid fluidized bed. Therefore, the purpose of this paper is to investigate CO2 photoreduction in the FBPR by MP-PIC modeling to understand the intrinsic mechanism of solid flow, species mass transfer, and CO2 photoreaction. The MP-PIC models for solid flow in the FBPR were validated by the bed expansion height and bubble size. The results showed the particle stress of the Lun model, the drag of the Ergun-WenYu (Gidaspow) model, and the coefficient of restitution e = 0.95 with the wall parameters e w = 0.9 and μw = 0.6 are the best fit to the experimental empirical correlations. The MP-PIC models developed in this work proved to be better than the Eulerian two-fluid modeling in the prediction of the bed expansion height and bubble size. It was also found from the simulation results that the maximum radiation intensity is in the half reactor height area, and the photocatalytic reaction mainly occurred around the inner wall. It showed that the gas velocity and catalyst loading were two crucial operating parameters to control the process. The results reported here can provide guidance for the operation and reactor design of the CO2 photoreduction process.

Introduction

CO2 photoreduction is attracting significant interest as it has the potential to decrease CO2 emissions as well as to provide chemicals and fuels.[1] However, this process has a huge challenge in practical applications due to a very low conversion rate.

Previous research demonstrated that CO2 photoreduction is a mass-transfer-controlled process.[2,3] Therefore, strategies to increase the mass-transfer rate are sought when using gas–liquid–solid three-phase photocatalytic reactors, such as a high CO2 pressure to increase CO2 solubility in the liquid;[4] mechanical stirring; or gas bubbling.[5] Additionally, some novel composite catalysts such as Cu2O nanoplatelets supported on multilayer graphene have been developed by Hurtado et al.[6] for improving CO2 reduction due to the enhanced electron migration from the photocatalyst into the graphene. Ola and Maroto-Valer[7] and Khan and Tahir[8] summarized the photoreactors for CO2 photoreduction. They divided the photoreactors into three groups: (i) gas–liquid–solid three-phase slurry photoreactors, such as twin reactors; (ii) fixed-bed photoreactors, such as packed bed, thin-film bed, optical fiber bed, and monolith honeycomb bed reactors; and (iii) membrane photoreactors, such as slurry bed membrane and fixed-bed membrane reactors. Traditional gas–liquid–solid three-phase slurry photocatalytic reactors are not conducive to CO2 photoreduction because CO2 mass transfer has to take place in two difficult stages: from gas to liquid, where CO2 gas needs to be dissolved into the liquid and from liquid to the solid, where the soluble CO2 needs to be adsorbed on the surface of the photocatalyst. Disadvantages of the fixed bed photoreactors are the low mass-transfer rate, temperature gradient, the catalyst is difficult to replace, and the active life of the catalyst must be considered. Mozia[9] presented that when the suspended photocatalysts are used, the efficiency of photocatalysis is usually much higher than that in the fixed-bed photoreactors, and the catalysts fixed on a carrier material, such as glass, quartz, wire mesh, or membrane, cause some shortcomings, such as low surface-to-volume ratios, less photocatalytically active surface, and inefficiencies introduced by the absorption and scattering of light through the reaction medium. The drawback of membrane photoreactors is a relatively low diffusion process, which is rate-determining, with membrane fouling, and a limited lifetime of the membrane.[9,10]

It is suggested that the annular photocatalytic gas–solid fluidized bed, which has been used for other photocatalytic processes,[11,12] would be suitable for CO2 photoreduction because it presents the following intrinsic advantages: (i) CO2 is directly adsorbed onto the surface of the photocatalyst with a one-step mass-transfer process; (ii) highly efficient gas–solid contact; (iii) high mass- and heat-transfer rates; (iv) flexible operation, such as CO2 or photocatalyst recycling into the reactor column; and (v) convenient UV lamp arrangement. However, the application of the annular fluidized bed photoreactor (FBPR) for CO2 photoreduction is rarely reported.[13,14]

Computational fluid dynamic (CFD) modeling is a robust tool to numerically investigate multiphase flow and reactions. Specifically, CFD simulation relies on the experimental validation. Furthermore, it can also effectively analyze the physical mechanism and reaction behaviors to provide theoretical support for understanding the flow and reactions in the reactor. The successful simulation of gas–solid flow in a fluidized bed relies on solid particle modeling. Currently, there are two CFD methods, namely, Eulerian and Lagrangian modeling for researching the photocatalytic process. For Eulerian modeling, the solid phase is simplified and treated as the continuous phase. Thus, the kinetic theory of granular flow (KTGF) is introduced as the closure model to derive the solid pressure and solid viscosity. For example, Jing et al.[15] simulated photocatalytic hydrogen production from the refinery gas of H2S in an annular fluidized bed by Eulerian gas–solid two-fluid modeling (TFM); Pareek et al.[16] investigated pilot-scale annular bubble column photocatalytic reactors by the Eulerian gas–liquid–solid three-phase models. Generally, the Eulerian method is suitable for the simulation of the flow in large pilot-scale reactors with the number of particles being more than 109.[17] This method can save computing time and has low requirements for computer resources. The disadvantages of this method include the following: (i) solid particles are simplified as a continuous interpenetrating phase and averaged flow properties in each cell; (ii) the assumption of KTGF is only for perfectly smooth spheres, and only the normal impact and rebound velocities of colliding particles are taken into account neglecting the tangential impact; (iii) the particle velocity distribution is assumed only to match the Maxwell function equation; and (iv) it is not possible to trace individual particles or small-particle clusters in the computational domain.[18,19] For Lagrangian modeling, Braham and Harris[20] simulated the gas–solid flow in an annular fluidized bed photoreactor (FBPR) by discrete particle modeling (DPM). For DPM, the force between two particles is from elastic soft sphere collision, just like a spring. The dominated contact forces comprise the normal force and the tangential force. DPM is very close to the real physical particle collision process, but it is only suitable for simulation of the flow in small geometries because it is limited by the number of particles in the system, generally less than 106 due to the computer memory size.[20] In addition, DPM generally uses the explicit central difference time integration scheme to solve the particle motion, which is only conditionally stable with the stability determined by the size of the time-step.[21,22] The required very small time-step, for example, 10–6 s, which is dependent on the particle stiffness and particle collision contact state, increases the computational time cost and makes the simulation expensive.

CFD Eulerian–Lagrangian multiphase particle-in-cell (MP-PIC) is a newly developed robust modeling method for solid flow, which adopts the concept of discrete particles or parcels in DPM. For this method, the motion of each particle or parcel follows the equation of Newton’s second law. The collision between particles is derived from the particle stress model. The MP-PIC method can achieve greater stability than other simulation methods by treating the solid particles as computational discrete particles or as a continuum phase. In this method, the particle properties are mapped from the Lagrangian coordinates to a Eulerian grid using interpolation functions, and after the evaluation of the continuum derivative terms, the particle properties are mapped back to the individual particles.[23] The advantages of this method are as follows: (i) the MP-PIC for modeling particle collision is suitable for a dense solid flow without significant loss of accuracy;[19] (ii) the MP-PIC method can be applied for the solid flow from small reactors including several hundreds of clouds with each cloud only consisting of one particle to industrial reactors using 106 clouds/parcels to represent about 1013 particles;[24,25] and (iii) the time step for the particle motion does not rely on the particle size because the particle interforce is modeled by the particle stress model without resolving individual particle collision using the averaged Eulerian field, so the MP-PIC modeling is significantly faster than DPM. Nowadays, this method has successfully simulated the heterogeneous gas chemistry of ozone decomposition and CO2 capture by dry amine-grafted sorbents in the fluidized bed with the particle motion predicted accurately.[26,27] However, MP-PIC modeling has not been employed to numerically investigate the photocatalytic reactors. Therefore, in this work, we explore the use of the FBPR for CO2 photoreduction by investigating the flow behavior and reaction process with MP-PIC modeling. Accordingly, this paper focuses on this parcel-in-cell method to understand the physical behavior and reaction mechanism of CO2 photoreduction in the FBPR.

FBPR for the CO2 Photoreduction Process

A typical photocatalytic gas–solid fluidized bed is shown in Figure . For the annular column, the inner and outer diameters are 26 and 40 mm, respectively, as reported in ref (28). The gas distributor is made from a sintered glass plate. The photocatalyst is TiO2 coated on the particle of γ-Al2O3. The particle size is 107 μm and the particle density is 3246 kg/m3. The minimum fluidization velocity is 1.16 cm/s. An UV lamp is in the bed center. The details of the UV lamp are described in Section . The photocatalysts are packed in the annular column at the initial state before fluidization commences. The gas is composed of an inert carrier gas such as argon, and the reactant gases such as CO2 and water (H2O) moisture. The molar percentages of Ar, CO2, and H2O for the gas inlet are 72.36, 25.05, and 2.59%, respectively, which are calculated from the Ar, CO2, and H2O partial pressures.[29] The H2O moisture was generated from a temperature-controlled saturator.[29] This setting is following Thompson et al.’s work.[29] The gas mixture is introduced into the column from the bottom at a constant velocity. Then, the photocatalyst particles are suspended fully if the gas velocity exceeds the minimum fluidization velocity and the bubbles flow from the bottom, rising through the bed and causing particles to be highly mixed and agitated. CO2 and H2O are adsorbed on the surface of TiO2 and activated to be reduced by the irradiation of UV light. Finally, the product gases, such as CH4, H2, and CO, and unreacted gases flow out of the column from the top.

Figure 1

Schematic diagram for the FBPR.

Simulation Theory

In this work, the Lagrangian MP-PIC modeling was compared with Eulerian TFM. For MP-PIC modeling, the closure model is the Lun equation as the particle stress model with e = 0.95. The drag model is the Ergun-WenYu (Gidaspow) equation, and the isotropy model is applied if there is no special mention or explanation in the following sections. Here, the isotropy model is employed because we consider that the additional effect of particle collision drives the particle velocity distributions toward isotropy in the solid dense phase.[30] For TFM, the closure model is the solid pressure using the equation from the Lun model and the solid viscosity using the equation from the Gidaspow model, and the drag model using the Ergun-WenYu (Gidaspow) equation. The detailed equations are presented below.

CFD Model

For Eulerian TFM, the Navier–Stokes equations were used, where the continuity equation iswhere ε is the k phase volume fraction, dimensionless; ρ is the phase density, kg/m3; u is the velocity, m/s; R is the reaction mass source, kg/(m3·s); t is the time, s; and x is the ith coordinate direction, m.

The momentum equation is given below.where P is the pressure, Pa; Fdrag,gs is the drag force between gas phase and particles, N; S is the momentum source, N/m3; g is the acceleration due to gravity, 9.81 m/s2; and τ is the stress tensor, Pa.

The species equation iswhere X is the mass fraction of n species, dimensionless and D is the diffusion coefficient, kg/(m·s).

For MP-PIC modeling, the equations for the gas phase are the same as the Eulerian models, but for the solid particles, the equations are[31,32]where a is the acceleration, m/s2; Fdrag,gs is the drag force between the gas phase and particles, N; Fdrag,gs = Dgs·(ug – us); Dgs is the drag coefficient, kg/(m3·s); and τs is the particle stress, Pa. The drag forces in eqs and 6 are equal due to the momentum exchange between the gas phase and particles.

The particles in a cell are averaged by the dual-mesh Lagrangian procedure and a barycentric interpolation scheme is used.

Radiative Transport Equation

For UV light radiation, the radiative transport equation (RTE) was used to describe the light radiation intensity distribution in the bed. The RTE iswhere I is the radiative intensity at a given position following the Ω direction, W/m2; κs, absorption coefficient, 1/m; βs, extinction coefficient, 1/m; and σs, scattering coefficient, 1/m. Here, for TiO2 particles, the absorption coefficient is κs (1/m) = 3.598 × Wcat (g/m3) and the scattering coefficient is σs (1/m) = 0.2758 × Wcat (g/m3) and Wcat is the photocatalyst loading.[33] The physical meaning of RTE (eq ) is that the change of intensity is equal to the intensity enhancement due to emission, intensity reduction due to absorption and intensity enhancement due to scattering. Ib is the intensity from emission in the temperature field. The particle concentration influences the UV radiation through the absorption coefficient and scattering coefficient. During simulation, Wcat is equal to the solid density of the photocatalyst times solid fraction (ρs × εs). Alvarado-Rolon et al.’s work on kinetic modeling of paracetamol degradation by photocatalysis considered the UV light absorptions at 254 nm by both the photocatalyst of TiO2 and the reagent of paracetamol.[34] In this work, the UV absorption by CO2 is not considered because He et al.[35] demonstrated that no absorption by CO2 is observed in the wavelength range of 307–725 nm. Moreover, the UV absorption by H2 is omitted in this work because the concentration of H2 in the gas is very low with the maximum mass fraction of about 10–10.

Generally, there are two methods for solving this equation. One is the discrete ordinate method (DOM), also known as the P–N method.[36−38] The full angle can be divided into a series of number of discrete angular intervals and then resolve RTE along different angular directions. Another method is the P1, which is the simplest approximation of P–N models. It expands the RTE into a spherical isotropic model. Because the DOM method takes a long time and high computer resources, the P1 method is used in this work, which has been employed by other simulations.[39,40]

The intensity of radiation on the inner wall as the boundary condition can be described as the following equation.[41]where Rin is the radius of inner wall, m; Hf is the height of fluidized bed column, m; L is the height of the lamp, m; z is the size along the vertical direction, m; S1 is the UV lamp radiation function, S1 = 2π × Rlamp × Ilamp; Rlamp is the radius of the lamp; and Ilamp, radiation intensity from the lamp. In this work, Rin = 0.013 m; Hf = 0.1 m; Rlamp = 0.008 m; and L = 0.08 m.

The light intensity from the UV lamp is set as 4000 W/m2 (400 mW/cm2) at a wavelength of 365 nm, which is the same as that in Thompson et al.’s work.[29] The typical radiation intensity distributions with Ilamp = 4000 W/m2 are shown in Figure . The purpose of Figure is to compare the radiation intensity distribution between a nearly empty column and a specific fluidized bed as an example. Figure a shows the distribution in the nearly empty bed with a very low solid fraction of 10–10 to imitate the rare tiny photocatalyst particles may float in the bed and Figure b shows the distribution in the dense fluidized bed with solid particles under the condition of Ug = 0.08 m/s and initial packed bed height, Hpack = 26 mm. In the nearly empty bed, the highest radiation intensity is in the bed center area. Along the vertical direction, the radiation intensity decreases along the direction from the bed center to the bed top and bed bottom. Along the horizontal direction, the radiation intensity decreases from the inner wall to the outer wall. These observations agree with the prior simulation and experimental work.[42] In the fluidized bed with a bed expansion height of Hb = 40.2 mm and average solid fraction, εs = 0.26, the radiation intensity attenuates highly in the solid photocatalyst particles, where the attenuation degree is a function of the solid concentration. Therefore, it can be concluded that the CO2 photoreduction reaction mainly occurs in the thin area next to the inner wall.

Figure 2

Radiation intensity distribution in the FBPR: (a) bed with rare tiny photocatalyst particles; (b) dense fluidized bed.

Reaction Model

CO2 photoreduction is a complex reaction system, and the reaction mechanism and kinetics have been investigated.[43−46] Generally, CO2 photoreduction is a mass-transfer-controlled process and thus the kinetic equation is widely accepted to be based on the Langmuir adsorption equation, rather than the Arrhenius rate equation.[44] For the reaction model, the Thompson model based on the Langmuir equation is applied.[29] Thompson et al.’s publication[29] provided many informative data and parameters, which best fit to the simulation in this work and therefore, we employed parameters from Thompson et al.’s work. Furthermore, we are planning to carry further validation using cases and parameters from other publications and experiments. The Thompson model equation is the Weibull probability density function (PDF) multiplied with the Langmuir adsorption equation form as the reaction rate. In the Weibull PDF, the time and radiation intensity are introduced. Therefore, the reaction rate is the function of time, radiation intensity, and CO2 partial pressure and H2O partial pressure.

The reaction rates for generating CH4, CO and H2 arewhere r is the reaction rate on the surface, μmol/(gcat·h); k is the rate constant, μmol/(gcat·h); PDF(t) is the Weibull PDF; K is the equilibrium adsorption constant, 1/bar; and P is the partial pressure, bar.

The Weibull PDF equation iswhere ηd is the deactivation scale parameter, dimensionless and α is the reaction order of light intensity, dimensionless. The Weibull PDF equation introduced by Thompson et al.[29] considers that the active sites on the surface of the photocatalyst deactivate over time.

In this work, all the rection parameters for eqs –12 are picked up from ref (29). Generally, the CO2 photoreduction route is from CO2 to CO and finally to CH4. CO is a very important intermediate product and CH4 is a target product as the solar fuel. The mechanism of CO2 photoreduction to CH4 has been theoretically studied and demonstrated by Ji and Luo.[43] H2 is produced from water splitting. Other products, such as formaldehyde and methanol, are the intermediate products and possible final products.[47,48] It has been analyzed in the publication.[45] From Thompson et al.’s experiments,[29] CO, H2, and CH4 are the main final products and other publications have demonstrated it as well.[46] Therefore, in this work, we consider CO, H2, and CH4 as the main final products by CO2 photoreduction.

The volumetric reaction rate, rv, in the modeling is

Here, we consider only the half particle area with a glazed side absorbs the UV light. As = Aps/Vp. As is the solid surface area per volume, m2/m3; Aps is the particle surface area, m2; and Vp is the particle volume, m3.

Simulation Method

The simulation studies were performed using OpenFOAM, which is a free and open-source toolbox providing different CFD robust solvers, such as multiphase flow and reaction solver, particle tracking solvers, and so forth, and providing the accessible and modifiable code interface. The assumptions used in this work are as follows: (i) the gas–solid flow in the annular column is isotropic and therefore, the 2D axisymmetric model and mesh are used in this work; (ii) the reactor is isothermal and therefore, the energy equation is not considered and applied in the simulation; (iii) the gas is incompressible with constant density and constant viscosity (gas density: ρg = 1.54 kg/m3; gas viscosity: μg = 2.64 × 10–5 Pa·s); (iv) the gas flow is in the laminar regime (Reg = 30–100); (v) the gas velocity is uniform for the gas inlet area; (vi) the wall for gas flow is no slip wall; (vii) one parcel has only one particle and the particle size is a unique fixed value; and (viii) the CO2 photoreduction occurs after t = 0.5 s so that the particles are fully fluidized to prevent the impact of the flow transition from the fixed state to the fluidized state on the CO2 photoreduction reaction. Heat can be generated by the UV radiation, but one of the advantages of the fluidized bed is that very fast heat transfer occurs and thus, it can achieve a nearly uniform temperature distribution. Moreover, Thompson et al.[29] reported no significant temperature rise caused by UV lighting. Therefore, in this work, we assumed that the temperature (314 K) is uniformly distributed in the bed. During simulation, we assume that all the boundary walls are transparent quartz walls and their influence on the UV irradiation is negligible. Thus, the boundary walls are set as the same as the internal field.

In this work, the 2D axisymmetric geometrical model used in the simulation is shown in Figure . The height is 100 mm, while the annular length is 7 mm, as shown in Figure a. The half single-layer angular angle is 0.8°. The mesh size should be more than one particle size, otherwise it may lose the physical meaning. The mesh for Lagrangian MP-PIC modeling is 14 × 200 and the size for each cell is 5 times of the particle diameter, dp, as shown in Figure b. The mesh for Eulerian TFM is 28 × 400, as shown in Figure c. The ratio of the mesh size to particle size is suggested from a prior published work.[49,50]

Figure 3

2D axisymmetric model and meshes: (a) 2D geometrical model; (b) mesh for Lagrangian MP-PIC modeling; (c) mesh for Eulerian TFM.

The discretization scheme is the Gauss limited linear. The time step is 10–4 s. The courant number, C, representing a particle stay in one cell of the mesh, C = uΔt/Δh, is 0.016 < 0.7 as the gas velocity is 0.08 m/s, which satisfies the convergence condition. Δt is the simulation step time and Δh is the mesh height. Each step simulation tolerance is less than 10–5. The numbers of particles are 78,390, 114,660, and 150,930 for the initial packed bed heights of 0.026, 0.038, and 0.050 m with a solid fraction of 0.6, respectively, which follows the rule of (nparticle × Vp/Vpacked bed = 0.6). nparticle is the number of particles; Vp is the particle volume; and Vpacked bed is the volume of the initial packed bed. The maximum solid fraction for the packed bed is 0.63. It means that the local solid fraction in the dense flow would not exceed the maximum solid fraction and at this limitation, no particle collision occurs. The simulation time duration for evaluating the simulation parameters is 10 s. The simulation time duration for investigating operating parameters is 20 s.

Validation of Simulation

Validation Data

We validated the simulation using a bed expansion height, H, bubble size, Db, and mass fraction of CH4, YCH4, based on a gas superficial velocity of 0.08 m/s. The experimental values of the bed expansion height and bubble size are derived from the empirical correlations from ref (28).

The average bed voidage is presented belowwhere ε is the bed voidage, dimensionless; εe is the emulsion phase voidage, dimensionless; εb is the bubble phase voidage, dimensionless; and δ is the bubble phase fraction, dimensionless.

The bubble phase fraction iswhere U is the gas superficial velocity, m/s and Umf is the minimum fluidization velocity, m/s.

The emulsion phase voidage is given here

The bubble phase voidage is

The bed expansion height is presented herewhere Hb is the bed expansion height, m; Hpack is the packed bed height, m; εpack = 0.6; and Hpack = 26 mm.

Due to the above correlations, it is derived that when the gas superficial velocity is 0.08 m/s and the bed expansion height is 41.1 mm.

The bubble generation and bubble rising in the bed are typical phenomena of gas–solid bubbling-fluidized beds. The bubble size, Db, can be estimated by the following correlation

The bubble size, Db, is 4.1 mm when the gas superficial velocity is 0.08 m/s and the bed height, Hr, the height between a point on the reactor bed column and the bed bottom, is 30 mm.

The mass fraction of CH4, YCH4, is 10–7 to 10–4 through calculation from ref (29).

Effect of Particle Stress Models for MP-PIC Modeling

For MP-PIC modeling, the interaction force between particles is described by the particle stress. Therefore, the particle stress equation is crucial in accurately predicting the particle motion in the fluidized bed. In OpenFOAM, there are two models for describing the particle stress in the bed, which are the Lun model and Harris–Crighton model listed in the followings.

Lun model equationwhere τs is the interparticle stress, Pa; e is the coefficient of restitution, dimensionless; εs,pack is the solid volume fraction of the packed bed, dimensionless; and σ is the root mean square (RMS) of velocity fluctuation.

Harris–Crighton model equationwhere ps and β are the coefficients for this model. Generally, ps = 10 Pa and β = 2.

Figure shows the particle distribution and gas fractions predicted by Lagrangian MP-PIC modeling with the Lun model and the Harris–Crighton model. In the simulation, e = 0.95 which is referred to refs (51–53), and the drag model is the Ergun-WenYu (Gidaspow) equation. The isotropy model is applied. Figure a shows the particle motion behaviors in the reactor predicted by MP-PIC modeling and where each particle motion can be traced. The purpose of Figure a is to show the prediction capability of MP-PIC modeling for particle cloud motion. Because of massive particle gathering, the solid fraction distribution cannot be identified clearly from this figure. However, the figures of the gas fraction distribution can indirectly display the solid fraction distribution. Therefore, Figure b,c for the gas fraction distribution were drawn and indirectly show the solid fraction distribution. The predicted bed height value is achieved from the averaged gas fraction distribution, as shown in the left picture of Figure b,c; thus in case to prevent the measurement error caused by the fluctuation of the bed surface. The bubble motion in the bed can be observed in the right picture of Figure b,c at t = 10 s. Finally, the corresponding bed expansion heights are 40.2 and 44.1 mm for the Lun model and Harris–Crighton model, respectively. It is found that MP-PIC modeling with the Lun model estimates the bed height very well with an error of 2.2%. The error of bed height predicted by MP-PIC modeling with the Harris–Crighton model is 7.3%, much higher than the real value. It can be concluded that the MP-PIC modeling with the Lun model is a suitable model for describing gas and solid flow in the FBPR through comparing the simulation results with experimental values of the bed expansion height from the correlation equation, eq .

Figure 4

Comparison of particle stress models by Lagrangian MP-PIC modeling: (a) particle distribution with the Lun equation as the particle stress model; (b) gas fraction distribution with the Lun equation as the particle stress model; (c) gas fraction distribution with the Harris–Crighton equation as the particle stress model.

Effect of Drag Models for MP-PIC Modeling

In order to investigate the drag models acting on the simulation of the FBPR, three drag models, namely, Ergun-WenYu (Gidaspow) model, Plessis–Masliyah model, and WenYu model were employed. The Ergun-WenYu drag model combines the Ergun equation when εg ≤ 0.8 and WenYu equation when εg > 0.8, as shown in eqs and 23. The Plessis–Masliyah model is based on the Ergun equation of eq , but the coefficients A and B are estimated by eqs and 25, respectively.

Ergun equationwhere Dgs is the interphase momentum exchange coefficient, kg/(m3·s); dp is the particle size, m; and A and B are coefficients, A = 150 and B = 1.75.

WenYu equationwhere CD is the drag coefficient.

Plessis–Masliyah equation

During simulation, the Lun model with e = 0.95 as the particle stress model and the isotropy model were applied. The simulation results of gas fraction distributions are shown in Figure , in which Figure a for the Ergun-WenYu drag, Figure b for the Plessis–Masliyah drag and Figure c for the WenYu drag. The corresponding estimated bed expansion heights are 40.2, 42.5, and 40.2 mm, respectively. In contrast to the experimental correlation value of bed height, 41.1 mm, the prediction error by the Ergun-WenYu drag and WenYu drag is only 2.2% but that by the Plessis–Masliyah drag is 3.4%. Additionally, the WenYu drag is generally suitable for the dilute phase with solid volume fractions below 0.2 and the Ergun-WenYu drag is suitable for solid volume fractions up to the packed state. Thus, it can be concluded that the Ergun-WenYu drag model is the best for predicting the particle flow in the annular FBPR because the predicted bed expansion height in this case is closest to the experimental value in the dense phase.

Figure 5

Effect of different drag models on the gas fraction distribution: (a) Ergun-WenYu (Gidaspow) drag model; (b) Plessis–Masliyah drag model; (c) WenYu drag model.

In a summary, the bed expansion height predicted by MP-PIC modeling is highly influenced by the particle stress models and drag models. It is demonstrated that the Lun particle stress model and the Ergun-WenYu drag model are suitable for describing the particle motion in the FBPR of this work.

Effect of Wall Boundary for MP-PIC Modeling

The interaction between the particles and the wall can be described by hard sphere collisions with the wall.[18] The normal collision velocity and tangential velocity are presented as follows

For the normal directionwhere ew is the coefficient of restitution for the particle and wall, dimensionless and u is the velocity with the normal direction of the wall, m/s.

For the tangential directionwhere ut is the velocity with the tangential direction to the wall, m/s and μw is the coefficient of friction for the wall, dimensionless.

The coefficient of restitution for collision between the particles and wall, ew, and the coefficient of friction between particles and wall, μw, are key parameters, thus the effects of these parameters on the flow behaviors in the bed have been studied. The results of gas fractions with different ew and μw are shown in Figure . Figure a–c are shown under different ew with the same μw of 0.6. Figure d–f are shown under different μw with the same ew of 0.9. It indicates that all the bed expansion heights are the same, as 40.2 mm, but the bubble sizes are different. It is discovered that the wall boundary does not affect the bed expansion height but influences the bubble size. In this work, the bubble size is determined by averaging the bubble vertical size and horizontal size.

Figure 6

Effect of wall function parameters on the gas fraction distribution. Coefficient of restitution (μw = 0.6): (a) ew = 0.8, (b) ew = 0.9, and (c) ew = 0.95; coefficient of friction (ew = 0.9): (d) μw = 0.09, (e) μw = 0.35, and (f) μw = 0.6.

It can be seen from Figure a–c that when ew increases from 0.8 to 0.95, the bubble size decreases from 4.6 to 4.3 and 3.5 mm, respectively. It demonstrates that ew = 0.9 is the best for predicting the particle flow dynamics in the FBPR with the lowest prediction error. Figure d–f shows the particle flow under different μw = 0.6, 0.35, and 0.09. With the wall coefficient of friction decreasing, the bubble size increasing from 4.3 to 5.2 mm. As μw = 0.6, the predicted bed expansion height matches the experimental value with the lowest error of 2.2%.

It concluded that the wall function parameters affect the bubble size heavily. This is because the ratio of the bubble size to the narrow annular size is 4/7. Therefore, the wall zone influences the bubble formation and bubble coalescence significantly. Due to the bubble motion and size, it is concluded that ew = 0.9 and μw = 0.6 are the optimum values for describing the particle flow in the studied FBPR because the predicted bubble size in this case is closest to the experimental value.

Comparison between TFM and MP-PIC

It is well known that the particle flow dynamics in the bubbling fluidized bed can be efficiently simulated by TFM as well. For comparing the effectiveness of TFM and MP-PIC modeling on solid flow in the FBPR, the simulation by TFM was performed with the Ergun-WenYu drag model, e = 0.95 and the particle pressure using the Lun equation, and these conditions and models are the same with those by MP-PIC modeling. The comparison between TFM and MP-PIC modeling is shown in Figure . Figure a shows the solid distribution predicted by TFM. Figure b shows the particle distribution predicted by MP-PIC modeling and for observing the bubble clearly, the gas fraction estimated by MP-PIC modeling is shown in Figure c. The bed expansion heights predicted by TFM and MP-PIC modeling are 37.0 and 40.2 mm, respectively, with prediction errors of 10.0% and 2.2%, respectively. The bubble sizes predicted by TFM and MP-PIC modeling are 3.3 and 4.3 mm, respectively, with prediction errors of 19.5 and 4.9%, respectively. The bed expansion height and the bubble size estimated by TFM are much smaller than the real values. It expresses that MP-PIC modeling on the FBPR is better than TFM.

Figure 7

Comparison of simulations with Eulerian TFM and Lagrangian MP-PIC modeling: (a) solid fraction distribution predicted by TFM; (b) particle distribution predicted by MP-PIC modeling; (c) gas fraction distribution predicted by MP-PIC modeling.

In addition, the differences of particle distributions simulated by TFM and MP-PIC modeling are on the locations of the bed bottom and bed top, as shown in the average solid fraction for TFM and average gas fraction for MP-PIC modeling. The MP-PIC model predicts that the bottom of the bed is a trapezoidal air cavity, while TFM predicts that the trapezoidal air cavity at the bottom of the bed is divided into two parts by a particle strip. The MP-PIC model predicts that the top of the bed is in the shape of a shallow plate and the side wall is a high-concentration area of particles, while TFM predicts that the top of the bed is in the shape of a goblet. The specific particle distribution at the top and bottom of the bed needs to be further verified by experimental imaging technology.

In order to further inspect the parameter settings on MP-PIC modeling, the bed expansion heights at the gas superficial velocities of 0.08, 0.16, and 0.24 m/s were used to validate the MP-PIC simulation by the empirical correlation equations. From the correlation equations, the bed expansion heights are 41.1, 44.4, and 47.4 mm, respectively, and the corresponding expansion height predicted by MP-PIC modeling are 40.2, 43.3, and 48.8 mm at the velocities of 0.08, 0.16, and 0.24 m/s, respectively. The prediction errors are 2.2, 2.5, and 2.9%, respectively. These results further prove that the current settings for MP-PIC models are the best fit to the experiment.

CO2 Photoreduction Process

The produced CH4, CO, and H2 mass fractions from CO2 photoreduction predicted by MP-PIC modeling using the Ergun-WenYu drag model under ew = 0.9 and μw = 0.6 are shown in Figure . The gas velocity is 0.08 m/s. It can be seen that CO2 photoreduction occurs mainly in the area near the inner wall because it receives UV light. Figure also shows that the local reaction intensity and product concentration have a close relationship with the local UV light intensity and solid photocatalyst concentration. The maximum CH4, CO, and H2 mass fractions are 3.1 × 10–7, 8.5 × 10–7, and 7.6 × 10–10, respectively. The maximum CH4 mass fraction is in the range of the experimental values reported[29] and therefore, it validates the reaction model and parameter settings used in this work.

Figure 8

Production of CH4, CO, and H2 through CO2 photoreduction predicted by MP-PIC modeling: (a) CH4, (b) CO, and (c) H2.

Simulation Results and Discussion for Operation Parameters

In this section, the simulations were performed under the following conditions: Ergun-WenYu (Gidaspow) drag model, Lun particle stress model with e = 0.95, wall coefficient of restitution with ew = 0.9, and wall coefficient of friction with μw = 0.6.

Effect of Gas Velocity on Solid Flow and CO2 Photoreduction

Regulating gas velocity is commonly used for controlling solid flow and CO2 photoreduction. The simulations here were carried out with gas velocities of 0.08, 0.16, and 0.24 m/s. As the gas velocity rises, the bed expansion height increases significantly from 40.2 to 43.3 and 48.0 mm with the bubbles going through the bed center. With the gas velocity increasing, the bubbles coalesce forming a long-crooked strip in the bed center area. In order to investigate the solid fraction distribution along the radial position, the simulation data were averaged from t = 0.5 s to t = 20 s. The averaged solid fraction distributions at bed height Hr = 30 mm are shown in Figure a. It shows that the solid concentration near the wall around εs = 0.54 is higher than that in the bed center (εs = 0.1–0.3). Weber and Mei[54] demonstrated that the solid fraction near the wall is around 0.5–0.6 and the solid fraction in the bed center is about 0.3 when Ug/Umin = 6. This demonstrated that our current simulation with Ug/Umin = 6.9 when Ug = 0.08 m/s is reasonable, as the solid concentration in the bed center (εs = 0.1–0.3) and near the wall (εs = 0.54), which is in coincidence with the experimental results in ref (54) with similar conditions. It further indicates that the settings of simulation parameters in the MP-PIC model of this work are fitted with predicting the particle flow in the reactor. The solid concentration near the inner wall is slightly higher than that near the outer wall. For example, when the gas velocity is 0.08 m/s and εs are 0.55 and 0.51 for the inner wall and outer wall, respectively. This may be caused by the axisymmetric annular structure with a small inner wall area and large outer wall area. This phenomenon was also demonstrated by Khan and Shamim.[55]

Figure 9

Effect of gas velocity on CO2 photoreduction: (a) solid fraction distribution along the radial direction at Hr = 30 mm under different gas velocities; (b–d) generated CH4, CO, and H2 mass fractions with time at the outlet for different gas velocities of (b) 0.08, (c) 0.16, and (d) 0.24 m/s.

The influence of gas velocity on CO2 photoreduction relies on the solid photocatalyst distribution, especially near the inner wall and therefore, high gas velocity may increase CO2 photoreduction. With gas velocity increasing, the enhanced bed expansion enlarges the solid phase area receiving UV light. However, this also decreases the solid concentration near the wall and may reduce CO2 photoreduction. Therefore, the CH4 yield from CO2 photoreduction is a function of local radiation intensity and local photocatalyst concentration. This is due to the radiation intensity not being distributed uniformly, with the highest value in the middle of the reactor height and lowest in the top side and bottom side of the reactor. The solid photocatalyst concentration distribution is strongly determined and can be controlled by the gas velocity. Therefore, it seems that CO2 photocatalysis can be optimized by the gas velocity to change the photocatalyst concentration distribution along the bed height.

The variations of CH4, CO, and H2 mass fractions at the gas outlet with time are shown in Figure b–d. When the time started from t = 0.5 s, the CH4, CO, and H2 mass fractions increase until t = 20 s. However, the simulation results indicate that the mass fractions of CH4, CO, and H2 at the gas outlet are affected by the gas velocity very slightly. For example, at t = 20 s, CH4 mass fractions at the gas outlet for Ug = 0.08, 0.16, and 0.24 m/s are 0.49 × 10–7, 0.44 × 10–7, and 0.48 × 10–7, respectively; CO mass fractions at the gas outlet are 1.35 × 10–7, 1.19 × 10–7, and 1.32 × 10–7, respectively; H2 mass fractions at the gas outlet are 1.16 × 10–10, 1.06 × 10–10, and 1.12 × 10–10, respectively. The possible reason for this situation is that the initial packed bed height remains unchanged. As the gas velocity increases, the bed expansion height increases with enhancing the possibility of receiving UV light, but the corresponding particle concentration decreases locally and the total photocatalyst surface area correspondingly decreases for receiving light, resulting in the concentration of CH4 generated by CO2 photoreduction at the outlet basically unchanged with above different gas velocities. Furthermore, the oscillation of the mass fraction is intensified by higher gas velocity. The oscillation of the mass fraction, as shown in Figure b–d, is caused by the variation of the solid concentration and pressure under the gas flow, which resulted from the nonlinear dynamic motion of heterogeneous flow structures in gas–solid fluidization. Cui et al.[56] demonstrated that with increasing gas velocity, the oscillation frequency raises, and the amplitude enlarges.

For further demonstrating the effect of gas velocity on CO2 photoreduction, the simulations were carried out based on the initial packed bed height of 38 mm with Ug = 0.08, 0.16 and 0.24 m/s. For different Ug, the corresponding CH4 mass fraction is 0.72 × 10–7, 0.66 × 10–7 and 0.66 × 10–7, respectively. Similarly, the corresponding CO mass fraction is 2.09 × 10–7, 1.84 × 10–7 and 1.85 × 10–7, respectively. The H2 mass fraction is 1.59 × 10–10, 1.50 × 10–10, and 1.52 × 10–10, respectively. This further indicates that the gas velocity influences the product yield very slightly. Only at a low velocity of 0.08 m/s, the CH4 yield is slightly higher (8%) than that at higher gas velocity.

In brief, CH4 generated by CO2 photoreduction in the FBPR indicates that CO2 photoreduction process mainly occurs near the inner wall. The gas velocity influences the generation of CH4, CO, and H2 very slightly.

Effect of Initial Packed Bed Height on Flow and CO2 Photoreduction

The packed bed height is the height between the top and the bottom of the filled photocatalysts in the reactor column at Ug = 0 m/s, which physically means the photocatalyst loading in the bed. Catalyst loading can affect CO2 photoreduction by: (1) increasing photocatalyst concentration in the bed and (2) increasing the area of the solid phase exposing to the UV light. The simulations were performed with the initial packed bed heights of 26, 38, and 50 mm.

The simulation results display that the bed expansion height increases sharply from 40.2 to 64.6 and 82.7 mm with the initial bed height increasing from 26 to 38 and 50 mm under the same gas velocity, Ug = 0.08 m/s. The bed expansion height is over the half bed height which is the position with the maximum radiation intensity. Figure shows the simulation results for the solid flow and photoreaction. The solid fraction distributions at the bed height Hr = 30 mm are shown in Figure a. Although the initial packed bed height rises from 26 to 50 mm, the solid fraction distribution at Hr = 30 mm is nearly the same with εs ≈ 0.54 near the wall and εs ≈ 0.24 in the bed center.

Figure 10

Effect of the initial packed bed height on CO2 photoreduction: (a) solid fraction distribution along the radial direction at Hr = 30 mm with different initial packed heights; (b–d) generated CH4, CO, and H2 mass fractions with time at the outlet for different initial packed bed heights of (b) 26, (c) 38, and (d) 50 mm.

Figure b–d shows the variation with the time of mass fractions of CH4, CO, and H2 at the gas outlet. It can be seen clearly that the CH4 mass fraction obviously increases with the initial packed bed heightened. For example, under Ug = 0.08 m/s and different initial packed heights of 26, 38, and 50 mm, the CH4 mass fractions at the gas outlet are 0.49 × 10–7, 0.72 × 10–7, and 1.02 × 10–7, respectively; the CO mass fractions at the gas outlet are 1.35 × 10–7, 2.09 × 10–7, and 2.96 × 10–7, respectively. The H2 mass fractions at the gas outlet are 1.16 × 10–10, 1.59 × 10–10, and 2.21 × 10–10, respectively. This indicates that the higher initial packed bed height, the higher product yield.

It can be concluded that increasing the initial packed bed height or photocatalyst loading enhances CO2 photoreduction to generate more CH4 in the FBPR. This is due to (i) the increase of the photocatalyst loading extending the surface area for CO2 adsorption and photoreaction and (ii) the increase of the bed height leading to more particles receiving the highest intensity of UV light because the radiation intensity is not distributed uniformly with the highest value in the middle of the reactor height.

Conclusions

This work investigated solid photocatalyst flow and CO2 photoreduction in a FBPR by CFD Lagrangian MP-PIC modeling. The simulation on the hydrodynamics of particle flow in the fluidized bed has been validated by the empirical correlations. These empirical correlations for predicting the particle flow and the reaction kinetic equations have been demonstrated to be effective, as shown in refs (28) and (29), respectively. In this work, in addition to validating the simulation by empirical correlations or experiments, we also demonstrated that CFD can be used to analyze the flow and reaction behaviors. The results reported here can provide guidance for the reactor design and operation of CO2 photoreduction using the FBPR.

First, different particle stress models, drag models, and simulation parameters, such as coefficient of restitution, ew, and coefficient of friction, μw, have been studied and evaluated for the simulation of particle flow in the FBPR. The results show that the Lun particle stress model with e = 0.95 and the Ergun-WenYu (Gidaspow) drag model with ew = 0.9 and μw = 0.6 are the best selection to describe the solid flow, particle interforce, and particle intercollision, which is validated by the experimental values of the bed expansion height and bubble size from the empirical correlation equations. It also indicates the model selections discovered and matched the intrinsic mechanism and rules of photocatalyst flow in the annular fluidized bed. Moreover, the Lagrangian MP-PIC modeling is compared with the Eulerian TFM in this work. It is concluded that the MP-PIC method is better than the TFM method in the field of the bed expansion height and bubble size. TMF treats the solid flow as a continuous phase using an uniform average flow characteristic under each cell and ignores the heterogeneity of particle clustering in a cell. This may cause the distortion of the simulation results. This has been discussed in detail in ref (57). The MP-PIC model treats the particles as the discrete sphere clouds and considers the difference of particle movement in a cell. Therefore, MP-PIC modeling is advantageous over TFM with considering a heterogeneous particle motion under one cell. The CH4 mass fraction by CO2 photoreduction was validated by Thompson et al.’s experiments.

Furthermore, the effects of operating parameters, such as the gas velocity and initial packed bed height on CO2 photoreduction in the FBPR, have also been investigated. The variation of the gas velocity influences the production of CH4, CO, and H2 very slightly, but the oscillation of the mass fraction increases with the gas velocity. The mass fraction of CH4, CO, and H2 increases at the gas outlet when the solid catalyst load is increased.

It can be seen that CO2 photoreduction mainly occurred in the area next to the inner wall. The maximum radiation intensity is located on the inner wall at the half of the bed height. This suggests that the solid catalyst loading should be more than half of the reactor column height. In addition, the reactor design and operating parameters for increasing the solid concentration on the inner wall can be optimized for improving CH4 productivity due to the increase of the photocatalyst particle concentration on the wall and enlargement of UV light adsorption and reaction surface area.

The CO2 photoreduction in the annular FBPR is a complex system, which is influenced by many factors and parameters from material physical properties, reactor dimensions, UV lighting, as well as operation conditions and reaction kinetics. Nevertheless, this work has successfully demonstrated that theoretical CFD models in combination with the UV radiation model and reaction kinetic model can help us to unravel the physical and reaction mechanisms of the CO2 photoreduction process.

The MP-PIC modeling developed here is an effective robust tool not only to understand the physical behavior of photocatalysts but also to explore the reaction mechanism of CO2 photoreduction process in the FBPR. The simulation method and results reported here can help optimize the design and operation of the FBPR, where the following factors should be considered: (i) the catalyst loading or initial packed bed height enhances the CO2 photoreduction to CH4; (ii) the gas velocity only increases the CH4 generation oscillation; (iii) the internal wall structures, such as the internal ring, are set up to increase the photocatalyst concentration near the wall and then enhance the CO2 photoreduction; and (iv) the optimum arrangement of the UV lamp around the inner and outer walls enhances the UV absorption by the photocatalyst.