**Simulation-Based Design and Economic Evaluation of a Novel Internally Circulating Fluidized Bed Reactor for Power Production with Integrated CO2 Capture**

#### **Jan Hendrik Cloete 1,\*, Mohammed N. Khan 2, Schalk Cloete <sup>1</sup> and Shahriar Amini 1,2,\***


Received: 10 September 2019; Accepted: 30 September 2019; Published: 11 October 2019

**Abstract:** Limiting global temperature rise to well below 2 ◦C according to the Paris climate accord will require accelerated development, scale-up, and commercialization of innovative and environmentally friendly reactor concepts. Simulation-based design can play a central role in achieving this goal by decreasing the number of costly and time-consuming experimental scale-up steps. To illustrate this approach, a multiscale computational fluid dynamics (CFD) approach was utilized in this study to simulate a novel internally circulating fluidized bed reactor (ICR) for power production with integrated CO2 capture on an industrial scale. These simulations were made computationally feasible by using closures in a filtered two-fluid model (fTFM) to model the effects of important subgrid multiphase structures. The CFD simulations provided valuable insight regarding ICR behavior, predicting that CO2 capture efficiencies and purities above 95% can be achieved, and proposing a reasonable reactor size. The results from the reactor simulations were then used as input for an economic evaluation of an ICR-based natural gas combined cycle power plant. The economic performance results showed that the ICR plant can achieve a CO2 avoidance cost as low as \$58/ton. Future work will investigate additional firing after the ICR to reach the high inlet temperatures of modern gas turbines.

**Keywords:** chemical looping combustion; power production; carbon capture; internally circulating reactor; reactor design; fluidization; techno-economics; computational fluid dynamics; filtered two-fluid model; coarse-grid simulations

#### **1. Introduction**

Several high-profile studies have shown that carbon capture and storage must play a central role in the future energy mix to reach the goal of limiting the global temperature increase to well below 2 C above preindustrial limits at a reasonable cost [1–3]. A low-cost pathway to limiting global CO2 emissions will be essential to prevent the negative consequences of climate change, while allowing for continued development in developing nations where billions of people still live in poverty.

Many different technologies have been proposed to capture CO2 from fossil-fuel power plants, after which the CO2 can either be stored or utilized in other industrial processes. However, a major challenge of such processes is the energy penalty associated with CO2 capture. An increased energy penalty requires more fuel to be used to achieve the same power output, increasing operating and capital costs, but also increasing the amount of CO2 that must be dealt with.

A promising group of technologies for capturing CO2 are those based on chemical looping combustion (CLC) [4], as they can essentially eliminate the energy penalty of CO2 and potentially even offer efficiency improvements in comparison to unabated plants [5]. Traditionally, the CLC process is performed in a dual circulating fluidized bed (CFB) configuration. In the oxidation reactor, a metallic oxygen carrier is oxidized, providing large amounts of heat. The thermal energy in the gas phase is used for power production, whereas the hot particles are transported to the fuel reactor where the oxidized particles are reduced by a fuel, producing CO2 and steam. The CLC process therefore keeps the CO2 stream separate from the nitrogen-containing air stream, allowing an almost pure CO2 stream for storage to be obtained simply by knocking out the water.

A drawback of the dual circulating fluidized bed CLC approach is that efficient power production with CO2 storage requires high pressure operation. However, progress on the scale-up of pressurized CLC systems has been limited [6] due to the complexity of pressurizing the two reactors, loop seals, and cyclones, and in maintaining the required solids circulation between the reactors. Consequently, several alternative CLC configurations have been proposed to overcome the challenges of the pressurized dual CFB CLC system. These include gas switching technologies [7,8], rotating bed reactors [9,10], packed bed chemical looping [11,12], and internally circulating reactors (ICRs) [13–15], which will be the focus of the present study.

The internally circulating reactor concept replaces the loop seals and cyclones that separate reactors in the dual CFB with simple ports connecting two sections of a reactor vessel, allowing the oxygen carrier to circulate between the reducing and oxidizing sections. This allows the CLC process to take place within a single unit, significantly simplifying pressurization and scale-up. The disadvantage is that gas will leak through the ports along with the circulating solids, reducing the CO2 capture efficiency and the purity of the captured CO2. However, it has been shown that the detrimental effect of gas leakage can be limited by controlling the fluidization velocity ratio of the two sections and the bed loading [13], achieving CO2 capture efficiencies greater than 95% and purities greater than 92%.

Academia has been prolific in proposing novel processes and reactors for CO2 capture. Unfortunately, implementation of new technologies in the process and energy industry has traditionally been slow, requiring several decades from process conception to commercial reality. The urgency of climate change will require very rapid scale-up and industrialization of these novel CO2 capture technologies, starting once governments start imposing strong policies to reduce carbon emissions (the IEA Sustainable Development Scenario assumes CO2 prices of \$63/ton and \$140/ton in 2025 and 2040, respectively [16]). This is also valid for other industries—rapid innovation and implementation of new process technologies will be necessary in a world that is increasingly environmentally and resource-constrained.

Simulation-based engineering will be an essential tool in enabling such rapid innovation by decreasing the number of costly and time-consuming experimental scale-up steps, and computational fluid dynamics (CFD) is the most suitable tool for investigating the chemical reactors common in the process and energy industries. However, although CFD has proven extremely useful in better understanding flow processes on the lab-scale, a common challenge to industrial simulation is the fact that important phenomena may occur on time- and length-scales that are several orders of magnitude smaller than those associated with the industrial processes [17]. This is especially relevant in multiphase processes and for the fluidized beds used in the ICR reactor studied here, where gas bubbles and particle clusters of length-scales in the order of ten particle diameters play an important role on the overall fluidized bed behavior. Using small enough grid cells and time steps to resolve these small-scale phenomena remains impossible for parametric studies of industrial-scale devices, even with large, modern computational clusters.

Multiscale methods are necessary to overcome this challenge—allowing the use of coarse computational grids to achieve reasonable computational times by using closures for unresolved subgrid effects to maintain acceptable accuracy. The filtered two-fluid model (fTFM) [18] is a common approach for multiscale modeling of fluidized beds. In the fTFM, the governing equations of the two-fluid model (TFM) closed by the kinetic theory of granular flow, where the solids phase is assumed to behave as a continuum and closures capture the effects of random particle collisions and translation, is spatially averaged, revealing subgrid terms that require closure. Several groups have strived to develop such closures. Most of the work has focused on the subgrid correction to the drag [19–26], which substantially reduces the drag compared to the drag law evaluated at the resolved quantities, although several other closures are necessary for fluidized bed hydrodynamics [27]. Research on closures for reactive flow has been limited. Most studies have investigated the influence of subgrid effects on the effective reaction rate of first-order solids-catalyzed reactions [28–30], where mass-transfer limitations imposed by the bubbles and clusters drastically reduce the effective reaction rates. Closures are also required for the dispersion of scalars (such as species and enthalpy) due to subgrid velocity fluctuations and for the effective interphase heat transfer rate [31,32].

The present study aimed to demonstrate how multiscale CFD simulations can be used to assist the evaluation of novel reactor concepts on an industrial scale, focusing on an internally circulating fluidized bed reactor for power production with CO2 capture. Firstly, some improvements were proposed for existing fTFM closures, improving the accuracy and simplicity of existing hydrodynamics closures [27,33] and, most importantly, proposing a generalized reactive fTFM closure. The latter is important, since existing closures [28–30] are only valid for simple first-order solids-catalyzed reaction equations. Next, an fTFM accounting for all important subgrid effects in reactive flows was used to evaluate the effect of several design and operating parameters on the ICR behavior. It can be noted that, to the best of the authors' knowledge, this is the most complete implementation—in terms of the number of subgrid effects accounted for—of a reactive fTFM to date. Then, results from the reactor simulations were combined with previously published power plant simulations by the same authors [34] to conduct an economic assessment of the ICR concept for low carbon power production from natural gas. Finally, the results are used to discuss the future of virtual prototyping of novel reactors using multiscale CFD simulations, as well as the potential of the ICR to combat climate change.

#### **2. Materials and Methods**

The present study utilizes both multiscale CFD reactor modeling and process modeling to inform the economic evaluation of power production with CO2 capture using the ICR concept. Figure 1 shows how information flows between these three parts of the study, and the subsequent sections describe each part in detail.

**Figure 1.** Information flow between different parts of the present study.

#### *2.1. Reactor Modeling*

#### 2.1.1. The Filtered Two-Fluid Model (fTFM)

The fTFM solves the spatially-averaged (or filtered) form of the governing equations for the two-fluid model closed by the kinetic theory of granular flow [35,36]. This section briefly presents the filtered governing equations, as well as the closures that are used for the subgrid terms. The interested reader may find a more complete discussion of the derivation of the filtered equations in earlier studies [33,37].

The filtered continuity equations are given below. *SR* is a source term due to the mass transfer during the reduction and oxidation of the oxygen carrier. Closures for the filtered reaction rates, which can be used to calculate the source term, are discussed in Section 2.1.2.

$$\frac{\partial}{\partial t} \left( \overline{a\_{\mathcal{S}}} \rho\_{\mathcal{S}} \right) + \nabla \cdot \left( \overline{a\_{\mathcal{S}}} \rho\_{\mathcal{S}} \, \overline{\overset{\textstyle \cdot}{\mathcal{V}}}{\mathcal{V}\_{\mathcal{S}}} \right) = S\_{R\prime} \tag{1}$$

$$\frac{\partial}{\partial t}(\overline{\alpha\_{\sf s}}\rho\_{\sf s}) + \nabla \cdot \left(\overline{\alpha\_{\sf s}}\rho\_{\sf s}\overline{\overrightarrow{\boldsymbol{\nu}}^{\sf s}}\right) = -S\_{\sf R}.\tag{2}$$

Next, the filtered momentum equations are shown in Equations (3) and (4):

$$\frac{\partial}{\partial t} \left( \rho\_{\mathcal{S}} \overline{\alpha\_{\mathcal{S}}} \overline{\stackrel{\textstyle \cdot}{\boldsymbol{\nu}}} \overline{\boldsymbol{\nu}} \right) + \nabla \cdot \left( \rho\_{\mathcal{S}} \overline{\alpha\_{\mathcal{S}}} \overline{\stackrel{\textstyle \cdot}{\boldsymbol{\nu}}} \overline{\stackrel{\textstyle \cdot}{\boldsymbol{\nu}}} \right) = -\overline{\alpha\_{\mathcal{S}}} \nabla \overline{p} - \nabla \cdot \left( \rho\_{\mathcal{S}} \overline{\alpha\_{\mathcal{S}}} \overline{\stackrel{\textstyle \cdot}{\boldsymbol{\nu}}} \overline{\stackrel{\textstyle \cdot}{\boldsymbol{\nu}}} \overline{\stackrel{\textstyle \cdot}{\boldsymbol{\nu}}} \right) + \nabla \cdot \overline{\overline{\frac{\pi}{4}}} + \overline{\alpha\_{\mathcal{S}}} \rho\_{\mathcal{S}} \overline{\stackrel{\textstyle \cdot}{\boldsymbol{\xi}}} + \overline{K\_{\mathcal{S}} \rho\_{\mathcal{S}} \overline{\stackrel{\textstyle \cdot}{\boldsymbol{\xi}}}} + \overline{K\_{\mathcal{S}} \left( \overline{\stackrel{\textstyle \cdot}{\boldsymbol{\nu}}} \stackrel{\textstyle \cdot}{\boldsymbol{\nu}} \right)} - \overline{\alpha\_{\mathcal{S}} \prime} \overline{\nabla p'}, \tag{3}$$

$$\frac{\partial}{\partial t} \left( \rho\_{\delta} \overline{\alpha\_{\delta}} \stackrel{\rightharpoonup}{\upsilon \delta} \right) + \nabla \cdot \left( \rho\_{\delta} \overline{\alpha\_{\delta}} \stackrel{\rightharpoonup}{\upsilon \delta} \stackrel{\rightharpoonup}{\upsilon \delta} \right) = -\overline{\alpha\_{\delta}} \nabla \overline{p} - \nabla \cdot \left( \rho\_{\delta} \overline{\alpha\_{\delta} \stackrel{\rightharpoonup}{\upsilon \delta} \stackrel{\rightharpoonup}{\upsilon \delta} \stackrel{\rightharpoonup}{\upsilon \delta}} \right) - \nabla \overline{p\_{\delta}} + \nabla \cdot \overline{\overline{\overline{\pi}}\_{\delta}} + \overline{\alpha\_{\delta}} \rho\_{\delta} \stackrel{\rightharpoonup}{\mathcal{G}} + \overline{K\_{\mathcal{G}} \left( \overline{\upsilon \mathcal{I}}\_{\mathcal{X}} - \overline{\upsilon}\_{\mathcal{I}} \right)} - \overline{\alpha\_{\mathcal{I}} \! \!/ \!/ \!/ \!} \overline{p\_{\mathcal{I}}}.\tag{4}$$

In both filtered momentum equations, the second term on the right-hand side represents the stresses due to subgrid velocity fluctuations (arising from gas bubbles and solids clusters), which add to the diffusive momentum transport. The gas-phase subgrid stress is usually relatively small compared to the solids stresses due to the large difference in the phase densities [38] and can safely be neglected. However, the solids subgrid stresses are accounted for by means of an anisotropic stress closure [39], which has been shown to offer significant improvements compared to Boussinesq approximation-based closures using isotropic independent variables [39,40]. In the filtered solids momentum equation, the filtered kinetic theory stresses (third and fourth terms on the right-hand side) are small at the grid sizes that are relevant for industrial-scale fluidized beds [39], which was used in the present study. Therefore, the filtered kinetic theory stresses were estimated on the basis of the unfiltered granular temperature equation, as it was previously shown to be sufficient [41].

The second-to-last term on the right-hand side of both momentum equations represents the filtered drag force, where subgrid effects generally reduce the drag compared to that in a homogenous suspension. This is due to the tendency of fluidized particles to form solids clusters and gas bubbles, which are not resolved on a coarse grid. These meso-scale structures vary in size and shape due to local flow conditions. Gas will tend to pass through dilute regions, reducing the effective drag on the solids clusters, the effect of which must be accounted for in a closure. A modified version of a 3-marker anisotropic closure published previously [27] was used to close the filtered drag force. It was found that the 3-marker closure could be simplified significantly, while maintaining similar accuracy, by eliminating the filtered slip velocity as a marker. More information about the development and verification of this new closure can be found in the Supplementary Material.

Finally, the last term on the right-hand side of both momentum equations is due to subgrid pressure gradient fluctuations and is referred to here as the meso-scale interphase force. This contribution arises from the redistribution of the pressure gradient over subgrid gas bubbles and solids clusters [42] and tends to add to the effective drag force [33]. For the present study, an older anisotropic closure [33] was improved on by drawing an analogy to the closure for the meso-scale solids stresses [39], where it was found that a filtered co-variance term can be accurately closed as a function of the relevant gradients. The Supplementary Material also details the development and verification of this new closure.

Next, Equation (5) gives the filtered species transport equation for reactant A, which is consumed in an nth order reaction.

$$\frac{\partial}{\partial t}(\rho\_{\mathcal{S}}\overline{\alpha\_{\mathcal{S}}}\widetilde{X\_{A}}) + \nabla \cdot \left(\rho\_{\mathcal{S}}\overline{\alpha\_{\mathcal{S}}}\widetilde{X\_{A}}\overline{\overline{\nu\_{\mathcal{S}}}}\right) = \nabla \cdot \left(\rho\_{\mathcal{S}}\overline{D a\_{\mathcal{S}}\nabla X\_{A}}\right) - \nabla \cdot \left(\rho\_{\mathcal{S}}\overline{\alpha\_{\mathcal{S}}X\_{A}\prime\prime\prime\_{\mathcal{S}}}\right) - \overline{k\_{A}a\_{\mathcal{S}}C\_{A}^{\
u}M\_{A}}.\tag{5}$$

*Processes* **2019**, *7*, 723

The species dispersion due to the filtered microscopic diffusion (the first term on the right-hand side) was expected to be small relative to meso-scale dispersion, as well as convective transport. Therefore, in line with previous work regarding scalar dispersion in fTFMs [31], it was simply evaluated as ∇ · ρ*gD*α*g*∇*X <sup>A</sup>* in the present study. The subgrid species dispersion rate (the second term on the right-hand side) tends to disperse the species due to sub-grid velocities arising from unresolved gas bubbles and solids clusters. This effect was accounted for using the closures of Agrawal et al. [31] but has been shown to only have a minor effect on the overall reaction rate [28]. The filtered reaction rate (third term on the right-hand side) is typically substantially reduced by subgrid bubbles and clusters and is essential to model [28]. This is because, for gas-solid reactions, the reactant will be consumed faster inside dense regions, creating a mass transfer limitation due to the finite rate at which reactants are transported to these dense regions. A limited number of studies have investigated reactive fTFM closures [28,30,41], but they have all focused on reactions that are solids catalyzed and first-order with respect to the gaseous reactant. The next section of the present study therefore proposes a novel, simplified approach for accounting for different reaction orders and for reactions where the solids phase participate in the reaction. Finally, it can be noted that the filtered solids species equations are similar to those of the gas-phase species and are thus treated in a similar way. Consequently, they are not discussed separately.

Finally, Equation (6) shows the filtered enthalpy transport equation for the gas-phase.

$$\frac{1}{2}\frac{\partial}{\partial t}\Big(\rho\_{\mathcal{S}}\overline{a\_{\mathcal{S}}}\overline{h\_{\mathcal{S}}}\Big) + \nabla\cdot\Big(\rho\_{\mathcal{S}}\overline{a\_{\mathcal{S}}}\overline{h\_{\mathcal{S}}}\overline{\widehat{\boldsymbol{\nu}}\_{\mathcal{S}}}\Big) = \nabla\cdot\Big(\overline{\kappa\_{\mathcal{S}}a\_{\mathcal{S}}\nabla T\_{\mathcal{S}}}\Big) - \nabla\cdot\Big(\rho\_{\mathcal{S}}\overline{a\_{\mathcal{S}}h\_{\mathcal{S}}\mathcal{T}}\overline{\widehat{\boldsymbol{\nu}}\_{\mathcal{S}}\prime}\Big) + \overline{\gamma\Big(T\_{\mathcal{S}}-T\_{\mathcal{S}}\Big)} + \overline{k\_{A}a\_{\mathcal{S}}C\_{A}^{\eta}M\_{\mathcal{A}}\Delta H\_{r,\mathcal{A}}}.\tag{6}$$

Here, as with the species diffusion, the filtered microscopic conductivity (first term on the right-hand side) is small compared to the enthalpy dispersion from sub-grid velocity fluctuations and is simply approximated as κ*g*α*g*∇*T <sup>g</sup>*. The subgrid enthalpy dispersion rate (second term on the right-hand side) and the filtered heat transfer rate (third term on the right-hand side) were modeled using the closures of Agrawal et al. [31]. The physical behavior of these contributions is analogous to that of the species dispersion rate and filtered reaction rate, due to the similarity of mass and heat transfer. The enthalpy source term due to reaction (fourth term on the right-hand side) was evaluated at the filtered reaction rate modeled in Equation (5), assuming that the heat of reaction is uniform for each cell in the coarse-grid simulations. This is a reasonable assumption based on the good mixing and fast heat transfer in fluidized beds. Finally, it can be noted that the solids-phase filtered enthalpy equation was treated in a similar way, and it is therefore not discussed here separately.

#### 2.1.2. Reaction Modeling

In the present study, Ni/NiO (supported on Al2O3) was used as oxygen carrier due to its high reactivity [43] and its ability to tolerate high operating temperatures [44]. In the fuel section, the oxygen carrier was reduced by the fuel according to:

$$4\text{NiO}(\text{s}) + \text{CH}\_4(\text{g}) \rightarrow 4\text{Ni}(\text{s}) + 2\text{H}\_2\text{O}(\text{g}) + \text{CO}\_2(\text{g}).\tag{7}$$

In the air side of the ICR, the oxygen carrier re-oxidizes as:

$$2\text{Ni}(\text{s}) + \text{O}\_2(\text{g}) \rightarrow 2\text{NiO}(\text{s}).\tag{8}$$

The reactions are implemented in the filtered species conservation equations, as follows, for species *i* taking part in reaction *k*, where *vi* is the stoichiometric constant:

$$
\frac{\partial}{\partial t} \left( \rho\_{\mathcal{S}} \overline{a\_{\mathcal{S}}} \overline{X\_{i}} \right) + \nabla \cdot \left( \rho\_{\mathcal{S}} \overline{a\_{\mathcal{S}}} \overline{X\_{i}} \overline{\overset{\leftarrow}{\boldsymbol{\nu}}\_{\mathcal{S}}} \right) = \nabla \cdot \left( \rho\_{\mathcal{S}} D \overline{a\_{\mathcal{S}}} \nabla \overline{X\_{i}} \right) - \nabla \cdot \left( \rho\_{\mathcal{S}} \overline{a\_{\mathcal{S}} X\_{i} \prime \prime} \overline{\overset{\leftarrow}{\boldsymbol{\nu}}\_{\mathcal{S}}} \right) + \sum\_{k=0}^{n} v\_{\mathcal{S}} R\_{k}^{H} M\_{i}.\tag{9}
$$

The effective reaction rate, *RH <sup>k</sup>* , in units of *mol*/ *m*3*s* was calculated as shown in Equation (10), where gas species A reacts with solids species B (here, B represents NiO in reduction and Ni in oxidation).

$$R\_{\mathbf{k}}^{H} = \eta \overline{\alpha\_{s}} \rho\_{s} a\_{0} \widehat{(\mathcal{X}\_{\text{Ni}} + \hat{X}\_{\text{NiO}})} k\_{s} \widehat{\mathcal{C}\_{A}}^{n} \left( 1 - \frac{\widehat{X\_{B}}}{\widehat{X\_{\text{Ni}}} + \widehat{\overline{X\_{\text{Ni}}}}} \right)^{\frac{2}{3}}.\tag{10}$$

The solid particles are porous, and the reaction can be considered to be kinetically controlled following the shrinking core model applied to microscopic grains inside the porous particle [45]. Application of the shrinking core model with reaction rate control [46] is evident in the final factor of Equation (10). The reaction rate constant, *ks*, is expressed as follows, where the detrimental effect of increasing pressure is accounted for in the pre-exponential factor:

$$k\_s = \frac{k\_0}{\overline{p}^q} e^{-\frac{E\_0}{kT}}.\tag{11}$$

The kinetic parameters for the reduction and oxidation reactions, as well as the oxygen carrier properties, were obtained from the experimental work of Abad et al. [45]. It can be noted that the aforementioned study found no intraparticle mass transfer limitations, as may be expected for such small, porous particles.

In the fTFM, the subgrid bubbles and clusters impose an additional mass transfer limitation on the reactions, since the gaseous reactants have to be transported into the dense solid clusters for the reactions to occur. This effect is modeled in Equation (10) by means of an effectiveness factor, η. In the present study, η was first modeled for a reference first-order solids catalyzed reaction with a fixed reaction rate constant, as in previous studies [28,41]. It was then found that the reference closure can be effectively scaled to different reaction rate constants and reaction orders by drawing an analogy with packed bed theory and defining a cluster-scale Thiele modulus, as follows:

$$
\phi = \sqrt{\frac{n+1}{2} \frac{k' \left(d\_p \mathcal{L}\right)^2}{D}}.\tag{12}
$$

Here, L is the average ratio of the cluster diameter to the particle diameter, which requires closure, and φ is the Thiele modulus [47]. The effective reaction rate constant, *k* , was obtained by re-writing the reaction equation as first-order with respect to the gaseous reactant and the solids volume fraction. This approach has previously been shown to be useful to extend effectiveness factors from intraparticle mass transfer theory to various reaction orders [48]. For the example of Equation (10), the effective reaction rate constant becomes:

$$k'=\rho\_5 a\_0 \widetilde{(X\_{Ni}} + \widehat{X\_{NiO}}) k\_s \widetilde{C\_A}^{n-1} \left(1 - \frac{\widetilde{X\_B}}{\widetilde{X\_{Ni}} + \widetilde{X\_{NiO}}}\right)^{\frac{2}{5}}.\tag{13}$$

The effectiveness factor for a spherical particle can then be written as follows [49]:

$$
\eta = \frac{1}{\phi} \left( \frac{1}{\tanh(3\phi)} - \frac{1}{3\phi} \right). \tag{14}
$$

This relation is exact for a first-order reaction in a porous particle with no convective transport. Relatively small discrepancies arise for reactions of different order, but the largest uncertainty in this application is the constant deformation of the clusters in the fluidized bed, as well as the convective species transport taking place inside the cluster.

The basic premise of the approach proposed in the present study is that the effectiveness factor in Equation (14) is analogous to the effectiveness factor of a particle cluster at the largest achievable mass transfer resistance (smallest effectiveness factor). This will typically occur at intermediate filtered solids volume fractions when maximum phase segregation is achieved and clusters are relatively large. As the filtered solids volume fraction tends to the limits of zero or maximum packing, clustering disappears, and the mass transfer resistance tends to zero (η = 1).

A hypothesis can then be formulated that, for different cluster-scale Thiele moduli, the minimum effectiveness factor (η*min*) can be scaled by using Equation (15) when a filtered effectiveness factor closure (η*re f*) is derived from resolved simulation data for a reference Thiele modulus (φ*re f*).

$$\eta\_{\text{new},\text{min}} = \eta\_{\text{ref},\text{min}} \frac{\frac{1}{\phi\_{\text{new}}} \left( \frac{1}{\tanh(3\phi\_{\text{new}})} - \frac{1}{3\phi\_{\text{new}}} \right)}{\frac{1}{\phi\_{\text{ref}}} \left( \frac{1}{\tanh\left(3\phi\_{\text{ref}}\right)} - \frac{1}{3\phi\_{\text{ref}}} \right)}. \tag{15}$$

Then, the new effectiveness factor (η*new*) can be calculated as follows, assuming that the tendency towards η = 1 will be proportional to the tendency of η*min* to unity (no subgrid correction):

$$\eta\_{\text{new}} = 1 - \left(1 - \eta\_{\text{ref}}\right) \frac{(1 - \eta\_{\text{new},\text{min}})}{\left(1 - \eta\_{\text{ref},\text{min}}\right)} \,. \tag{16}$$

It was found that the suggested hypothesis holds well and that this approach is essential to accurately model reactions that are not simple first-order solids catalyzed reactions in the fTFM. Consequently, the proposed approach was used to model the reactions in the present study. The complete development and verification of the generalized reactive fTFM closure is presented in the Appendix A.

Finally, it can be noted that the effectiveness factor closure presented here does not account for the Stefan flow (one mole of methane produces three moles of gas products) occurring in the fuel section of the ICR and investigation of this topic is recommended for future work. However, considering that the reduction reactions are extremely fast (see Section 3.1.1) and occur only near the inlet, it is not expected to have a large impact on the overall reactor behavior.

#### 2.1.3. Simulation Geometry and Mesh

Figure 2 shows the reactor geometry that was considered for the ICR. In the base case, the reactor consists of a cylinder with a height of 6.92 m and a diameter of 3.46 m. These sizes were selected to yield a fluidization velocity of roughly 1 m/s in the freeboard, which is a typical value for vigorous bubbling fluidization. An aspect ratio of 2, typical of fluidized beds, was chosen. A thin wall separates the reactor into the reduction and oxidation sections, consisting ofa2m high vertical section at the center of the bed and a section sloping at an angle of 30◦ with the vertical axis to the reactor wall to ensure that solids will not deposit on this surface. Two ports allow the oxygen carrier to circulate between the sections. Reduced oxygen carrier travels through the bottom port (height of 0.6 m) to the air section, whereas oxidized oxygen carrier is carried through the top port (terminating 0.7 m above the bottom of the reactor) to the reduction section. In the base geometry, the width of the square ports (see Figure 2b) is 20 cm. The gas outlet from each section was sized to yield a velocity of roughly 50 m/s, accounting for the much larger flow rate in the air section, which is a typical value for gas transport.

A cut-cell mesh was used to mesh the complex ICR geometry. Long simulations, in the order of 2500 simulated seconds, were necessary to achieve steady reactor behavior; therefore, the average grid size was chosen to yield a coarse mesh of approximately 50,000 cells. A minimum of five cells across the gaps in the ports and the outlets were specified to resolve the most important flow gradients.

**Figure 2.** Simulated geometry for the internally circulating reactor (ICR): (**a**) isometric view; (**b**) details of the ports connecting the reactor sections, where x is the port width; (**c**) front view; (**d**) side view.

#### 2.1.4. Reactor Operating Conditions

NiO particles supported on Al2O3 were used as oxygen carrier. The oxygen carrier particles were considered to have a diameter of 150 μm, a typical value for bubbling fluidization, a density of 3446 kg/m3, and an active mass fraction of 0.4 [45]. It can be noted that the reactor model assumes monodisperse particles, due to the complexity of accounting for particle size distributions in fTFMs and due to the limited state of development of subgrid closures accounting for polydispersity [50]. Additionally, the simulations assume the particle density to be constant during the reactor operation,

since the density changes will be small (mostly less than 5%) due to the high inert content of the oxygen carrier. The loading of the bed corresponds to an initial bed height of 1 m at a solids volume fraction of 0.6.

Fluidization gas was added uniformly at the bottom of the reactor to the two reactor sections, assuming a perfect gas distributor. Additionally, the effect of the inlet conditions on the subgrid behavior were not accounted for in the fTFM closures. These simplifications are necessary, since none of the state-of-the-art fTFMs have thus far included these effects in their closures. However, considering the large dimensions of the reactors considered, inlet effects are expected to have a relatively small influence on the overall reactor behavior, thereby minimizing the error associated with these simplifications. The other inlet boundary conditions are listed in Table 1 (note that the natural gas used in the process simulations was replaced with an equivalent amount of methane in the reactor simulations).

**Table 1.** Summary of the conditions for the inlets of the two reactor sections.


Uniform pressure outlet boundary conditions were considered for the fuel and air section outlets. For the air section outlet, a pressure of 18 bar (absolute) was considered, which results from the air compressor pressure ratio of 18 [51] employed in the process simulations, which is a typical value for standard, large-scale, F-class gas turbines [52]. For the fuel section outlet, a relatively small overpressure relative to the air section outlet was employed to achieve a target flow rate. This is discussed in more detail in Section 3.1.2.

A no-slip boundary condition was specified for the gas at the walls, whereas partial slip boundary conditions with a specularity coefficient of 0.1 was employed at the walls, based on the model of Johnson and Jackson [53]. It can be noted that, technically, a subgrid closure is required for the particle–wall interaction. However, such closures have not yet been developed in the fTFM community and were therefore neglected in the present study. The effect of the particle–wall boundary condition is expected to be small for the large reactor dimensions considered in the present study; therefore, neglecting the sub-grid effects is a reasonable assumption.

#### 2.1.5. Solver

The reactor simulations were performed in the commercial CFD solver, ANSYS FLUENT 19.2, using user defined functions to implement the subgrid closures of the fTFM. The phase-coupled SIMPLE algorithm [54] was used for pressure-velocity coupling, and all other equations were discretized based on the QUICK scheme [55].

#### *2.2. Process Modeling*

This study conducted an economic assessment of the ICR integrated into a natural gas combined cycle (NGCC) power plant, as recently evaluated for CO2 capture using CLC [34]. The interested reader is referred to that study for details about the process modeling methodology. One important change from this previous work was the inclusion of the gas leakage between reactor sections in the ICR. The mixing between the outlet streams of the oxidation and reduction reactors was adjusted to yield 95% CO2 capture and purity (molar percentage and dry basis), based on the reactor simulations (see Section 3.1.3) for the conditions considered.

The layout of the simulated plant is shown in Figure 3, where 20 parallel ICR reactors were needed to accommodate the required air throughput. Natural gas is pre-heated and fed to the fuel section of the ICR reactors where it is converted to CO2 and H2O, which is expanded to generate some power. After H2O is condensed out, the remaining CO2 is compressed and pumped to 110 bar for transport and storage. The air section of the ICR replaces the combustor for the main gas turbine. Air from the main compressor reacts with the reduced oxygen carrier in a highly exothermic reaction and is heated to 1150 ◦C in the base case. This temperature was selected based on material limitations, and a sensitivity analysis of this value was performed in the economic evaluation in Section 3.2. The hot depleted air stream is then expanded in the main gas turbine before being sent to a heat recovery steam generator for extra power production using a steam cycle. The results of this plant were compared to the reference NGCC plant detailed in Khan et al. [34].

**Figure 3.** Process flowsheet of the ICR integrated into a combined cycle. It can be noted that the reactor in the flowsheet represents a cluster of ICR reactors and that the outputs from these reactors were combined in stream 3.

#### *2.3. Economic Assessment*

Capital costs: The total cost of the combined power cycle was obtained directly from the PEACE component in Thermoflex. This includes direct component costs and several additional cost components accounting for construction, engineering, contingencies, and other cost components.

Costs related to the CO2 compressors and intercoolers were estimated using installed cost data from Aspen Plus. This cost was increased by approximately 74% to account for engineering, contingencies, and owner's costs, based on the methodology of Gerdes et al. [56].

The ICR capital costs were estimated based on cost correlations for process vessels from Turton et al. [57]. Each ICR was composed of two process vessels: (1) an inner vessel to carry the temperature, attrition, and corrosive loads constructed from an expensive Ni-alloy, and (2) a thick pressure shell carrying the pressure load constructed from carbon steel. An insulation layer of 0.4 m thickness was inserted between these two vessels. To account for the relatively complex ICR geometry, the cost of the inner vessel was increased by a factor of three. This was a somewhat arbitrary adjustment, and a sensitivity of total plant economics to ICR cost is therefore presented later. Costs for auxiliaries and contingencies were subsequently added according to Turton et al. [57] to yield the total reactor cost. A breakdown of the different components of the cost of the 20 ICR units required for the base case is shown in Figure 4. It can be noted that the capital cost associated with the oxygen carrier was only for the initial loading (replacement of the oxygen carrier was considered under operating and maintenance costs). Further, the number of ICR reactors selected to deliver the required process throughput at the reactor operating conditions is specified in Section 2.1.4.

**Figure 4.** Breakdown of the cost of 20 ICR units.

The total costs of these three main components were then added together to yield the total plant cost. All costs are reported in 2019 US dollars using the chemical engineering plant cost index (CEPCI).

Operating and maintenance (O&M) costs: Fixed O&M labor costs were calculated assuming 11 personnel per shift for the NGCC reference plant and 13 personnel per shift for the ICR plant using the methodology of Peters and Timmerhaus [58]. A \$45/hour rate was used with appropriate increases for benefits, maintenance labor, and overheads. In addition, 1.5% of total plant costs per year was added for insurance and property taxes. The key assumptions for variable O&M costs are summarized in Table 2. Costs for the oxygen carrier [59] and water [60] were taken from the literature, whereas the oxygen carrier lifetime was specified based on discussions with catalyst suppliers. CO2 transport and storage costs can vary widely based on the transport distance and the type of transport and storage, but a reasonable average value was selected based on costs provided in two IEA reports [61,62]. Natural gas prices are known to vary widely, and a value representative of Europe was assumed for this study.

**Table 2.** Variable O&M cost assumptions.


Capital and O&M costs were then used to calculate the levelized cost of electricity (Equation (17)) and the CO2 avoidance costs (Equation (18)) using a discount rate of 8%, a plant economic lifetime of 30 years, and a construction period of 2 years for NGCC and 3 years for the ICR plants (investment is assumed to be linear over the construction period).

$$\text{LCOE} \left( \\$/\text{MWh} \right) = \frac{\sum\_{t=1}^{m} \frac{I\_t + M\_t + F\_t}{\left(1 + r\right)^t}}{\sum\_{t=1}^{m} \frac{E\_t}{\left(1 + r\right)^t}} \text{.} \tag{17}$$

$$\text{CAC} \left( \\$/\text{ton}\right) = \frac{\text{LCOE}\_{\text{CCS}} - \text{LCOE}\_{\text{ref}}}{\varepsilon\_{\text{nf}} - \varepsilon\_{\text{CCS}}}. \tag{18}$$

Here, the summations are done for each year during construction and operation (*t*) up to the end of the plant economic lifetime (*m*). *I* is the investment expenditures, *M* is the O&M expenditures, *F* is the fuel expenditure, *E* is the electricity generation, *r* is the discount rate, and *e* is the plant-specific emissions (ton/MWh).

#### **3. Results**

This section outlines how the fTFM described in Section 2.1 was first used to optimize and size an industrial-scale ICR reactor for power production with integrated CO2 capture. Subsequently, the ICR process was then evaluated economically using the reactor size and performance suggested by the simulations.

#### *3.1. Reactor Optimization*

#### 3.1.1. Characteristics of ICR Operation

In this study, plots and animations from the reactor simulations were used to introduce important characteristics of ICR operation for a typical case. Firstly, Figure 5 (as well as the associated Video S1 in the Supplementary Material) demonstrates the circulation of the oxygen carrier between the two reactor sections. In the reduction section, the relatively low fluidization velocity from the fuel feed results in a dense bubbling bed. Due to the very fast reaction of the oxygen carrier with the methane, most of the conversion takes place near the inlet, leading to a highly reduced oxygen carrier. However, the mixing is very fast in the fluidized bed and a relatively uniform distribution of the oxygen carrier is rapidly attained in the rest of the bed on the fuel side. Oxygen carrier particles, reduced by the fuel, pass through the bottom port to the air (oxidation) section, where they are rapidly oxidized and mixed into the rest of the particles. Owing to the much larger molar flow rate on the air side, a more vigorous fluidization occurs, lifting the particles to the freeboard, including a diameter expansion (which helps to reduce particle elutriation), and allowing them to pass back to the fuel section through the top port. Again, the oxidized particles mix rapidly into the reduction side bed, where they are reduced by the fuel.

**Figure 5.** Particle plot of the instantaneous NiO mass fraction. A corresponding animation is provided in Video S1. Note that the particles pictured are tracers following the continuous solids flow for visualization purposes and do not influence the simulation solution.

Some further comments can be made about the nature of the solids flow through the ports. In the bottom port, the flow is quite dense, with a time-average solids volume fraction of about 0.5. The flow through the top port is more dilute, with a time-averaged solids volume fraction between 0.3 and 0.4. The animations show some transient fluctuations of solids in both ports; therefore, the flow in not completely steady and there is a risk of backflow, which might reduce the reactor performance. It may be noted that no problems with blockage of the ports have been experienced during extensive experimental evaluations of the ICR concept [14,15].

Figure 6 (and Video S2) shows that the solids circulation between the reactor sections is associated with undesired gas leakage—CO2 leaks from the fuel section to the air section, reducing the CO2 capture efficiency of the reactor, and N2 leaks from the air section to the fuel section, reducing the purity of the CO2. One of the most important criteria for designing and operating the ICR is therefore to minimize the amount of gas leakage between the reactor sections, while maintaining sufficient oxygen carrier circulation to ensure that the fuel is completely converted in the reduction sector.

**Figure 6.** Contour plot of the instantaneous CO2 and N2 mole fractions at the outer wall of the ICR showing the undesired gas leakage between the two sections of the reactor. A corresponding animation is provided in Video S2.

Many design and operating parameters can influence the ICR performance. These include, but are not limited to, the solids loading, the particle size and distribution, the gas flow rates to the reactor sections, the operating pressure and temperature, and several dimensions of the reactor and internals. Due to the complexity of simultaneously optimizing these parameters, the scope of the present study was limited to three important factors that will be investigated in the subsequent sections:


#### 3.1.2. Reduction Section Overpressure

During ICR operation, the most practical way to control the reactor behavior would be to tune the pressure difference between the reduction and oxidation section outlets. For example, applying an overpressure at the reduction side outlet will lead to more gas exiting the reactor on the oxidation side. To satisfy the mass balance of the reactor, this means that relatively more gas must pass through the bottom port to reach the oxidation section. As is shown in this section, this gas flow influences the solids circulation between the sections, as well as the distribution of solids in the two sections. Both these factors have a critical effect on the reactor performance.

To better understand this behavior, ICR simulations were performed at different overpressures applied at the reduction section outlet. Specifically, reduction outlet flow ratios (ROFR) from 0.96 to 1 were investigated. The ROFR is defined as the ratio of the reduction outlet molar flow rate to the ideal outlet molar flow rate that would occur in case of no gas leakage between the reactor sections and complete fuel conversion. A lower ROFR implies a higher overpressure in the fuel section, with the ROFR = 0.96 case corresponding to an overpressure of 0.34 bar. It can also be noted that the figures and animations presented in the previous section are for the case of ROFR = 0.98.

Figure 7 (and Video S3) shows the effect of the reduction section outlet overpressure on the ICR behavior. Firstly, it is interesting to note from the solids volume fraction values that, despite the large grid sizes employed, a substantial amount of phase segregation is still resolved in the more vigorously fluidized air section. In the slowly fluidized fuel section, the resolved solids distribution is nearly homogenous, and the effects of particle clusters and bubbles are therefore nearly completely accounted for in the subgrid closures of the fTFM.

Increasing the overpressure (corresponding to lower ROFR values), more gas will pass through the bottom port to the air section, which also increases the solids flow rate through the bottom port. This causes the bed loading on the oxidation side to increase, which creates the hydrostatic pressure buildup required to achieve solids flow through the top port against the overpressure imposed in the fuel section.

Therefore, cases with a lower ROFR will reach a pseudo-steady state (where the time-averaged solids flow rates through the top and bottom ports are equal) with a larger fraction of the oxygen carrier on the oxidation side, as shown in Figure 8b. The solids elutriation rate (Figure 8a), which occurs almost entirely from the oxidation side, is therefore greatest at lower ROFR values.

Furthermore, Figure 8a shows that a maximum solids circulation rate occurs at an ROFR of 0.98. At first, when decreasing the ROFR from 1, the solids circulation rate increases due to more gas flow through the bottom port, thereby entraining more solids, as well as more solids flowing through the top port as a result of the higher bed loading on the air side. However, as the ROFR further decreases, the bed on the fuel side becomes so low that the solids barely reach the top of the bottom port, thus limiting the achievable circulation rate through the bottom port. Regular backflow through the bottom port also starts to happen during these cases, since the low bed height on the fuel side does not provide sufficient hydrostatic pressure to maintain a steady flow through the bottom port. If the solids circulation rate becomes too low, significant fuel slip, that is, incomplete fuel conversion, starts to occur (Figure 8c), since not enough oxygen is available in the fuel section to convert all of the fuel and the low bed height reduces the gas residence time in the bed.

**Figure 7.** Instantaneous particle plot of the cases with varying reduction outlet flow ratios, where the particles are colored by the particle volume fraction. A corresponding animation is provided in Video S3. Note that the particles pictured are tracers following the continuous solids flow for visualization purposes and do not influence the simulation solution.

Figure 8c also shows that the CO2 capture percentage decreases as the ROFR decreases, due to the increased gas flow through the bottom port, allowing more CO2 to exit with the depleted air at the oxidation section outlet. The CO2 purity behavior is more complex—it remains relatively constant when lowering the ROFR from 1 to 0.98, despite the increasing solids flow rate through the top port. This indicates that the gas-to-solids leakage ratio through the top port increases at high ROFR values, resulting in more gas leakage per unit of solids circulation. The purity increases in the ROFR = 0.97 case due to the lowering solids flow rate but decreases in the ROFR = 0.96 case due to backflow through the bottom port.

Based on the results in this section, the ROFR = 0.98 case (corresponding to a reduction outlet overpressure of 0.26 bar) was chosen for further investigation, primarily due to the high solids circulation rate that was obtained. Furthermore, this case showed a good compromise between decreasing solids elutriation and decreasing CO2 capture with decreasing ROFR.

**Figure 8.** Summary of important time-averaged properties of the ICR as a function of the reduction outlet flow ratio: (**a**) solids circulation and elutriation rates; (**b**) fraction of solids and NiO mass fraction in the reduction side; (**c**) CO2 capture, CO2 purity (molar percentage on a dry basis), and fuel conversion.

#### 3.1.3. Port Size

The previous section revealed that the primary criterion for achieving complete fuel conversion is a sufficient solids circulation rate to transport enough oxygen to the fuel reactor to oxidize the methane. However, since the ROFR = 0.98 case from the previous section had a much higher than necessary solids circulation rate, there is the potential to further increase the reactor performance by decreasing the size of the ports (dimension *x* in Figure 2). This will decrease the solids circulation rate between the reactor sections, but also decrease the associated undesired gas leakage.

Simulations were therefore performed at an ROFR of 0.98 while decreasing the port size dimensions, as shown in Figure 9. As expected, the results showed a decreasing solids circulation rate with decreasing port size. Consequently, the average NiO mass fraction in the fuel section was also reduced. In the

case with a port size of 16 cm, not enough NiO was present in the reduction section to fully oxidize the fuel, leading to a large fuel slip, which would be unacceptable in practice.

**Figure 9.** Plot of the solids circulation rate and the NiO mass fraction in the fuel section with changing port size.

The case with ROFR = 0.98 and a port size of 18 cm was consequently chosen as the best ICR case. Compared to the case with ROFR = 0.98 and a port size of 20 cm, the CO2 capture efficiency increased from 95.0% to 95.7%, and the CO2 purity increased from 93.1% to 95.2% while still preserving complete fuel conversion. Consequently, CO2 capture efficiencies and purities of 95% are assumed as reasonable for the economic assessment in Section 3.2.

#### 3.1.4. Reactor Size

One of the main criteria for reactor sizing is usually to ensure complete reactant conversion, with the increasing reactor diameter and height both serving to increase the gas residence time (the latter by decreasing the superficial velocity of a fixed mass flow rate of gas) for complete conversion. However, as previously noted, the oxidation reaction of methane with NiO is very fast for the temperatures considered for ICR operation. Additionally, the height of the bed on the fuel side is limited to a minimum due to the height of the bottom port. Consequently, the solids circulation rate between the reactor sections is primarily responsible for ensuring complete fuel conversion, and the overall reactor size is of lesser importance.

Nonetheless, the overall reactor size has an important effect on the particle elutriation. In the best case investigated so far (ROFR = 0.98 and a port size of 18 cm), the predicted solids elutriation rate was 11 kg/s, corresponding to 2.0% of the total solids loading per hour. Consequently, a rather large amount of elutriated solids will have to be separated using cyclones and filters downstream of the reactor. The amount of solids elutriation can be reduced by increasing the reactor height (more space for solids to fall back down) and diameter (lower superficial gas velocities).

To investigate the effect of increasing the reactor size, all the reactor dimensions were multiplied by a scaling factor (SF) from 1 to 1.2 (compared to the base case geometry discussed in Section 2.1.3. with a port size of 18 cm). Figure 10 (and Video S4) visualize the change in reactor behavior with increasing scaling factor, showing clearly that a larger reactor leads to a slower fluidization and less solids elutriation on the air side.

**Figure 10.** Instantaneous particle plot of the cases with reactor scaling factors ranging from 1 to 1.2, where the particles are colored by the particle volume fraction. A corresponding animation is provided in Video S4. Note that the particles pictured are tracers following the continuous solids flow for visualization purposes and do not influence the simulation solution.

The decreasing elutriation rate is quantified in Figure 11a, with the solids elutriation (in percentage of solids loading per hour) decreasing from 2.0% to 0.49% at a scaling factor of 1.1, and 0.20% at a scaling factor of 1.2. However, it was also found that the CO2 capture and purity decreases with increasing scaling factor (Figure 11c). This is because the port sizes were also scaled along with the other reactor dimensions, leading to an increase in the solids circulation rate (Figure 11a) and the associated gas leakage. However, the results also showed an increase in the average NiO fraction in the fuel section (Figure 11b), indicating higher than necessary solids circulation rate. As a result, it is expected that CO2 capture efficiencies and purities of more than 95% can be obtained by again decreasing the port size to achieve a solids circulation rate close to the minimum required for complete fuel conversion, as was done in Section 3.1.3.

In summary, the preferred ICR reactor size will in practice be determined by the trade-off in increasing capital cost of the reactor and the decreasing cost of handling elutriated solids with increasing reactor size. The cost of the former is quantified in the subsequent economic assessment, whereas quantification of the latter is left for a future study.

**Figure 11.** Summary of important time-averaged properties of the ICR as a function of the reactor size scaling factor: (**a**) solids circulation and elutriation rates; (**b**) fraction of solids and NiO mass fraction in the reduction side; (**c**) CO2 capture, CO2 purity, and fuel conversion.

#### *3.2. Economic Evaluation*

In this section, an ICR-based NGCC plant is evaluated economically based on the CO2 capture and purity estimated in the reactor simulations, as well as the size of reactor that was found to be adequate (SF = 1). In addition, an unabated NGCC plant is assessed as a reference.

Three ICR cases were investigated with different assumptions about the maximum reactor temperature. In combined cycles, it is crucial to maximize the turbine inlet temperature to maximize net electric efficiency and minimize the required air flow rate through the system. However, the reactor, oxygen carrier, and downstream filters will impose constraints on the maximum allowable temperature. The oxygen carrier material may be the most important constraint. To date, the highest successfully demonstrated operating temperature has been a NiO-based oxygen carrier successfully operated at 1185 ◦C [44]. However, modern gas turbines can operate with combustor outlet temperatures well over 1600 ◦C, implying that this is a very important constraint. The economic impact of this constraint will be illustrated in this work.

To this end, the cases listed in Table 3 were evaluated in this study. It can be noted that the reactor simulations were performed only for the ICR-1150 case, since it was expected that the gas leakage would not vary much between the cases at different temperatures.

It is immediately evident from the results in Table 3 that the achieved plant efficiency is strongly dependent on the ICR outlet temperature. When the temperature can reach up to 1300 ◦C, an attractively low energy penalty of 3.6% points is attained, but this increases to an unacceptable 12.3% points in the case with an ICR outlet temperature of 1000 ◦C. Since fuel costs are typically the major cost component in a natural gas fired plant, this is expected to have a major impact on the economic performance.


**Table 3.** Performance of the different plants evaluated in this study.

The levelized cost of electricity and CO2 avoidance cost for the different cases are shown in Figure 12. Clearly, the ICR outlet temperature has a large impact on the economic performance of the plant. Interestingly, the relative increase in capital costs from the ICR-1300 case to the ICR-1000 case (50%) is greater than the relative increase in fuel costs (25%). This is due to the large increase in air flowrate required by the cases with lower ICR outlet temperatures. The larger air flowrate increases the costs of ICR reactors, gas turbines, and heat recovery steam generators, adding to the increase in specific capital costs caused by a reduction in net electric efficiency.

**Figure 12.** Levelized cost of electricity (LCOE) and CO2 avoidance costs (CAC).

CO2 avoidance costs range from \$58/ton to \$141/ton. For comparison, a thorough review by Rubin et al. [63] found the representative CO2 avoidance cost for NGCC plants with post-combustion

capture to be \$107/ton (after adjusting to 2019 currency and adding CO2 transport and storage costs). The ICR-1300 case, if it is technically possible, therefore appears attractive relative to this benchmark, whereas the \$91/ton CO2 avoidance cost of the more realistic ICR-1150 case offers marginal benefits.

As outlined earlier, all the reactor sizes investigated performed well in terms of fuel conversion, and the primary difference was the amount of particle elutriation. Preventing particle elutriation will require the addition of a cyclone (possibly an internal cyclone placed in the ICR freeboard), including smaller reactor sizes needing cyclones with a higher separation efficiency to keep escaped fines constant. Inclusion of such a cyclone was not studied here, but the effect of the change in reactor size was investigated.

Figure 13 shows that higher reactor costs have a relatively small impact on overall economic performance of the plant. An increase in reactor costs from 108 to 233 M\$ only increased the levelized cost of electricity by \$5.5/MWh and the CO2 avoidance cost by \$16.5/ton. The uncertainty in the reactor cost estimation was also highlighted in the methodology description, particularly the factor 3 that was used to account for the cost increase of including the ICR internal structure in a Ni-alloy process vessel. For perspective, the range of reactor costs shown in Figure 13 is equivalent to a wide range in this tuning factor of 1.7–6.4, indicating that the overall plant performance is not overly sensitive to uncertainties in the ICR cost assessment.

**Figure 13.** The effect of reactor scale factor (reactor diameter and height relative to the base case) on economic performance indicators.

#### **4. Discussion**

The energy and process industry is in a state of flux due to climate change concerns, increasingly stringent environmental regulations, and constraints on raw material availability. However, a swift response to this changing market environment is challenging due to the multidecade scale-up and demonstration timeframes typical in the industry. This study applied a simulation-based approach that could alleviate this challenge to the design and economic evaluation of a novel reactor concept for clean power production from natural gas.

Specifically, filtered two-fluid model (fTFM) simulations were used to design and further investigate a relatively complex internally circulating reactor (ICR) at an industrial scale. Naturally, there are many uncertainties with modeling such a complex multiphase reactor at industrial scales and operating conditions, but the simulation results served to increase confidence in the ICR technology and offered several insights on the practical operation of a large-scale ICR. Such insights would be very expensive and time-consuming to gain via gradual experimental scale-up and demonstration.

The ICR has many different design and operating parameters that can be adjusted to optimize performance. Simulations can allow for cost-effective completion of the very large number of design iterations that will be required to successfully optimize within this large parameter space. Model accuracy is the main uncertainty in this strategy, making validation against large-scale experimental data a high priority. Unfortunately, such data is scarce, often not publicly available, and generally not detailed enough for proper model validation. Given the large fundamental benefits of such a simulation-based design approach, efforts to collect such data from large-scale reactors for further fTFM development and validation is highly recommended. Furthermore, there is substantial scope for further development of fTFM closures for improved accuracy and generality, which would reduce the uncertainty of simulation-based investigations of industrial-scale fluidized beds. Accounting for polydispersity in fTFMs is an especially rich area for research, considering the importance of polydispersity to many industrial applications, the complexity of its subgrid modeling, and the limited amount of research on the topic to date.

Sizing of the ICR via simulations allowed for an economic assessment of the ICR integrated into a natural gas combined cycle (NGCC) power plant for clean power production using the chemical looping combustion (CLC) principle. The primary challenge of using CLC in combined cycles is the maximum temperature achievable in the reactor, which limits the turbine inlet temperature and, hence, the power cycle efficiency. Results from the economic assessment illustrated this with a decrease in CO2 avoidance costs from \$141/ton to \$58/ton if the ICR operating temperature could be increased from 1000 to 1300 ◦C. These results emphasize that one of the main focusses of oxygen carrier development effort should be the high temperature durability of the materials.

Since it is doubtful that ICR operating temperatures can increase far beyond 1200 ◦C and gas turbine technology keeps progressing beyond 1600 ◦C to allow for higher efficiencies, the maximum reactor temperature is a serious limitation. One potential solution is to combust additional fuel after the ICR to further increase the gas temperature, as outlined in Khan et al. [34]. The use of such a combustor can eliminate most of the energy penalty, although this is at the expense of increased CO2 emissions if natural gas is used for the extra firing, or at the expense of added costs if hydrogen is used. The economic implications of this trade-off will be explored in future work involving detailed modeling of a modern gas turbine with a combustor outlet temperature exceeding 1600 ◦C.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2227-9717/7/10/723/s1, Document: Hydrodynamic closure development; Video S1: Oxygen carrier circulation; Video S2: Gas leakage; Video S3: Effect of reduction section outlet overpressure; Video S4: Effect of reactor size.

**Author Contributions:** Conceptualization, J.H.C. and S.C.; data curation, J.H.C. and M.N.K.; formal analysis, J.H.C., M.N.K., and S.C.; funding acquisition, S.A.; investigation, J.H.C. and M.N.K.; methodology, J.H.C., M.N.K., and S.C.; project administration, S.A.; software, J.H.C.; supervision, S.C. and S.A.; validation, J.H.C. and M.N.K.; visualization, J.H.C. and M.N.K.; writing—original draft, J.H.C., M.N.K., and S.C.; writing—review and editing, J.H.C., S.C., and S.A.

**Funding:** This research was funded by the Research Council of Norway under the CLIMIT program, grant number: 255462.

**Acknowledgments:** The resolved simulations used for closure development and verification in this study were performed on computing resources provided by UNINETT Sigma2—the National Infrastructure for High Performance Computing and Data Storage in Norway.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **Appendix A. Reactive Closure Development**

#### *Appendix A.1. Resolved Two-Fluid Model (TFM) Simulations*

The subgrid reactive closure was developed from resolved TFM simulations using a similar setup, material properties, and methodology to what has been done in previous work [41]; therefore, the detailed description of the simulation setup and statistical analysis of the data is not repeated here. However, new resolved simulations were performed with some changes in place compared to the aforementioned study in order to collect better data for reaction modeling. Firstly, the height of the periodic simulation domain was doubled to 1.28 m to increase the region for collecting data. Secondly, the species boundary conditions for the reactants were set to a fixed value at the bottom boundary. This ensured that the reactants were never depleted in the periodic simulations, allowing data to be collected continuously. However, imposing this fixed boundary condition influenced the effectiveness factors predicted near the boundary. Therefore, data was only collected from heights between 0.4 and 0.88 m, where the effectiveness factor was found to be independent of the domain height. Four simulations were performed with domain-averaged solids volume fractions of 0.05, 0.1, 0.2, and 0.4 to provide continuous data over the entire range of filtered solids volume fractions.

In the present study, isothermal conditions were assumed, and five solids-catalyzed reactions were considered. In each of the reactions, a single product gas phase species was converted to another product species. All reactions were independent from each other and each species only took part in one reaction. Three first-order reactions were considered with different rate constants, as well as a 0.5th order and a 2nd order reaction. The different reactions are summarized in Table A1. Consequently, the species transport equation can be written as follows for the reactants and products, respectively.

$$\frac{\partial}{\partial t} \Big( a\_{\mathcal{S}} \rho\_{\mathcal{S}} X\_i \Big) + \nabla \cdot \Big( a\_{\mathcal{S}} \rho\_{\mathcal{S}} \overrightarrow{\boldsymbol{\upsilon}}\_{\mathcal{S}} X\_i \Big) = \nabla \cdot \Big( a\_{\mathcal{S}} \rho\_{\mathcal{S}} D\_i \nabla X\_i \Big) - k\_i a\_{\mathcal{S}} \mathbb{C}\_i^{\eta\_i} M\_{i\boldsymbol{\prime}} \tag{A1}$$

$$
\frac{\partial}{\partial t} \left( a\_{\mathcal{S}} \rho\_{\mathcal{S}} X\_i \right) + \nabla \cdot \left( a\_{\mathcal{S}} \rho\_{\mathcal{S}} \overrightarrow{\boldsymbol{\upsilon}}\_{\mathcal{S}} X\_i \right) = \nabla \cdot \left( a\_{\mathcal{S}} \rho\_{\mathcal{S}} D\_i \nabla X\_i \right) + k\_i a\_{\mathcal{S}} \mathbb{C}\_i^{w\_i} M\_i. \tag{A2}
$$


**Table A1.** Summary of the reactions that were considered.

#### *Appendix A.2. Closure Development*

It was previously established that a simple effectiveness factor closure is sufficient for use in reactive fTFMs [41]. Therefore, a simple one-marker closure, similar to the one described previously [41], was developed for the reference reaction, 1 mid. It was found that the following expression can accurately predict the effectiveness factor for the reference reaction as a function of the filter size and the filtered solids volume fraction:

$$-\log(\eta) = \left(\frac{2}{\pi}\right)^3 \text{atan}\left(\mathbf{x}\_1 \Delta\_f^{\*,\mathbf{x}\_2} \overline{a}\_s\right) \text{atan}\left(\mathbf{x}\_3 \Delta\_f^{\*,\mathbf{x}\_2} \text{max}\left(\mathbf{a}\_{\mathbf{s},\max} - \overline{a}\_s, 0\right)\right) \text{atan}\left(\mathbf{x}\_4 \Delta\_f^{\*}\right) \mathbf{x}\_5. \tag{A3}$$

An excellent fit (*R*<sup>2</sup> = 0.991) was obtained against the binned (conditionally-averaged) data from the resolved TFM simulations, using the following model coefficient values: *x*<sup>1</sup> = 2.162, *x*<sup>2</sup> = 0.3204, *x*<sup>3</sup> = 5.504, *x*<sup>4</sup> = 1.163, *x*<sup>5</sup> = 1.281, and αs,max = 0.5621.

Next, it was found that the length factor, L, in Equation (12) can be closed as a function of the filter size. This is understandable, as the average size of the subgrid clusters will increase as the grid size increases in the coarse-grid simulations. Using the methodology described in Section 2.1.2, the best fit to the data from all five reactions was obtained with the following length factor closure:

$$
\mathcal{L} = \left(\frac{2}{\pi}\right) \text{atan}\left(\mathbf{x}\_1 \Delta\_f^\*\right) \mathbf{x}\_2. \tag{A4}
$$

A reasonably good fit of *R*<sup>2</sup> = 0.864 was obtained over all the reaction data with the coefficients *x*<sup>1</sup> = 0.08768 and *x*<sup>2</sup> = 71.24, and when evaluating the minimum effectiveness factor of the reference closure (used in Equations (15) and (16)) at a filtered solids volume fraction of 0.3403.

The good fit against the binned data for such a large variety of reactions implies that applying an analogy to intraparticle mass transfer works sufficiently well for generalizing the reactive fTFM closure. This is further shown in Figure A1, which compares the binned reactions rates from the resolved simulations to the model predictions as a function of the reactant concentration. Both axes are shown on a base-10 log scale to show the entire range of data collected. It is clearly seen from Figure A1a that the importance of the effectiveness factor due to sub-grid effects increases as the reaction rate increases. Furthermore, the generalized effectiveness factor closure performs well in capturing this effect.

Next, Figure A1b shows the importance of accounting for the reaction rate order in the effectiveness factor closure, especially for the 0.5 order reaction. Although some deviations exist in the 0.5 order case, the intraparticle mass transfer analogy generally captures the effect of changing reaction order well. Furthermore, the generalized closure clearly presents a substantial improvement over simply using the closure for a first order reaction, noting that all previous fTFM closures for the reaction effectiveness factor only considered first order reactions.

**Figure A1.** The reaction rate plotted against the concentration, comparing the binned data from resolved simulations (symbols) to model predictions (solid lines) for (**a**) first order and (**b**) other order reactions. The dashed lines show the predictions when using the effectiveness factor from the reference reaction (1 mid) and the dotted line the predictions without an effectiveness factor. Data is shown for the largest filter size considered (Δˆ *<sup>f</sup>* = 10.4) and for an intermediate filtered solids volume fraction (α*<sup>s</sup>* = 0.30 ).

#### *Appendix A.3. Closure Verification*

The generalized reactive effectiveness factor closure was verified by comparing coarse-grid simulation predictions to that of resolved simulations for the fast bubbling case discussed in the Supplementary Material. Figure A2a shows the importance of subgrid modeling for the reactions, showcasing large overpredictions of the scaled conversion occurring when a closure is neglected, especially for large grid sizes and for the case with large subgrid corrections (fast reaction rate and/or low reaction order). Figure A2b shows that even though both the reference and generalized effectiveness factor closures offer a significant improvement over the case with no closure, the generalized closure consistently outperforms the reference closure for all reactions (except for the reference reaction, C, where the closures are identical and the performance is the same). It is interesting to note that an excellent prediction of the conversion could be obtained despite larger inaccuracies in the hydrodynamics

prediction (Figure S4b in the Supplementary Material), suggesting that the effectiveness factor closure is robust.

**Figure A2.** Evaluation of different approaches to modeling the filtered reaction rate in coarse-grid simulations of fast bubbling flow: neglecting subgrid effects on the filtered reaction rate (None), using an effectiveness factor for a reference first-order reaction (Reference), and using the generalized effectiveness factor closure described in this paper (General). Three different grid sizes were considered in the coarse-grid simulations. (**a**) Comparison of the scaled conversion for the different cases; (**b**) comparison of the percentage error of the generalized and reference closure predictions when compared to the benchmark resolved simulations.

The generalized closure is reasonably accurate in predicting the scaled conversion for all reactions, although some discrepancies occur, especially at large grid sizes. However, it is likely that these discrepancies are mainly due to the grid size being too coarse to resolve important macro-scale

hydrodynamic flow structures [27], and may therefore not be primarily due to shortcomings of the reactive closure. Future studies will focus on a more detailed verification of reactive fTFM closures.

The results of this section emphasize the benefit of the generalized effectiveness factor closure; therefore, such a closure is highly recommended for future studies of industrial-scale fluidized bed reactors using reactive fTFM simulations.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
