2.2.2 Non isothermal catalyst particle.
The numerical method applied for isothermal particle can be easily adopted for not isothermal process, especially for single not reversible reaction.
According to Arrhenius low the kinetic of n-order reaction is:
\(r=k\exp\left(-\frac{E}{\text{RT}}\right)C^{n}\) (23)
where: E is activation energy, R is gas constant, and T is temperature.
The dimensionless, steady state mass and heat transfer model can be written as follows:
\(\left(\frac{\alpha\left(R_{p}-R_{d}\right)}{R_{d}+\left(R_{p}-R_{d}\right)R}\frac{\partial y_{p}}{\partial R}+\frac{\partial^{2}y_{p}}{\partial R^{2}}\right)\frac{R_{p}^{2}}{\left(R_{p}-R_{d}\right)^{2}}=\left(1-\varepsilon_{p}\right)\cdot\Phi^{2}\cdot\text{Re}\left(y_{p},\theta\right)\)(24)
\(\left(\frac{\alpha(R_{p}-R_{d})}{R_{d}+(R_{p}-R_{d})R}\frac{\partial\theta}{\partial R}+\frac{\partial^{2}\theta}{\partial R^{2}}\right)\frac{R_{p}^{2}}{{(R_{p}-R_{d})}^{2}}=\left(1-\varepsilon_{p}\right)\cdot\Phi^{2}\cdot\beta\cdot Re\left(y_{p},\theta\right)\)(25)
and boundary conditions for mass transfer are:
for R=1
\(\frac{R_{p}}{(R_{p}-R_{d})}\frac{\partial y_{p}(R)}{\partial R}=Bi_{m}\cdot[y_{b}-y_{p}(R)\)] (26)
for R=0
\(y_{p}(R)=0\mathrm{\text{\ \ \ \ and\ \ }}\frac{\partial y_{p}(R)}{\partial R}=0\)(27)
and for heat transfer
for R=1
\(\frac{R_{p}}{(R_{p}-R_{d})}\frac{\partial\theta}{\partial R}=Bi_{T}\cdot[\theta_{b}-\theta(R)\)] (28)
for R=0
\(\frac{\partial\theta(R)}{\partial R}=0\) (29)
The kinetic term can be rewritten as follows:
\(Re(y_{p},\theta)=y_{p}^{n}\cdot exp(\gamma*(1-\frac{1}{\theta}))\)(30)
where: \(\theta=T/T_{r}\) , Tr is reference temperature, \(\beta=\text{ΔH}D_{\text{eff}}y_{b}/(\lambda T_{r})\) is Prater parameter, yb is a concentration in a bulk phase,\(\text{Bi}_{T}=R_{p}\beta/\lambda\), \(\lambda\) is a heat conductivity, ΔH is enthalpy of reaction and\(\gamma=E/(RT_{r})\).
The system of two differential equations can be simplified to one differential equations for mass transfer. Namely it can be shown that following relationship between dimensionless temperature \(\theta\), concentration in the particle yp and concentration on the surface of particle ys, (\(y_{s}=y_{p}\left(R=1\right))\) exist[9]:
\(\theta=1+\beta*((1-y_{s})\frac{\text{Bi}_{m}}{\text{Bi}_{T}}+y_{s}-\ y_{p})\)(31)
The above model can be solved with technique presented before for isothermal process. However, this time estimation of critical value of Thiele modulus is more complicated because two parameters have to be simultaneously estimated: Thiele modulus and unknown value of surface concentration ys. This estimation can be done effectively with Levenberg – Marquardt method[27]. Time of estimation takes only a fraction of second.
The accuracy of numerical calculation as before systematically decrease with increase of Thiele modulus above about 10.
The numerical inaccuracies resulted from the mentioned in the previous section types of errors are characteristic not only for our method. The numerical solution of the discussed here problem was also presented in fundamental paper [9] devoted for dead zone occurrence – see Fig. 5 in the cited paper. The author solved non-isothermal model by the numerical integration in the MATLAB package. Unfortunately, for parameter n increasing to 1 the solution is not robust. For example, for β=0, Bim→∞ and n=0.8 the critical Thiele modulus is equal 10.49 according to analytical solution. On other hand from Fig. 5 and line 7 on this Figure follows that the critical Thiele modulus is equal about 14. The most likely, analytical and numerical solutions vary as a result of not enough significant digits in representation of real number in the MATLAB. There is no information on precision of calculations in [9].
Effectiveness factor versus Thiele modulus.
The catalyst pellet work is frequently analyzed on the base on dependency of effectiveness factor versus Thiele modulus. The effectiveness factor was calculated from the formula:
\(\eta=\frac{(\alpha+1)}{\Phi^{2}y_{s}^{n}}\)\(\left.\ \frac{\text{dy}}{\text{dR}}\right\rfloor_{R=1}\) (32)
With proposed numerical method the calculation of curves presented in Fig. 1 and Fig. 2 takes several seconds.
We will show some results (Fig. 1, Fig. 2a, b) that confirm robustness of the developed algorithm. All the next calculations were performed for sphere assuming Rp=1.
In Fig. 1 is presented dependency between effectiveness factor, ƞ and Thiele modulus, \(\Phi\) calculated analytically[28] and numerically according to proposed method for slab geometry. The both solutions overlap each other, what confirm the accuracy of tested numerical method.
[CHART]
Fig. 1. Dependency between effectiveness factor, ƞ and Thiele modulus,\(\Phi\) calculated analytically (solid line) and numerically (circles) for slab geometry.
More complex dependencies are expected for n < 0. Our method is able to detect irregular phenomena as multiple steady-states, including chaotic solutions – Fig. 2. Between Φcr and\(\Phi^{\prime}\) the multiple steady states are observed (see also[28]). Zoom of the region between Φ=1.50 and 2.00 shows that the structure of multiple steady-state region is very complicated. The visible in Fig. 2 irregular forms are successively repeated.
a)
[CHART]
b)
[CHART]
Fig. 2. Dependency between effectiveness factor \(\eta\) versus\(\Phi\). Calculation performed for sphere, n=-0.6, Bi=10000, Φcr=1.6769; a) full range of Φ and b) magnified region 1.55 .. 1.75.
Further analysis of irregular solutions is omitted, it is not the aim of the present work.
Numerical solution of PBRM.
In chemical reactor also can appears a dead zone, of course only when conditions of forming dead zone in pellet are fulfilled.
In the following will be presented a selected simulations for PBRM described with set of equations (10-16, 17a). However, at the beginning let’s assume that dead zone does not exist (Rd =0), so the PBRM consist of Eqs. (10-16, 17b). The model was solved with version of Orthogonal Collocation on Finite Element (OCFE) method described in[29]. After approximation of spatial derivatives with appropriate algebraic formulae according to OCFE a set of Ordinal Differential Equation (ODE) is obtained. This set of equation was solved with CVODE (https://computation.llnl.gov/projects/sundials) program with relative accuracy of calculation equal 1e-8 and absolute equal 1e-10. In the following 3 subdomains (finite elements) were used for catalyst particle and 50 – 150 for reactor. In all cases the number of collocation points inside subdomain was equal three.
When dead zone in catalyst particle can exists then its position Rd should be calculated. The position of Rd for any cross-section of reactor and any time of integration of PBRM was obtained from independent solution of Eq. (11) with method described above in section 2.2.1. The position Rd was found with bisection method for given bulk concentration \(y(\tau,x\)).
All performed in the following calculation were done reference reaction rate to volume of catalyst matrix – the kinetic term was multiplied by\((1-\varepsilon_{p})\). It means that the value\((1-\varepsilon_{p})F^{2}\) is equal to Thiele modulus, \(F^{2}\), when reaction rate is reference to volume of catalyst. It was also assumed Rp=1, Cr=1.
The accuracy of proposed numerical solution of PBRM was tested assuming stationary state, negligible mass transfer resistances and dispersion, inlet concentration \(y_{f}=1\) and Rp=1. After mathematical manipulation analogous to that presented in[30,31] the PBRM simplifies to equation:
\(\frac{\partial y}{\partial x}=-\frac{\left(1-\varepsilon_{e}\right)}{\varepsilon_{e}}\cdot\left(1-\varepsilon_{p}\right)\frac{\Phi^{2}\text{St}}{\left(\alpha+1\right)\text{Bi}}y^{n}\)(34)
The solution of Eq. (34) is:
\(y\left(x\right)=\left\{\left[\frac{1}{1-n}-\frac{\left(1-\varepsilon_{e}\right)}{\varepsilon_{e}}\cdot\left(1-\varepsilon_{p}\right)\frac{\Phi^{2}\text{St}}{\left(\alpha+1\right)\text{Bi}}x\right]\left(1-n\right)\right\}^{\frac{1}{1-n}}\)\((35)\)
The selected numerical solutions of time dependent PBRM for n=-0.6, 0, 0.6 and different St, Bi and Pe are presented and compared with analytical expression (Eq. (35). The other parameters were equal:\(\varepsilon_{e}=0.4\), \(\varepsilon_{p}=0.5\), \(\Phi=0.2\). Comparison of numerical solution of full model of the reactor and its simplified version is presented in Fig. 3. It confirms precision of the developed method.
[CHART]
Fig. 3. Dependency between dimensionless concentration, y, and dimensionless distance along column, x, for n = 0. The solid line represent solution of PBRM for St=3e5, Bi=1e3, Pe=1e7. The dash line represents Eq. (35).
Solving the PBRM the initial and inlet, dimensionless concentration in reactor was assumed to be equal 1.0. In Fig. 3 dependency of concentration along column is presented for n=0. The value of concentration equal zero indicates the beginning of dead zone inside column, in this case at x=1/3. The solution of PBRM is identical with analytical expression Eq. (35) at the moment when concentration decrease to zero, what means that stationary states was reached. Continuation of calculation with PBRM has not sense because negative “concentration” in the area x>1/3 appears as indicate Eq. (35).
In Figs. 4 the concentrations decrease along column for n = 0.6 or n = - 0.6 are presented.
a)
[CHART]
b)
[CHART]
Fig. 4. Dependency between dimensionless concentration, y, and dimensionless distance along column, x, a) for n = 0.6 and b) for n=-0.6. The solid line represents solution of PBRM for: St=3e5, Bi=1e3, Pe=1e7; dash-dot line for St=30, Bi=0.1, Pe=1e7; and dot line for St=30, Bi=0.1, Pe=10. The dash line represents Eq. (35).
For lack of mass transfer resistances and dispersion (we assumed St/Bi = 300, St=3e5, Bi=1e3, Pe=1e7, for larger values of these parameters non changes in concentration profiles was observed for presented calculations) the numerical solutions overlap the analytical. However, the decrease of concentration towards zero is in both cases different. In the first case concentration slowly approximate to zero when in second case it abruptly achieves zero. These phenomena are connected with reaction kinetic rate. In first case the kinetic rate decreases proportionally to yn, it means slower then decrease of concentration. Whereas, in second case the kinetic rate increases to infinity when concentration approaches to zero – kinetic rate is proportional to 1/yn. It has also consequences in attempt to numerically solve the problems with negative order of kinetic rate. The numerical procedure does not converge near the boundary of dead zone. When the problem is solved using time depended model (as it was done in this paper) the time step of the integration of the ODE is decreasing dramatically when concentration in some part of catalyst particle approaches to zero. It is indicator that practically steady state condition was reached and calculation can be stopped what was done for simulations presented in Fig. 4b.
In Fig. 4 the results for simulation with not negligible mass transfer resistances and dispersion are included. As should be expected the dead zone area in reactor is decreasing with increasing of mass transfer resistance and/or dispersion.
As it was mentioned at the beginning of this section the PBRM was solved with OCFE method with controlling the position of dead zone in catalyst pellet. We did several numerical experiments comparing the results of simulation with switched on and switched off the control of the position of the dead zone inside pellet. In all cases the obtained positions of the beginning of dead zone in reactor were practically the same. However, the computations without controlling of dead zone position inside pellet can be much faster and are recommended.