Application of Pseudohomogeneous and Heterogeneous Models in Assessing the Behavior of a Fluidized-Bed Catalytic Reactor

Abstract

Comparative analysis of the steady-state and transient properties of a bubbling fluidized-bed catalytic reactor obtained according to different mathematical models of the emulsion zone was performed to verify the commonly used assumption regarding the pseudohomogeneous nature of this zone. Four different mathematical models of the fluidized-bed reactor dynamics were formulated, based on different thermal and diffusional conditions at the gas-solid interface and within the catalyst pellet, namely the model based on the assumption of pseudohomogeneous character for the emulsion zone, and a group of two-scale models accounting for the heterogeneous character of this zone. It was demonstrated that, while the pseudohomogeneous model of the emulsion zone predicts almost identical behavior of the reactor at steady-state as the proposed heterogeneous models, it may fail in the prediction of the reactor start-up behavior, especially when dealing with highly exothermic processes run at relatively high fluidization velocity.

1. Introduction

Bubbling fluidized-bed reactors (BFBR) for solid-catalyzed gas phase processes are unquestionably, units of great importance in the chemical and petrochemical industries, as well as in energy conversion. Due to its hydrodynamic properties and tendency to equalize temperature, the fluidized bed provides an excellent environment for running highly exothermic solid-catalyzed chemical reactions since any hot spots are rapidly quenched. Another advantage is the fluid-like behavior of the bed that makes the BFBR attractive for running chemical processes that require continuous regeneration of the catalyst particles [1,2]. These features have favored a wide industrial application of the BFBR, among others, in the processes of olefin polymerization [2,3,4], oxidation of hydrocarbons [1,2,5], propylene ammoxidation to acrylonitrile or fluidized-bed cracking [1,2]. Besides solid-catalyzed chemical processes, fluidized-bed units are widely used in energy industry for non-catalytic combustion and gasification of solid fuels [6,7], and in material processing for powder granulation [8,9,10] or drying [11].

Effective design, operation and control of these apparatuses require a comprehensive analysis of their steady-state and transient behavior made with the aid of an adequate mathematical model. Many predictive models of catalytic BFBR [12,13,14,15,16] are based on a classical two-phase concept of fluidization, whose origins date back to the 1950s [17,18]. According to the simplest two-phase model, the BFB consists of two phases (zones): a bubble phase and an emulsion phase (also called a dense phase), which is usually treated as a pseudohomogeneous medium, i.e., both intraparticle and interphase gradients of the concentration and temperature are neglected in the emulsion model. Mass and heat transport resistances are usually disregarded a priori when modeling the catalytic BFBR, and such a simplification is motivated by the small size of the particles used as the bed material. However, it should be borne in mind that too simplified description of intraparticle and interphase transport phenomena may lead to erroneous prediction of the performance of these solid-gas reactors. Some early studies concerning mathematical modeling of catalytic fluidized-bed reactors [19,20] take indeed into account the external resistance to mass and heat transport at the pellet surface, however most of the later studies concerning catalytic processes completely neglect both intraparticle mass and heat diffusion as well as external transport resistances at the particle surface [12,13,14,15,16]. The assumption of a pseudohomogeneous nature of the emulsion phase is usually not supported by any analysis of the transport resistances. One of the few exceptions are the works of Fernandes and Lona [3,4] concerning mathematical modeling of the ethylene polymerization in a BFBR. The authors compared the results obtained using a full heterogeneous model of the emulsion zone with the pseudohomogeneous model of this zone. They demonstrated that, for low values of the mean residence time of the polymer particles in the reactor, the emulsion phase should be described using a heterogeneous model. When the mean residence time of particles in the reactor is relatively long, then the results obtained using both models are comparable.

Considering the state of art of mathematical modeling of catalytic FBR operating in the bubbling regime, there is still a necessity of revision, especially in relation to the transient analysis, of the commonly adopted pseudohomogeneous model of the emulsion zone. Heterogeneous models of the emulsion zone incorporated in the reactor models based on two-phase theory are nowadays computationally affordable, even when dealing with real-time applications, such as control or transient behavior prediction. Aim of this study is to show that it is necessary to move away from the—usually “blind”—application of pseudohomogeneous models. While the effect of transport phenomena around and inside the catalyst pellets on the reactor behavior was analyzed quite intensively in the past in relation to fixed-bed reactors [21], there is a lack of similar studies referred to fluidized-bed reactors. Moreover, although the influence of transport resistances is much smaller, it still needs particular attention when dealing with highly exothermic processes. Other systems for which the necessity to employ fully heterogeneous model may arise are adsorptive reactors and catalytic reactors [22] integrating two or more types of active sites. When different functionalities are integrated within a structured hybrid pellet [23,24], then the intrapellet distribution of the functionalities needs to be taken into account to predict correctly the reactor behavior. Therefore, in this study, the steady-state and transient properties of a catalytic BFBR determined using mathematical models accounting for various diffusional and thermal conditions at the gas-solid interface and within the catalyst pellet, were compared, with the aim of evaluating the applicability of the pseudohomogeneous model for the prediction of the reactor behavior.

2. Materials and Methods

2.1. Mathematical Models of a Catalytic BFBR

Figure 1a shows a schematic diagram of a catalytic BFBR together with the basic notation used in the mathematical models formulated below. The models consider various diffusional and thermal conditions at the solid-gas interface and inside the catalyst pellet. They are based on a well-established theory of mathematical modelling of solid-catalyzed chemical reactions accompanied by diffusion in porous particles [21]. In all cases, a two-phase bubbling-bed model [1,25] which considers the bed consists of two phases or zones, i.e., a bubble zone and an emulsion zone, was adopted to describe the fluidized bed. To avoid ambiguity hereafter, the term zone is used to refer to the phase (i.e., bubbles or emulsion) present in the fluidized bed, whereas the term phase refers to gas or solid phase only.

The equations used to determine the hydrodynamic properties of the fluidized bed are summarized in

Table 1. Hydrodynamic properties of the BFBR.

Parameter	Equation
Fluidized-bed void fraction at minimum fluidization conditions, ${\epsilon }_{mf}$ [27]	${\epsilon }_{mf}=0.586\cdot \mathrm{A}{\mathrm{r}}^{-0.029}{\left(\frac{{\rho }_g}{{\rho }_p}\right)}^{0.021}$.	(1)
Superficial gas velocity at minimum fluidization conditions, $u_{mf}$ [1]	$\frac{1.75}{{\epsilon }_{mf}^3{\varphi }_p}{\mathrm{Re}}_{p,mf}^2+\frac{150(1-{\epsilon }_{mf})}{{\epsilon }_{mf}^3{\varphi }_p^2}{\mathrm{Re}}_{p,mf}=\mathrm{A}\mathrm{r}$ where $\mathrm{A}\mathrm{r}=\frac{d_p^3{\rho }_g({\rho }_p-{\rho }_g)g}{{\mu }_g}$ and ${\mathrm{Re}}_{p,mf}=\frac{d_p{u}_{mf}{\rho }_g}{{\mu }_g}$.	(2)
Terminal velocity of a falling particle, $u_t$ [28]	${\mathrm{Re}}_{p,t}=\frac{\mathrm{Ar}}{18+0.6\mathrm{A}{\mathrm{r}}^{0.5}}$ where ${\mathrm{Re}}_{p,t}=\frac{d_p{u}_t{\rho }_g}{{\mu }_g}$.	(3)
Effective bubble diameter, $d_b$ [29]	$d_b=d_{b0}+0.5\cdot \xi \cdot H$ where $H=\frac{H_{mf}}{1-\delta }$, $d_{b0}=0.8205{\left(\frac{u_0-u_{mf}}{n_0}\right)}^{0.4}$ and $\xi =0.14{\rho }_p{d}_p\frac{u_0}{u_{mf}}$.	(4)
Bubble fraction in a fluidized bed, $\delta$ [1]	$\delta =\frac{u_0-u_{mf}}{u_b-u_{mf}}$.	(5)
Bubble rise velocity, $u_b$ [1]	$u_b=u_0-u_{mf}+u_{b0}$ where $u_{b0}=2.227\sqrt{d_b}$	(6)
Superficial velocity of gas in the emulsion, $u_e$ [1]	$u_e=\frac{u_0-\delta u_b}{(1-\delta )\cdot {\epsilon }_{mf}}$.	(7)

. More details on the methodology of determination of the hydrodynamic characteristics of the catalyst bed can be found in reference [26], which, however, is limited to a steady-state analysis made with the use of a pseudohomogeneous model only.

In the simplest mathematical model of the catalytic BFBR evaluated in this work, the emulsion zone is treated as a pseudohomogeneous medium (Figure 1b). The adoption of such an assumption results in the omission of equations describing the interphase mass and heat exchange, and equations describing the distribution of concentration and temperature in the catalyst pellet. The second category of the models formulated within this work are the models of the BFBR which account for the heterogeneous nature of the emulsion zone (Figure 1c).

Due to the general character of the analysis, focused mainly on the comparison of reliability of different mathematical models of the emulsion zone and on the evaluation of the effect of the selected parameters onto the results discrepancy, a single first-order irreversible exothermic chemical reaction $A\stackrel{k}{\to }P$ was considered as a test case. The chemical reaction rate equation per unit volume of the catalyst particle with respect to reactant A is expressed as:

$$r_{\mathrm{A}}(C_{\mathrm{A}},T)=k_0\exp \left(-\frac{E}{RT}\right)\cdot C_{\mathrm{A}}$$

(8)

Other assumptions of mathematical models formulated below for the catalytic fluidized-bed reactor operating in the bubbling regime are the following:

• ideal mixing of gas in the emulsion zone;
• plug flow of gas in the bubble zone and quasi-steady state character of this zone, which is motivated by lower mass and thermal inertia of the bubble zone compared to the emulsion zone;
• gas bubbles are spherical, and they have a constant diameter along the bed height;
• values of mass and heat exchange coefficient between the emulsion zone and the bubble zone are calculated according to the assumptions of a three-phase model (

  Table 2. Mass and heat interchange coefficients.

  Parameter	Equation
  Overall coefficient of gas interchange between bubble and emulsion zone, ${\beta }_{g\mathrm{A}}^{be}$ [30]	$\frac{1}{{\beta }_{g\mathrm{A}}^{be}}=\frac{1}{{\beta }_{g\mathrm{A}}^{bc}}+\frac{1}{{\beta }_{g\mathrm{A}}^{ce}}$ where ${\beta }_{g\mathrm{A}}^{ce}=6.78\sqrt{\frac{{\epsilon }_{mf}{D}_{e\mathrm{A}}{u}_b}{d_b^3}}$ and ${\beta }_{g\mathrm{A}}^{bc}=4.5\frac{u_{mf}}{d_b}+10.353\frac{D_{b\mathrm{A}}^{0.5}}{d_b^{1.25}}$.	(9)
  Overall coefficient of heat interchange between bubble and emulsion zone, ${\alpha }_q^{be}$ [30]	$\frac{1}{{\alpha }_q^{be}}=\frac{1}{{\alpha }_q^{bc}}+\frac{1}{{\alpha }_q^{ce}}$ where ${\alpha }_q^{ce}=6.78\sqrt{\frac{{\epsilon }_{mf}{\lambda }_g{\rho }_g{c}_g{u}_b}{d_b^3}}$ and ${\beta }_{g\mathrm{A}}^{bc}=4.5\frac{u_{mf}}{d_b}+10.353\frac{D_{b\mathrm{A}}^{0.5}}{d_b^{1.25}}$.	(10)

  );
• the chemical reaction takes place only in the emulsion zone;
• the chemical reaction may be accompanied by physical sorption of the reactant on the inert support of the catalyst pellet; however, the heat of adsorption is neglected because it is much lower than the heat of chemical reaction;
• the sorption process is at equilibrium, described by a linear isotherm with temperature-dependent adsorption constant;
• contribution of the intraparticle gas phase and the adsorbed phase heat capacities to the rate of increase of the catalyst pellet thermal energy is neglected.

In a single apparatus with no solid recycle, physical adsorption does not influence the loci of steady-state branches; however, it may affect their stability and have a significant impact on the dynamic features of the reactor. The importance of the reactant adsorption on the support was recognized in theoretical studies concerning the transient behavior of fixed-bed reactors [31] and it was also confirmed by experimental observations regarding fluidized-bed catalytic cracking [32], therefore, the occurrence of this phenomenon was accounted for in the models formulated in this work.

2.2. Model of the Bubble Zone

Regardless of the form of the emulsion zone model, under the assumption of a plug flow of gas in the bubble zone and considering a quasi-steady state character of bubbles, the balance equations for this zone can be written as follows:

$$S\delta u_b\frac{d{C}_{\mathrm{A}}^b}{dh}=-S\delta {\beta }_{g\mathrm{A}}^{be}(C_{\mathrm{A}}^b-C_{\mathrm{A}}^e)$$

(11)

$$S\delta u_b{\rho }_g{c}_g\frac{d{T}^b}{dh}=S\delta ({\alpha }_q^{be}+{\beta }_p^{be}{\rho }_p{c}_s)\cdot (T^e-T^b)-S\delta a_q{k}_q(T^b-T_q)$$

(12)

where $h\in [0,H]$, with H being the total height of the fluidized bed and ${\rho }_p=(1-{\epsilon }_p){\rho }_s$ is the effective particle density, with ε_p being the catalyst porosity and ρ_s the solid phase density.

Ordinary differential equations (ODEs) (Equations (11) and (12)) are supplemented with the initial conditions (ICs) resulting from the concentration of reactant A and temperature of the feed stream:

$$C_{\mathrm{A}}^b(h=0)=C_{\mathrm{A}f},T^b(h=0)=T_f.$$

(13)

After the introduction of the following dimensionless variables:

$$X_{\mathrm{A}}=\frac{C_{\mathrm{A}}}{C_{\mathrm{A},ref}}\mathrm{where}{C}_{\mathrm{A},ref}=\frac{P\cdot y_{\mathrm{A}f}}{R\cdot T_{ref}},z=\frac{h}{H}\in [0,1],$$

(14)

Equations (11) and (12) with ICs given by Equation (13) take, respectively, the form:

$$\frac{d{X}_{\mathrm{A}}^b}{dz}=-B_3(X_{\mathrm{A}}^b-X_{\mathrm{A}}^e),$$

(15)

$$\frac{d{T}^b}{dz}=-B_4(T^b-T^e)-Q_2(T^b-T_q),$$

(16)

and:

$$X_{\mathrm{A}}^b(z=0)=X_{\mathrm{A}f},T^b(z=0)=T_f,$$

(17)

where:

$$B_3=\frac{H{\beta }_{g\mathrm{A}}^{be}}{u_b},B_4=\frac{H}{u_b}\left(\frac{{\alpha }_q^{be}}{{\rho }_g{c}_g}+{\beta }_p^{be}{a}_1\right),Q_2=\frac{H{a}_q{k}_q}{u_b{\rho }_g{c}_g},a_1=\frac{{\rho }_p{c}_s}{{\rho }_g{c}_g}.$$

(18)

Separation of variables in Equations (15) and (16) followed by closed-form integration with conditions defined by Equation (17) gives the following expressions describing the distributions of dimensionless concentration of reactant A, $X_{\mathrm{A}}^b(z)$, and temperature, $T^b(z)$:

$$X_{\mathrm{A}}^b(z)=X_{\mathrm{A}}^e\left[1-\exp (-B_{3}z)\right]+X_{\mathrm{A}f}\cdot \exp (-B_{3}z),$$

(19)

$$T^b(z)=\frac{B_4{T}^e+Q_2{T}_q}{B_4+Q_2}+\left[T_f-\frac{B_4{T}^e+Q_2{T}_q}{B_4+Q_2}\right]\cdot \exp \left[-(B_4+Q_2)z\right].$$

(20)

2.3. Pseudohomogeneous Model of the Emulsion Zone (Model P)

Beyond the general assumptions formulated at the beginning of Section 2.1, basis for the construction of the simplest model of the emulsion zone, i.e., the pseudohomogeneous model (Figure 1b), is the assumption that there is no transport resistance between interstitial emulsion gas and catalyst particles, and within catalyst particles (

Table 3. Main assumptions of the mathematical models.

Model of the Emulsion Zone	Internal Transport Resistance		External Transport Resistance		Model Equations
	Mass	Heat	Mass	Heat
P	-	-	-	-	(25), (26), (34)
H1	✓	✓	✓	✓	(42)–(45), (48), (49), (53)–(55)
H2	✓	-	✓	✓	(42), (44), (48), (49), (53), (55), (58), (59)
H3	-	-	✓	✓	(48), (49), (55), (63), (64), (66)

).

Under the assumption that the emulsion zone can be considered as a pseudohomogeneous medium, i.e., when the intensity of transport processes is so high that the chemical process is controlled by kinetics, the dynamic mass and energy balances are:

$$\begin{gathered} SH(1-\delta ){\epsilon }_{mf}\frac{d{C}_{\mathrm{A}}^e}{dt}+SH(1-\delta )(1-{\epsilon }_{mf})\left({\epsilon }_p\frac{d{C}_{\mathrm{A}}^e}{dt}+{\rho }_p\frac{d{q}_{\mathrm{A}}}{dt}\right)=S(1-\delta ){\epsilon }_{mf}{u}_e(C_{\mathrm{A}f}-C_{\mathrm{A}}^e)+ \\ +S\delta {\beta }_{g\mathrm{A}}^{be}{\displaystyle \underset{0}{\overset{H}{\int }}(C_{\mathrm{A}}^b-C_{\mathrm{A}}^e)}dh-SH(1-\delta )(1-{\epsilon }_{mf})\cdot r_{\mathrm{A}}(C_{\mathrm{A}}^e,T^e), \end{gathered}$$

(21)

$$\begin{gathered} SH(1-\delta )\left[(1-{\epsilon }_{mf}){\rho }_p{c}_s+{\epsilon }_{mf}{\rho }_g{c}_g\right]\frac{d{T}^e}{dt}=S(1-\delta ){\epsilon }_{mf}{u}_e{\rho }_g{c}_g(T_f-T^e)- \\ -S\delta ({\alpha }_q^{be}+{\beta }_p^{be}{\rho }_p{c}_s){\displaystyle \underset{0}{\overset{H}{\int }}(T^e-T^b)}dh+SH(1-\delta )(1-{\epsilon }_{mf})\cdot r_{\mathrm{A}}(C_{\mathrm{A}}^e,T^e)\cdot (-\Delta h)- \\ -SH(1-\delta )a_q{k}_q(T^e-T_q). \end{gathered}$$

(22)

The derivative $d{q}_{\mathrm{A}}/dt$ in Equation (21) can we written as:

$$\frac{d{q}_{\mathrm{A}}}{dt}=\frac{d{q}_{\mathrm{A}}}{d{C}_{\mathrm{A}}^e}\frac{d{C}_{\mathrm{A}}^e}{dt},$$

(23)

where $d{q}_{\mathrm{A}}/d{C}_{\mathrm{A}}^e$ is the slope of the sorption equilibrium curve. For a linear isotherm with a temperature-dependent equilibrium constant, $K_a$, the slope expression can be rewritten as:

$$\frac{d{q}_{\mathrm{A}}}{d{C}_{\mathrm{A}}^e}=K_a(T^e)=K_{a0}\exp \left(\frac{E_a}{R{T}^e}\right).$$

(24)

Introduction of the above defined dimensionless variables (Equation (14)) and expressions describing the distribution of reactant A concentration and temperature in the bubbles (Equations (19) and (20)) into Equations (21) and (22), followed by closed-form integration of integrals leads to the following final form of the pseudohomogeneous model of the emulsion zone (hereinafter referred to as model P):

$$\frac{d{X}_{\mathrm{A}}^e}{dt}=a_2(X_{\mathrm{A}f}-X_{\mathrm{A}}^e)+{\phi }_1(X_{\mathrm{A}}^e)-a_3{R}_{\mathrm{A}}(X_{\mathrm{A}}^e,T^e),$$

(25)

$$\frac{d{T}^e}{dt}=a_5(T_f-T^e)-{\phi }_2(T^e)+a_4{R}_{\mathrm{A}}(X_{\mathrm{A}}^e,T^e)\cdot (-\Delta h)-Q_1(T^e-T_q).$$

(26)

where ${\phi }_1(X_{\mathrm{A}}^e)$, ${\phi }_2(T^e)$ and $R_{\mathrm{A}}(X_{\mathrm{A}}^e,T^e)$ are:

$${\phi }_1(X_{\mathrm{A}}^e)=\frac{B_1}{B_3}\left[\exp (-B_3)-1\right]\cdot (X_{\mathrm{A}}^e-X_{\mathrm{A}f}),$$

(27)

$${\phi }_2(T^e)=\frac{B_2}{B_4+Q_2}\left[Q_2(T^e-T_q)+\left(T_f-\frac{B_4{T}^e+Q_2{T}_q}{B_4+Q_2}\right)\cdot \left[\exp (-(B_4+Q_2))-1\right]\right],$$

(28)

$$R_{\mathrm{A}}(X_{\mathrm{A}}^e,T^e)=k_0\exp \left(-\frac{E}{R{T}^e}\right)\cdot X_A^e,$$

(29)

whereas the parameters appearing in Equations (25)–(28) are defined as:

$$a_1=\frac{{\rho }_p{c}_s}{{\rho }_g{c}_g},a_2=\frac{{\epsilon }_{mf}{u}_e}{H[{\epsilon }_{mf}+(1-{\epsilon }_{mf})({\epsilon }_p+{\rho }_p{K}_a)]},a_3=\frac{1-{\epsilon }_{mf}}{{\epsilon }_{mf}+(1-{\epsilon }_{mf})({\epsilon }_p+{\rho }_p{K}_a)},$$

(30)

$$a_4=\frac{(1-{\epsilon }_{mf})C_{\mathrm{A},ref}}{(1-{\epsilon }_{mf}+{\epsilon }_{mf}/a_1){\rho }_p{c}_s},a_5=\frac{{\epsilon }_{mf}{u}_e}{H[(1-{\epsilon }_{mf})a_1+{\epsilon }_{mf}]},Q_1=\frac{a_q{k}_q}{(1-{\mathsf{\epsilon }}_{mf}+{\mathsf{\epsilon }}_{mf}/a_1){\rho }_p{c}_s},$$

(31)

$$B_1=\frac{\delta {\beta }_{g\mathrm{A}}^{be}}{(1-\delta )[{\epsilon }_{mf}+(1-{\epsilon }_{mf})({\epsilon }_p+{\rho }_p{K}_a)]},B_2=\frac{\delta }{(1-\delta )(1-{\epsilon }_{mf}+{\epsilon }_{mf}/a_1)}\left(\frac{{\alpha }_q^{be}}{{\rho }_p{c}_s}+{\beta }_p^{be}\right),$$

(32)

$$B_3=\frac{H{\beta }_{g\mathrm{A}}^{be}}{u_b},B_4=\frac{H}{u_b}\left(\frac{{\alpha }_q^{be}}{{\rho }_g{c}_g}+{\beta }_p^{be}{a}_1\right),Q_2=\frac{H{a}_q{k}_q}{u_b{\rho }_g{c}_g}.$$

(33)

The system of ODEs (Equations (25) and (26)) must be supplemented by ICs that determine the values of initial concentration and temperature in the emulsion zone, namely:

$$X_{\mathrm{A}}^e(t=0)=X_{\mathrm{A}0}^e,T^e(t=0)=T_0^e.$$

(34)

2.4. Heterogeneous Models of the Emulsion Zone

The second category of models are the models of a catalytic BFBR with a heterogeneous model of the emulsion zone, according to which gas and catalyst particles are considered as two separate phases (Figure 1c). Assuming that the chemical reaction takes place in the emulsion zone only, the balance equations for the bubbles are identical to balance Equations (11) and (12) formulated for the reactor model with pseudohomogeneous model of the emulsion zone.

The form of the mass and heat balance equations for the emulsion depends on the diffusional and thermal conditions prevailing in the catalyst particle and its surroundings (Table 3). To limit the number of analyzed models, all formulated below heterogeneous models of the emulsion zone account both for external mass and heat transfer resistance.

2.4.1. Heterogeneous Model with a Distributed-Parameter Model of the Pellet (Model H1)

Considering both the resistance to mass and heat exchange between catalyst pellet and gas flowing through the emulsion, as well as intraparticle resistance to mass and heat transport (Figure 2) and accounting for the assumptions formulated in Section 2.1, the behavior of the emulsion zone is described by the following equations of a single pellet:

$${\epsilon }_p\frac{\partial C_{\mathrm{A}}^p}{\partial t}+{\rho }_p\frac{\partial q_{\mathrm{A}}}{\partial t}=D_{ef}\left(\frac{{\partial }^2{C}_{\mathrm{A}}^p}{\partial r^2}+\frac{2}{r}\frac{\partial C_{\mathrm{A}}^p}{\partial r}\right)-r_{\mathrm{A}}(C_{\mathrm{A}}^p,T^p),$$

(35)

$${\rho }_p{c}_s\frac{\partial T^p}{\partial t}={\lambda }_{ef}\left(\frac{{\partial }^2{T}^p}{\partial r^2}+\frac{2}{r}\frac{\partial T^p}{\partial r}\right)+r_{\mathrm{A}}(C_{\mathrm{A}}^p,T^p)\cdot (-\Delta h),$$

(36)

with boundary conditions (BCs) in the center ($r=0$) and at the pellet surface ($r=L_p$):

$$\frac{\partial C_{\mathrm{A}}^p(0,t)}{\partial r}=0,D_{ef}\frac{\partial C_{\mathrm{A}}^p(L_p,t)}{\partial r}=k_m\left[C_{\mathrm{A}}^e(t)-C_{\mathrm{A}}^p(L_p,t)\right],$$

(37)

$$\frac{\partial T^p(0,t)}{\partial r}=0,{\lambda }_{ef}\frac{\partial T^p(L_p,t)}{\partial r}={\alpha }_q\left[T^e(t)-T^p(L_p,t)\right],$$

(38)

coupled with mass and heat balance equations of the gas phase in the emulsion zone:

$$\begin{gathered} SH(1-\delta ){\epsilon }_{mf}\frac{d{C}_{\mathrm{A}}^e}{dt}=S(1-\delta ){\epsilon }_{mf}{u}_e(C_{\mathrm{A}f}-C_{\mathrm{A}}^e)+S\delta {\beta }_{g\mathrm{A}}^{be}{\displaystyle \underset{0}{\overset{H}{\int }}(C_{\mathrm{A}}^b-C_{\mathrm{A}}^e)}dh- \\ -\frac{6SH(1-\delta )(1-{\epsilon }_{mf})}{d_p}{k}_m\left(C_{\mathrm{A}}^e-C_{\mathrm{A}}^p(L_p)\right), \end{gathered}$$

(39)

$$\begin{gathered} SH(1-\delta ){\epsilon }_{mf}{\rho }_g{c}_g\frac{d{T}^e}{dt}=S(1-\delta ){\epsilon }_{mf}{u}_e{\rho }_g{c}_g(T_f-T^e)-S\delta ({\alpha }_q^{be}+{\beta }_p^{be}{\rho }_p{c}_s){\displaystyle \underset{0}{\overset{H}{\int }}(T^e-T^b)}dh+ \\ +\frac{6SH(1-\delta )(1-{\epsilon }_{mf})}{d_p}{\alpha }_q\left(T^p(L_p)-T^e\right)-SH(1-\delta )a_q{k}_q(T^e-T_q). \end{gathered}$$

(40)

After the introduction of the previously defined dimensionless variables (Equation (14)) and dimensionless concentration of reactant A in the catalyst pellet, ${\beta }_{\mathrm{A}}$, and coordinate, ζ:

$${\beta }_{\mathrm{A}}=\frac{C_{\mathrm{A}}^p}{C_{\mathrm{A},ref}},\zeta =\frac{r}{L_p}\in [0,1],$$

(41)

the partial differential equations (PDEs) (Equations (35) and (36)) describing the concentration and temperature profiles within the catalyst pellet take the following form:

$$\frac{\partial {\beta }_{\mathrm{A}}}{\partial t}={\Gamma }_m\cdot \left[\frac{{\partial }^2{\beta }_{\mathrm{A}}}{\partial {\zeta }^2}+\frac{2}{\zeta }\frac{\partial {\beta }_{\mathrm{A}}}{\partial \zeta }-{\Phi }^2\frac{R_{\mathrm{A}}({\beta }_{\mathrm{A}},T^p)}{R_{\mathrm{A},ref}}\right],$$

(42)

$$\frac{\partial T^p}{\partial t}={\Gamma }_q\cdot \left[\frac{{\partial }^2{T}^p}{\partial {\zeta }^2}+\frac{2}{\zeta }\frac{\partial T^p}{\partial \zeta }+\delta \cdot {\Phi }^2\frac{R_{\mathrm{A}}({\beta }_{\mathrm{A}},T^p)}{R_{\mathrm{A},ref}}\right],$$

(43)

with boundary conditions:

$$\frac{\partial {\beta }_{\mathrm{A}}(0,t)}{\partial \zeta }=0,\frac{\partial {\beta }_{\mathrm{A}}(1,t)}{\partial \zeta }={\mathrm{Bi}}_m\left[X_{\mathrm{A}}^{\mathrm{e}}(t)-{\beta }_{\mathrm{A}}(1,t)\right],$$

(44)

$$\frac{\partial T^p(0,t)}{\partial \zeta }=0,\frac{\partial T^p(1,t)}{\partial \zeta }={\mathrm{Bi}}_q\left[T^e(t)-T^p(1,t)\right],$$

(45)

where:

$${\mathrm{Bi}}_m=\frac{L_p{k}_m}{D_{ef}},{\mathrm{Bi}}_q=\frac{L_p{\alpha }_q}{{\lambda }_{ef}},{\Gamma }_m=\frac{D_{ef}}{L_p^2}\frac{1}{{\epsilon }_p+{\rho }_p{K}_a},{\Gamma }_q=\frac{{\lambda }_{ef}}{L_p^2{\rho }_p{c}_s},$$

(46)

$$\delta =\frac{D_{ef}{C}_{\mathrm{A},ref}(-\Delta h)}{{\lambda }_e},{\Phi }^2=\frac{L_p^2{R}_{\mathrm{A},ref}}{D_{ef}},R_{\mathrm{A}}=k_0\exp \left(-\frac{E}{R{T}^p}\right)\cdot {\beta }_{\mathrm{A}},R_{\mathrm{A},ref}=k_0\exp \left(-\frac{E}{R{T}_{ref}}\right).$$

(47)

Performing similar mathematical operations as in the case of the pseudohomogeneous model of the emulsion zone, the balance equations of interstitial emulsion gas (Equations (39) and (40)) for the heterogeneous model may be transformed into the following form:

$$\frac{d{X}_{\mathrm{A}}^e}{dt}={a^{\prime }}_2(X_{\mathrm{A}f}-X_{\mathrm{A}}^e)+{\phi }_1(X_{\mathrm{A}}^e)-{a^{\prime }}_3\left(X_{\mathrm{A}}^e-{\beta }_{\mathrm{A}}(1)\right),$$

(48)

$$\frac{d{T}^e}{dt}={a^{\prime }}_2(T_f-T^e)-{\phi }_2(T^e)+{a^{\prime }}_4\left(T^p(1)-T^e\right)-{Q^{\prime }}_1(T^e-T_q),$$

(49)

where the functions ${\phi }_1(X_{\mathrm{A}}^e)$ and ${\phi }_2(T^e)$ are given, respectively, by Equations (27) and (28), however, for the heterogeneous model parameters $B_1$ and $B_2$ are replaced, respectively, by ${B^{\prime }}_1$ and ${B^{\prime }}_2$ defined as:

$${B^{\prime }}_1=\frac{\delta {\beta }_{g\mathrm{A}}^{be}}{(1-\delta ){\epsilon }_{mf}},{B^{\prime }}_2=\frac{\delta }{(1-\delta ){\epsilon }_{mf}}\left(\frac{{\alpha }_q^{be}}{{\rho }_g{c}_g}+{\beta }_p^{be}{a}_1\right).$$

(50)

The other parameters appearing in Equations (48) and (49) are defined as follows:

$${a^{\prime }}_2=\frac{u_e}{H},{a^{\prime }}_3=\frac{6(1-{\epsilon }_{mf})}{{\epsilon }_{mf}{d}_p}{k}_m,{a^{\prime }}_4=\frac{6(1-{\mathsf{\epsilon }}_{mf}){\alpha }_q}{{\mathsf{\epsilon }}_{mf}{d}_p{\rho }_p{c}_s},{Q^{\prime }}_1=\frac{a_q{k}_q}{{\epsilon }_{mf}{\rho }_g{c}_g},$$

(51)

where the interphase mass and heat transfer coefficient are determined from the relations:

$$k_m=\frac{\mathrm{Sh}\cdot D_{e\mathrm{A}}}{d_p},{\alpha }_q=\frac{\mathrm{Nu}\cdot {\lambda }_g}{d_p}.$$

(52)

Resolving the equations of the heterogeneous model requires the formulation of appropriate ICs for the catalyst pellet and the interstitial emulsion gas, namely:

$${\beta }_{\mathrm{A}}(\zeta ,t=0)={\beta }_{\mathrm{A}0},\zeta \in [0,1],$$

(53)

$$T^p(\zeta ,t=0)=T_0^p,\zeta \in [0,1],$$

(54)

$$X_{\mathrm{A}}^e(t=0)=X_{\mathrm{A}0}^e,T^e(t=0)=T_0^e$$

(55)

The above formulated heterogeneous model of the emulsion zone, consisting of PDEs (Equations (42) and (43)) with BCs (Equations (44) and (45)) and ICs (Equations (53) and (54)), and ODEs (Equations (48) and (49)) with ICs (Equation (55)), hereinafter is referred to as full heterogeneous model or model H1 (Table 3). In the following subsection, selected modifications of the full heterogeneous model are formulated based on some simplifications resulting from the omission of internal resistances to mass and heat transfer.

2.4.2. Heterogeneous Model with Simplified Models of the Pellet (Models H2 and H3)

The temperature gradient within a catalyst pellet is often negligible, thus in process calculations, it is usually assumed that the intraparticle temperature distribution is uniform. Neglecting the temperature distribution leads to the transformation of the PDE (36) into time-dependent ODE yielding the lumped-thermal model (referred to as model H2, Table 3).

Assuming that the resistance to heat transport is concentrated in the boundary layer around the particle, and at the same time the temperature gradient inside the particle is negligibly small, the heat balance equation of the particle for the lumped-thermal model becomes:

$$V_p{\rho }_p{c}_s\frac{d{T}^p}{dt}=A_p{\alpha }_q(T^e-T^p)+V_p{r}_{\mathrm{A}\upsilon }(C_{\mathrm{A}}^p,T^p)\cdot (-\Delta h),$$

(56)

where $r_{\mathrm{A}\upsilon }(C_{\mathrm{A}}^p,T^p)$ is an overall process rate, per unit volume of the catalyst, defined as:

$$r_{\mathrm{A}\upsilon }=\frac{1}{V_p}{\displaystyle \underset{V_p}{\int }{r}_{\mathrm{A}}(C^p(r),T^p)}d\upsilon ,$$

(57)

where $d\upsilon =4\pi r^{2}dr=4\pi L_p^3{\zeta }^{2}d\zeta$.

After the introduction of the dimensionless concentration of reactant A, ${\beta }_{\mathrm{A}}$, Equation (56) takes the following form:

$$\frac{d{T}^p}{dt}={a^{″}}_1(T^e-T^p)+{\delta }^{\prime }\cdot R_{\mathrm{A}\upsilon }({\beta }_{\mathrm{A}},T^p),$$

(58)

with initial condition:

$$T^p(t=0)=T_0^p,$$

(59)

where:

$${a^{″}}_1=\frac{A_p{\alpha }_q}{V_p{\rho }_p{c}_s}=\frac{6{\alpha }_q}{d_p{\rho }_p{c}_s},{\delta }^{\prime }=\frac{C_{\mathrm{A},ref}(-\Delta h)}{{\rho }_p{c}_s}.$$

(60)

Assuming additionally that both internal and external resistance to mass transfer is significant, the concentration profile of A in the pellet is quantitatively described, as in the full heterogeneous model, by PDE (Equation (42)) with corresponding BCs defined by Equation (44). Balance equations of the gas in the emulsion are also identical to the equations formulated for the full heterogeneous model, i.e., Equations (48) and (49), while functions defined by Equations (19) and (20) describe the distributions of the state variables in the bubbles.

The third heterogeneous model tested in this work is the so-called lumped-particle model (model H3, Table 3). It was formulated based on the assumption that the resistance to mass and heat transport is centered in the boundary layer around the particle, while both concentration and temperature gradients inside the catalyst particle are negligible. As in the case of the model H1 and model H2, the concentration of reactant A and temperature of the gas in the bubbles and interstitial emulsion gas are described by Equations (19), (20), (48) and (49). The mass and heat balance equations of a single catalyst pellet for the lumped-particle model can be written as:

$$V_p\left({\epsilon }_p\frac{d{C}_{\mathrm{A}}^p}{dt}+{\rho }_p\frac{d{q}_{\mathrm{A}}}{dt}\right)={\mathrm{A}}_p{k}_m(C_{\mathrm{A}}^e-C_{\mathrm{A}}^p)-V_p{r}_{\mathrm{A}}(C^p,T^p),$$

(61)

$$V_p{\rho }_p{c}_s\frac{d{T}^p}{dt}=A_p{\alpha }_q(T^e-T^p)+V_p{r}_{\mathrm{A}}(C^p,T^p)\cdot (-\Delta h),$$

(62)

where $C^p$ and $T^p$ are lumped variables characterizing the entire volume of the catalyst pellet. After making use of the dimensionless variables, Equations (61) and (62) take the form:

$$\frac{d{\beta }_{\mathrm{A}}}{dt}={a^{″}}_2(X_{\mathrm{A}}^e-{\beta }_{\mathrm{A}})-{a^{″}}_3{R}_{\mathrm{A}}({\beta }_{\mathrm{A}},T^p),$$

(63)

$$\frac{d{T}^p}{dt}={a^{″}}_1(T^e-T^p)+{\delta }^{\prime }\cdot R_{\mathrm{A}}({\beta }_{\mathrm{A}},T^p),$$

(64)

where ${a^{″}}_1$ is defined by Equation (60) and:

$${a^{″}}_2=\frac{A_p{k}_m}{V_p({\epsilon }_p+{\rho }_p{K}_a)},{a^{″}}_3=\frac{1}{{\epsilon }_p+{\rho }_p{K}_a}.$$

(65)

Equations (63) and (64) must be supplemented with ICs determining the initial values of concentration and temperature in the entire volume of the catalyst pellet, namely:

$${\beta }_{\mathrm{A}}(t=0)={\beta }_{\mathrm{A}0},T^p(t=0)=T_0^p.$$

(66)

2.5. Model Parameters and Numerical Methods

To compare qualitatively and quantitatively the steady-state and transient properties of the bubbling fluidized-bed catalytic reactor determined using different mathematical models, some representative values of the physicochemical and operating parameters were selected based on the analysis of several solid-catalyzed chemical processes of industrial importance [1,5,33,34] and adsorption isotherms from the literature [35,36]. The values of the model parameters used in the numerical simulations are given in

Table 4. Model parameters used in numerical simulations.

Parameter	Value	Unit	Parameter	Value	Unit
aqkq	0	kW·m−3·K−1	P	101,325	Pa
cg	1.0	kJ·kg−1·K−1	S	1	m2
cs	0.8	kJ·kg−1·K−1	Sh	2	-
dp	2 · 10−4	m	Tf	300	K
DbA = DeA	2 · 10−5	m2·s−1	yAf	0.1	-
Def	2 · 10−6	m2·s−1	yPf	0	-
E	7 · 104	kJ·kmol−1	εp	0.5	-
Ea	3 · 104	kJ·kmol−1	λef	10−4	kW·m−1·K−1
Δh	−2 · 105 to −6 · 105	kJ·kmol−1	λg	2 · 10−5	kW·m−1·K−1
Hmf	1.0	m	μg	2.6 · 10−5	Pa·s
k0	106, 5 · 106	s−1	ρg	0.7	kg·m−3
Ka0	5 · 10−5	m3·kg−1	ρp	1600	kg·m−3
Nu	2	-

. To simplify the analysis, it was assumed that the reactor is operated adiabatically, i.e., the product of the heat exchange area per unit volume of the fluidized bed, a_q, and the overall wall heat transfer coefficient, k_q was set to zero. However, some results concerning the influence of the heat removal from the fluidized bed onto the reactor steady-state behavior can be found in [26].

The influence of the following parameters, expected to affect strongly the transport resistances, onto the steady-state and transient properties of the reactor was investigated:

• fluidization ratio, $l_f=u_0/u_{mf}$, determining the mean residence time of gas and hydrodynamic conditions in the reactor;
• enthalpy of the chemical reaction, Δh, and frequency coefficient in the Arrhenius equation, k₀;
• initial temperature of the fluidized bed, $T_0$ (in case of transient simulations only).

In the steady-state analysis, solution branches with respect to l_f were obtained by parametric continuation of the equations describing the steady-state operation of the reactor, equations derived by setting to zero the time derivatives in the dynamic models P, H1, H2 and H3 (Table 3). A local parametrization algorithm [37] implemented in FORTRAN was employed to calculate successive points of the steady-state branches. In the local parameterization method, consecutive points of the solution branch are determined by numerical resolution of the so-called extended system of nonlinear algebraic equations, consisting of the model equations and an additional parameterizing equation [37]. At each step of the continuation procedure, the model equations were solved using Newton’s method, which in case of models H1 and H2 was combined with the shooting algorithm and the fourth-order Runge-Kutta method for the integration (with a step size Δζ = 0.01) of the differential equations describing the concentration and temperature profiles within the catalyst pellet.

For the transient analysis, the dynamic equations were integrated in time using numerical codes developed in FORTRAN and employing a variable-coefficient ordinary differential equation solver (DVODE) [38] based on implicit backward formulas for numerical differentiation. In the case of models H1 and H2, the PDEs describing the distribution of the state variables within the particles were previously transformed into ODEs by the method of lines with N = 101 discrete nodes along particle radius [39]. Depending of the model (P or H1-H3) the number of ODEs to be integrated in time was from two (for model P) to 204 (for model H1, with two ODEs describing the gas phase and 202 ODEs describing the discretized catalyst pellet). It is worth underlining that the computational time of the integration of the most complex model, i.e., the so-called full heterogeneous model H1, over t = 10⁵ s was of the order of 1 s. The computational time of the integration of less complex and thus lower-dimensional models P, H2 and H3 was about 0.2–0.5 s.

3. Results and Discussion

Figure 3 shows the steady-state branches of $X_{\mathrm{A}}^e$ and T^e determined with respect to the fluidization ratio, l_f, for several values of the chemical reaction enthalpy, $\Delta h$, obtained using a pseudohomogeneous model P and heterogeneous models: H1, H2 and H3. Solid lines indicate stable steady states, whereas dashed lines indicate unstable steady states.

Regardless of the adopted value of $\Delta h$ and l_f, the differences between the values of $X_{\mathrm{A}}^e$ and $T^e$ determined according to different mathematical models of the reactor are indistinguishable in the scale of Figure 3. Despite the rather complex steady-state characteristics of the examined system, which consists of multiple solutions, the results obtained suggest that the steady-state operation of the analyzed catalytic BFBR can be described using even the simplest model with the assumption of pseudohomogeneity made for the emulsion zone (model P). It is worth underlining here that the complex structure of the steady-state diagrams is not surprising: such kind, and even more complex behavior (e.g., isolated solution branches) was observed also experimentally both in fixed-bed and fluidized-bed catalytic reactors [40,41].

Representative values of the state variables derived from the different models collected in

Table 5. Representative values of state variables corresponding to the high-conversion steady state obtained using different models for $\Delta h=-6\cdot {10}^5$ kJ·kmol⁻¹, $k_0=5\cdot {10}^6$ s⁻¹ and $T_f=300$ K.

State Variable	${\mathit{\beta }}_{\mathbf{A}}(\mathbf{0})$	${\mathit{\beta }}_{\mathbf{A}}(\mathbf{1})$$\mathbf{or}{{\mathit{\beta }}_{\mathbf{A}}}^1$	${\mathit{X}}_{\mathbf{A}}^{\mathit{e}}$	${\mathit{T}}^{\mathit{p}}(0),\mathbf{K}$	${\mathit{T}}^{\mathit{p}}(1)$$\mathbf{or}{{\mathit{T}}^{\mathit{p}}}^{\mathbf{2}}\mathbf{,}\mathbf{K}$	${\mathit{T}}^{\mathit{e}}\mathbf{,}\mathbf{K}$	${\mathit{W}}_{\mathbf{P}}\mathbf{,}\mathbf{kmol}\mathbf{·}{\mathbf{m}}^{-\mathbf{2}}{\mathbf{s}}^{-\mathbf{1}}$
Model	lf = 2 (point P1 in Figure 3)
P	–	–	5.808 · 10−6	–	–	1145.88	4.312 · 10−5
H1	1.501 · 10−6	1.034 · 10−5	1.346 · 10−5	1145.87	1145.87	1145.87	4.312 · 10−5
H2	1.501 · 10−6	1.034 · 10−5	1.346 · 10−5	–	1145.87	1145.87	4.312 · 10−5
H3	–	5.808 · 10−6	8.926 · 10−6	–	1145.88	1145.87	4.312 · 10−5
Model	lf = 10 (point P2 in Figure 3)
P	–	–	1.825 · 10−3	–	–	686.05	9.839 · 10−5
H1	1.803 · 10−3	1.838 · 10−3	1.846 · 10−3	686.06	686.06	686.04	9.839 · 10−5
H2	1.803 · 10−3	1.838 · 10−3	1.846 · 10−3	–	686.06	686.04	9.839 · 10−5
H3	–	1.824 · 10−3	1.831 · 10−3		686.06	686.04	9.839 · 10−5

¹ In case of lumped-parameter model of the catalyst pellet. ² In case of lumped-thermal or lumped-parameter model of the catalyst pellet.

confirm these visual observations. It can be seen that treating the emulsion zone as a pseudohomogeneous medium (model P) results in the values of the emulsion temperature, $T^e$, to be nearly identical to those calculated using the most complex model, i.e., the so-called full heterogeneous model H1, both for $l_f=2$ and $l_f=10$. Moreover, minimal differences between temperatures $T^p(0)$, $T^p(1)$ and $T^e$ obtained from model H1 indicate that both external and internal resistance to heat transfer is negligibly small. Nevertheless, some discrepancies are observed between the values of dimensionless concentration of reactant A in the emulsion zone, $X_{\mathrm{A}}^e$, determined with the aid of models P and H3 (i.e., lumped particle model) and the models accounting for the non-uniform distribution of the concentration in the pellet (i.e., models H1 and H2). These differences are much higher at lower values of the fluidization ratio and are due to the presence of the significant gradients of concentration both inside the catalyst pellet and at the solid-gas interface: for instance at $l_f=2$ the ratio between the surface concentration of A, ${\beta }_A(1)$, and ${\beta }_A(0)$ is as high as 6.89, whereas the ratio between the concentration of A in the interstitial emulsion gas, $X_{\mathrm{A}}^e$, and ${\beta }_{\mathrm{A}}(1)$ is around 1.30 (applies to the values calculated from models H1 and H2, Table 5). However, the discrepancies between the values of $X_{\mathrm{A}}^e$ calculated from models of different complexity, and the presence of intraparticle gradients of ${\beta }_{\mathrm{A}}(\zeta )$ (in the case of model H1 and H2) do not influence the value of the reactor yield with respect to product P, W_P, defined as:

$$W_{\mathrm{P}}=\frac{{\nu }_{\mathrm{P}}}{\left|{\nu }_{\mathrm{A}}\right|}{u}_0{C}_{\mathrm{A}f}\cdot {\overline{\alpha }}_{\mathrm{A}}(1)=\frac{{\nu }_{\mathrm{P}}}{\left|{\nu }_{\mathrm{A}}\right|}{u}_{mf}{C}_{\mathrm{A},ref}\frac{T_{ref}}{T^e}\cdot {\overline{\alpha }}_{\mathrm{A}}(1),$$

(67)

where:

$${\overline{\alpha }}_{\mathrm{A}}(1)=a_6{\alpha }_{\mathrm{A}}^b(1)+(1-a_6){\alpha }_{\mathrm{A}}^e,a_6=\frac{\delta u_b}{\delta u_b+(1-\delta ){\epsilon }_{mf}{u}_e}.$$

(68)

The values of W_P at the reactor outlet resulting from models P and H1-H3 (Table 5) are identical up to four digits, and this confirms the capability of the simplified model P to provide a correct prediction of the steady-state productivity of the reactor.

As already mentioned, the resistance to internal and external mass transfer increases when the fluidization ratio, l_f, determining the mean residence time of the gas in the apparatus, decreases (Figure 4). Furthermore, as expected, both internal and solid-gas interface concentration gradients get higher as the enthalpy of chemical reaction, $\Delta h$, increases.

Another parameter that strongly influences the steady-state characteristics of the reactor and determines the importance of resistance to mass transport is the frequency coefficient, k₀. As can be seen in Figure 5a, again the values of the concentration of A in the interstitial emulsion gas, $X_{\mathrm{A}}^e$, obtained using models characterized by different complexity, are almost identical, and the branches of steady states determined with respect to l_f again practically overlap. However, according to the full heterogeneous model, for faster chemical reaction and for longer mean residence time of the gas in the reactor, i.e., at lower values of l_f, a significant amount of component A reacts in the narrow outer shell of the pellet which is manifested by a very low concentration of A at the pellet center (Figure 5b). As before, the temperature distributions (not reported here; applies to model H1) that develop within the pellet are practically uniform, regardless of the adopted value of k₀; also, the temperature gradient at the solid-gas interface is negligibly small.

The above results of numerical simulations were confronted with criteria used in practice for assessing the significance of internal and external transport resistances. For the m^th order chemical reaction described by power law kinetics, Weisz and Prater [42] proposed the following criterion to assess the significance of resistance to internal mass diffusion:

$$\frac{m+1}{2}\frac{r_{\mathrm{A}\upsilon }{L}_p^2}{D_{ef}{C}_{\mathrm{A}s}}<<1,$$

(69)

In the above inequality (Equation (69)), $C_{\mathrm{A}s}$ denotes the concentration of component A at the pellet surface, i.e., $C_{\mathrm{A}s}=C_{\mathrm{A}}^p(L_p)$. If the inequality is satisfied, it means that the process is controlled by chemical kinetics and that intraparticle mass diffusion does not play a significant role.

To evaluate the presence of temperature gradients in the pellet, Anderson [43] formulated the following criterion:

$$\frac{4}{3}\frac{(-\mathsf{\Delta }h)r_{\mathrm{A}\upsilon }{L}_p^{2}E}{{\lambda }_{ef}R{T}_{bulk}^2}<1.$$

(70)

Again, fulfilling the inequality (Equation (70)), referring to the spherical particle, means that the temperature field is practically uniform in the entire volume of the catalyst pellet.

The criteria that permit to assess the significance of interphase concentration and temperature gradients proposed by Mears [44] are as follows:

$$\frac{r_{\mathrm{A}\upsilon }{L}_{p}m}{0.15\cdot k_m{C}_{\mathrm{A},bulk}}<1,$$

(71)

$$\frac{(-\mathsf{\Delta }h)r_{\mathrm{A}\upsilon }{L}_{p}E}{0.15\cdot {\alpha }_{q}R{T}_{bulk}^2}<1.$$

(72)

The importance of diffusion and reaction limitations is also commonly evaluated in terms of the internal, η, and the overall, η₀, effectiveness factor of the catalyst [45]:

$$\eta =\frac{r_{\mathrm{A}\upsilon }}{r_{\mathrm{A}s}},{\eta }_0=\frac{r_{\mathrm{A}\upsilon }}{r_{\mathrm{A},bulk}},$$

(73)

where $r_{\mathrm{A}s}$ is the rate of chemical reaction evaluated at the particle surface conditions, $r_{\mathrm{A},bulk}$ is the rate of chemical reaction evaluated at the bulk conditions, i.e., at the interstitial emulsion gas conditions, whereas $r_{\mathrm{A}\upsilon }$ is the overall reaction rate defined by Equation (57).

Figure 6 shows the values of the moduli specifying the significance of resistance to mass and heat transport, and the values of internal and overall effectiveness factor calculated based on the results obtained from model H1 for the same parameters as the steady-state branches presented in Figure 3. The values of the moduli characterizing the mass transport resistances (MTRs) (Figure 6a) confirm the previous observations: both external and internal MTRs become higher when l_f decreases and $\mathsf{\Delta }h$ increases. Regardless of the thermal effect of the process, both the values of moduli characterizing internal and external heat transport resistances (HTRs) are much lower than one, i.e., inequalities given in Equations (70) and (72) are fulfilled, which is also in line with the previous observations.

Analysis of the internal, η, and of the overall, η₀, effectiveness factor of the catalyst pellet shown in Figure 6c also confirms the very strong effect of external and internal mass transport on the overall process rate at lower values of l_f and higher values of $\mathsf{\Delta }h$. Moreover, due to negligible gradients of the temperature both within the pellet and at the solid-gas interface (representative values are given in Table 5) both the values of η and η₀ do not exceed unity, even when the reactor is operated at relatively high l_f.

Despite the occurrence of significant resistance to mass transport and large intraparticle gradients of the concentration, it can be concluded that it is sufficient to use the pseudohomogeneous model of the emulsion zone when analyzing the steady-state operation of the reactor. However, the application of the more complex heterogeneous model may be necessary to simulate the transient behavior of the reactor. The transient analysis consisted of numerical simulations of the reactor start-up and was conducted in order to emphasize the differences that emerge as a results of application of different mathematical models of the emulsion zone. It was motivated by the findings concerning highly exothermic processes conducted in fixed-bed reactors, for which it was demonstrated that pseudohomogeneous models may lead to erroneous evaluation of the transient behavior of the reactor [46].

Figure 7a,b show, respectively, the time trajectories of $X_{\mathrm{A}}^e$ and $T^e$ during the reactor start-up calculated using models P and H1-H3. In all simulations it was assumed that the bed was preheated to the initial temperature of 450 K, while the initial concentration of A was set to zero.

In the light of the results obtained for steady-state conditions, the results of the dynamic simulations are quite surprising: while for $l_f=2$ the time trajectories of the concentration (Figure 7a) and temperature (Figure 7b) obtained using all models practically overlap, the decrease of the residence time of reacting gas in the reactor ($l_f=5$) results in significant discrepancies between the trajectories of $X_{\mathrm{A}}^e$ and $T^e$ calculated using the models accounting for the intraparticle and interphase concentration gradients (i.e., models H1 and H2) and the models in which the intraparticle or both the intraparticle and interphase concentration gradients were neglected (models H3 and P, respectively). This phenomenon is all the more surprising, because as in the case of steady-state behavior, the instantaneous values of the moduli characterizing the importance of MTRs shown in Figure 8a suggest that the resistance to mass transport is less significant at higher values of the fluidization ratio. Therefore, one would instead expect some discrepancies in the results at lower values of l_f.

The differences between the results obtained, respectively, using models H1 and H2 and models H3 and P, may be explained by the steep gradients of the state variables (Figure 7), also reflected in the MTR and HTR curves (Figure 8) around $t\approx 8\cdot {10}^3\mathrm{s}$ for $l_f=5$ and resulting from a sudden ignition of the reacting mixture (Figure 7b). At higher value of l_f one initially observes a temporary reduction of the emulsion temperature (Figure 7b) due to the short residence time of the gas in the reactor, resulting in a very low global process rate. This leads to a significant increase of the concentration of reactant A (Figure 7a), followed by sudden ignition and rapid heat release.

In some cases, the simplest models of the reactor, i.e., models H3 and P, not only fail to predict the time of ignition correctly, but they fail completely in the prediction of such a phenomenon. As shown in Figure 9, for $l_f=5.06$ both models P and H3 are unable to predict the sudden drop of reactant concentration (Figure 9a) due to rapid ignition (Figure 9b). On the other hand, model H2, that is the model that takes into account the nonuniform distribution of the reactant in the pellet but assumes uniform intraparticle temperature, faithfully reflects—both qualitatively and quantitatively—the solution obtained from model H3.

The nature of ignition is strictly connected to the occurrence of multiple steady states in the analyzed system, and in particular to basins of attraction of the stable steady states. Figure 10 shows the influence of the fluidization ratio, l_f, and the start-up temperature, T₀, onto the final steady state. This is represented graphically in terms of so-called reduced basins of attraction [47]. Unstable solution branches $T^e$ from Figure 3 are plotted in Figure 10 (denoted by MU), to compare their loci with the position of the basin’s boundaries. The boundaries (also called separatrices, denoted in Figure 10 by Sep.) of the upper (US) and lower (LS) steady state were determined using the algorithm proposed in Reference [47] using models P and H1, and under the assumption ${\beta }_{\mathrm{A}}(\zeta )=X_{\mathrm{A}}^b(z)=X_{\mathrm{A}}^e=X_{\mathrm{A}0}$ and $T^p(\zeta )=T^b(z)=T^e=T_0$ at $t=0$. Although the unstable solution branches of the emulsion temperature $T^e$ obtained using model H1 and model P practically overlap, the separatrices obtained from model P are always located above those determined using the model H1. The differences between the results derived using these two models increase, even up to a few Kelvins, as l_f and $\Delta h$ increase. All pairs of the fluidization ratio, l_f, and reactor start-up temperature, T₀, located in the narrow area in between Sep. H1 and Sep. P (Figure 10b) correspond to the situation in which the pseudohomogeneous model P fails to predict correctly the reactor behavior, i.e., while the trajectories determined in this parameter range using model H1 converge to the upper stable steady state (i.e., ignited state), model P predicts reactor extinction.

4. Conclusions

The commonly used assumption regarding the pseudohomogeneous nature of the emulsion zone in a catalytic BFBR was verified both for steady-state and transient conditions using four different models of the reactor, including the model with the pseudohomogeneity assumption made for the emulsion zone and three models accounting for the heterogeneous character of this zone.

Even in the case of significant resistance to mass transfer, the pseudohomogeneous model of the emulsion zone proved to be sufficiently accurate for the prediction of the reactor steady state. However, the selection of the appropriate mathematical model is of major importance in case of transient simulation. Adoption of the pseudohomogeneous model of the emulsion zone to simulate chemical processes characterized by relatively high thermal effect, when run at higher values of fluidization ratio may lead to erroneous results. While it is possible to assume a uniform distribution of the temperature within the pellet, assumption of uniform intraparticle concentration for such processes may results in a wrong prediction of the reactor behavior, especially when simulating reactor start-up. In the region of occurrence of multiple steady states, both the pseudohomogeneous (P) and lumped-particle model (H3) fail to predict correctly the sharp time gradients of the state variables related to a sudden ignition characterizing the process run at high l_f and started-up with T₀ slightly higher than the temperature corresponding to the boundary of domain of attraction of the ignited steady state.

Comparison of the results of the numerical simulation with the criteria used in practice to evaluate the significance of internal and external transport resistances demonstrates that empirical criteria are not a sufficient tool to assess the significance of transport resistance. Although the values of the moduli characterizing transport resistances are in line with the observations made based on the parametric analysis of steady states and of the transient behavior of the apparatus, the consistency is only of qualitative, not quantitative character. Thus, considering the results presented in this work, it is inappropriate to adopt a priori the pseudohomogeneous model for the emulsion zone. Such an assumption, especially in case of highly exothermic processes, must be preceded by a detailed analysis made with the aid of a full heterogeneous model.

The modelling approach presented here may be extended to analyze real processes of industrial importance. However, it has to be borne in mind that, when analyzing the influence of internal and external mass transport resistances, particular attention must be paid to a proper description of the diffusive mass transport with the aid of mathematical models suitable for the description of multicomponent gas mixtures.