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.