Abstract
For two processes of large importance, catalysis and biocatalysis, were
reported zones without reactants, so called dead zone (DZ). They results
from diffusional transport limitations, when apparent reaction order is
between (-1..1). Formation of DZ reduces effectiveness of catalyst and
influence packed bed reactor productivity. For simple reaction kinetic
model, a DZ width inside a pellet can be calculated analytically solving
appropriate differential mass balance model. However, generally the
analytical solution is unknown and only with using numerical method the
position of DZ can be established.
The problem with DZ appearance belongs to problems with moving
boundaries. Its solution requires application of special numerical
procedure and relatively long CPU time. In this work it was proposed a
simple, very fast numerical method for calculation of DZ position inside
pellet. The method proposed combined with orthogonal collocation on
finite elements can be applied for analysis of work of packed bed
reactor.
Keywords : dead zone; catalyst particle; packed bed reactor;
numerical algorithms.
Introduction
Heterogeneous catalytic processes have for years been of crucial
importance in the chemical industry. Nowadays a biocatalitic processes
gain also in importance. For both types of the processes the existence
of a zones without key reactant\souts around center of symmetry for
catalyst particle or inside biofilm were reported. It means that the key
component was exhausted before reaching center of particle.
The dead zone can appear in real processes relatively often. Lack of
substrata in the zone decreases potential reactor or bioreactor
productivity and influences the process economy.
When dead zone develops an internal part of catalyst pellet does not
work and the part of catalyst is wasted. In a consequence, a designed
reactor volume would be oversized.
The existence of zone without reactants was predicted by
Aris,[1] and Temkin,[2] for
heterogeneous catalysis. Next many authors theoretically analyzed
problems of existence and formation of dead core Diaz and
Hernandez,[3], Bandle and
Stackgold,[4], Bobisud,[5].
Different generation terms were tested Bobisud and
Royalty,[6], Fedotov,[7],
Garcia-Ochoa and Romero,[8] and it led to
formulation conditions which have to be satisfied for dead zone
formation for selected kinetic equations (so-called necessary conditions
of dead zone formation). The sufficient conditions of dead zone
formation for the simplest kinetic equation were presented
in[8]. Extended results of this topic were
presented also by Andreev,[9]. The theoretical
prediction of dead zone formation was confirmed experimentally both for
heterogeneous catalysis processes (e.g. for methanol steam reforming
over Cu/ZnO/Al2O3 catalyst Lee et al,[10], for
hydrogenation of benzene over nickel-alumina catalyst Jiracek et
al,[11], for acetic acid oxidation Levec et
al,[12] and Look and
Smith,13], for hydrogenation of propylene reaction
carried out on commercial Ni catalyst pellets Szukiewicz et
al,[14]) and bioprocesses (for cephalosporin C
production processes Araujo et al,[15], Cruz et
al,[16], for process of Penicillin G enzymatic
hydrolysis Cascaval et al,[17], for process of
3-chloro-1,2-propanediol degradation Konti, et
al,[18], in an anaerobic fixed-bed reactor Zaiat
et al,[19] in catalytic particles containing
immobilized enzymes for the Michaelis-Menten kinetics Pereira and
Oliveira,[20] ).
The probability of dead zone appearance is higher if catalyst activity
increases[1]. Recently, a lot of effort are put
into improving activity of catalysts, in most cases by adding promoters.
It is such a widely discussed topic that it would be difficult for to
present it here in more detail. In only one scientific journal (Applied
Catalysis B: Environmental) this topic is mentioned in 113 research
works in 2018 and 95 in the current year. For this reason, we present
few research and review articles only to indicate fields of interest and
application of promoters. E.g. reverse water-gas shift reaction was
analyzed in Zhang et al,[21]; authors show that Cu
as a promoter boost the CO2 conversion due to higher
active sites concentration. Cs is an excellent dopant to enhance CO
selectivity, especially at low temperatures. CO2methanation process was discussed by Quindimil et
al,[22]; there were explained that impregnation of
La2O3 provides basicity and enhances Ni
dispersion, significantly improving the
CO2-to-CH4 activity. Mukherjee et
al,[23] investigated importance of development of
alkali metal amide/imide as catalyst for ammonia decomposition.
Promoters are also applied for catalysts of biomass conversion Pang et
al,[24]; there were found that transition metal
carbides showed unique activity in hydrodeoxygenation reactions and the
Ni promoter lowered the carburization temperature and protected carbides
surface from phase changes.
The analytical expression for calculation of a real size of a dead zone
inside a catalyst pellet or the real depth of penetration reagents in a
biofilm can be obtained for relatively simple reaction kinetic. However,
generally those calculations are possible only with the use of numerical
methods. It should be noted that solution for the mathematical model of
dead zone problem is a harder task than the regular model solving
(regular model = model without dead zone, commonly applied and discussed
in literature). From a mathematical point of view the regular model is
typical boundary value problem while the dead zone model belongs to the
free boundary value problems. The solution of the free boundary problem
(see problem formulation in the next section) needs to use an iterative
procedure to establish position of the boundaries, and for each
iteration step the mass balance for catalyst pellet have to be solved.
Additionally, many of popular algorithms of solution of differential
equation fail. Tests of few popular algorithms employed to the dead-zone
problem solution have been presented in the source already
cited,[20] . Troubles reported are: divergence,
low rate of convergence, high complexity of programming. The recommended
by authors,[20] orthogonal collocation method on
finite elements is in their opinion time-consuming and complex to code
in a programming platform. Summarizing the considerations concerning
complexity of free boundary problem solution, we conclude that time
necessary can be much longer than those needed for regular model
solution. And follows that developing of fast, new algorithm be very
helpful for further investigations in the mentioned areas. One more
problem, that was not discussed in [20] appears.
The problems concerns numerical precision of calculations. Mainly
truncation and round-off errors can be a source large errors. We will
show the problem and discuss it as well as present a solution in the
sections 2.2.1 and 2.2.2.
In the previous our work[14] the dead zone
investigations were model-based and verified experimentally. The
model-based method consists of two stages: (i) experimental
determination of kinetic equation and some process parameters (ii)
numerical solution of the model (we use a numerical method described
in[25]. The efficiency of applied numerical method
was similar to that from literature[20] however, a
CPU time of solution was rather long.
In this work we make some effort to improve numerical side of
investigations and develop simple and very fast algorithm that can find
solution of a problem with moving boundary. The method presentation is
accompanied by verification of calculations precision on the basis of
analytic solution of the proper problems simplified.
The natural enhancement of the problem of the dead zone inside a
catalyst pellet is model of reactor packed with catalyst bed where dead
zone in the pellets can be formed or not (operating conditions in the
reactor determine dead zone formation, so dead zones can be present in
the pellets in the part of reactor, while in another part cannot). This
task is also realized by us for heterogeneous fixed bed reactor and
analysis of results obtained were done. We also present simple
verification of usefulness of the presented method basis on our own
results of experiments.
Note that an application of the method presented to bioreactors requires
only simple modification of the equations that describe mobile phase,
mathematical approach remains the same.
- Mathematical models and their solutions. Verification of a
method precision.
- Packed Bed Reactor Model (PBRM)
The model is formulated with following assumptions:
- Process is isothermal.
- Flow rate of mobile phase is constant.
- Column is packed with catalyst with the same radius
Rp.
- The radius of dead zone in catalyst is Rd and can be
function of time and position inside column.
- Changes of physico – chemical parameters in radial direction a
negligible.
- Dispersion coefficient (DL) is constant.
On the base of this assumptions the two-mass balance equation for each
component are formulated. First for mobile phase percolating through a
bed and second for fluid inside catalyst pores. The model is coupled
with appropriate initial and boundary conditions and a model of reaction
kinetic.
Mass balance for component in mobile phase:
\(\varepsilon_{e}\cdot\frac{\partial C}{\partial t}+u\cdot\frac{\partial C}{\partial z}=\varepsilon_{e}D_{L}\cdot\frac{\partial^{2}C}{\partial z^{2}}-(1-\varepsilon_{e})k_{\exp}\frac{\left(\alpha+1\right)}{R_{p}}\cdot\left[C-C_{p}\left(r=R_{p}\right)\right]\)(1)
Mass balance for component in catalyst pellet:
\(\varepsilon_{p}\cdot\frac{\partial C_{p}}{\partial t}=D_{\text{eff}}\cdot\frac{1}{r^{\alpha}}\cdot\frac{\partial}{\partial r}(r^{\alpha}\frac{\partial C_{p}}{\partial r})-(1-\varepsilon_{p})\cdot k_{\text{rea}}{(C_{p})}^{n}\)(2)
In above equations the C is concentration of component in mobile phase,
Cp is concentration in catalyst pores, u is superficial
velocity, DL is dispersion coefficient,
kext is external mass transfer coefficient, r is radial
coordinate, z is axial coordinate, t is time, Rp is
radius of the catalyst pellet, Deff is effective
diffusion coefficient, krea is reaction kinetic rate, n
is apparent reaction order, α is particle geometry: spherical,
cylindrical, slab (2,1,0 respectively), εe is bulk
porosity, and εp is particle porosity.
The Eqs. (1) and (2) have to be coupled with initial and boundary
conditions.
Initial conditions:
for 0<z<L and 0
<r<Rp
\(C\left(0,z\right)=C^{o}\) (3)
\(C_{p}(0,r,z)=C_{p}^{o}(r,z)\) (4)
Boundary conditions for Eq. (1):
for t>0 and z=0
\(u_{f}C_{f}-u(0)C(0)=-\varepsilon_{e}D_{L}\cdot\frac{\partial C}{\partial z}\)(5)
for t >0 and z=L
\(\frac{\partial C}{\partial z}=0\ \) (6)
Here L can be regarded as total column length or the position at which
the dead zone in column begins.
Boundary conditions for Eq. (2):
for t >0 and r=Rp
\(D_{\text{eff}}\cdot\frac{\partial C_{p}(t,r)}{\partial r}=k_{\text{ext}}\cdot[C-C_{p}(t,r)]\)(7)
for t>0 and r=Rd
\(C_{p}(t,r)=0\mathrm{\text{\ \ \ and\ }}\frac{\partial C_{p}(t,r)}{\partial r}=0\)(8a)
when Rd >= 0
or
\(\frac{\partial C_{p}(t,r)}{\partial r}=0\) (8b)
when dead zone does not exist.
The meaning of the symbols is: L - reactor length, Rd -
radius of dead zone, Cf - concentration of component at
the inlet to the column, uf - superficial velocity at
the inlet to the column.
We introduce following dimensionless parameters:
\(x=\frac{z}{L}\) \(\tau=\frac{\text{tu}}{L\varepsilon_{e}}\)\(R=\frac{r-R_{d}}{R_{p}-R_{d}}\)
\(y=\frac{C}{C_{r}}\) \(y_{p}=\frac{C_{p}}{C_{r}}\)\(D_{\text{eff}}=\frac{D_{m}}{\chi}\varepsilon_{p}\) (9)
\(\text{Pe}=\frac{\text{uL}}{D_{L}\varepsilon_{e}}\)\(St=\frac{k_{\text{ext}}a_{p}L\varepsilon_{e}}{u}\)\(\text{Bi}=\frac{k_{\text{ext}}R_{p}}{D_{\text{eff}}}\),\(\Phi^{2}=\frac{R_{p}^{2}k_{\text{rea}}}{D_{\text{eff}}}{(C_{r})}^{n-1}\)
where: Cr is reference concentration, Bi is Biota
number, Pe is Peclet number, St is Stanton number and Φ is Thiele
modulus. In this work we assume Cr = 1.
Equations (1) and (2) of mass balance with dimensionless parameters can
be rewritten as follows:
\(\frac{\partial y}{\partial\tau}+\frac{\partial y}{\partial x}=\frac{1}{\text{Pe}}\cdot\frac{\partial^{2}y}{\partial x^{2}}-\frac{(1-\varepsilon_{e})}{\varepsilon_{e}}\cdot\text{St}\cdot[y-y_{p}(R=1)]\)(10)
\(\frac{\partial y_{p}}{\partial\tau}=\frac{1}{\varepsilon_{p}}\frac{St}{\left(\alpha+1\right)Bi}\cdot\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}}-\frac{\left(1-\varepsilon_{p}\right)}{\varepsilon_{p}}\cdot\frac{\Phi^{2}St}{\left(\alpha+1\right)Bi}\left(y_{p}\right)^{n}\)(11)
Initial conditions for \(\tau=0\):
\(y(0,x)=y^{o}\) (12)
\(y_{p}(0,R,x)=y_{p}^{0}(R,x)\) (13)
Boundary conditions for Eq. (10):
For \(\tau>0\); x=0
\(y_{f}-y(\tau,0)=-\frac{1}{\text{Pe}}\cdot\frac{\partial y(\tau,0)}{\partial x}\),
(14)
For \(\tau>0\); x=1
\(\frac{\partial y(\tau,1)}{\partial x}=0\) (15)
When dead zone in reactor appears than Eq. (15) is fulfilled at the
beginning of reactor dead zone.
Boundary conditions for Eq. (11):
for \(\tau>0\); R=1
\(\frac{R_{p}}{(R_{p}-R_{d})}\frac{\partial y_{p}(\tau,R)}{\partial R}=Bi\cdot(y-y_{p}\left(\tau,R\right))\)(16)
for \(\tau>0\); R=0
\(y_{p}(\tau,R)=0\mathrm{\text{\ \ \ \ and\ \ }}\frac{\partial y_{p}(\tau,R)}{\partial R}=0\)(17a)
when dead zone exists
or
\(\frac{\partial y_{p}(\tau,R)}{\partial R}=0\) (17b)
when dead zone does not exist.
Please notice that for given definition of dimensionless radius R the R
value change between 0 and 1 despite the fact that the radius of dead
zone can change with decreasing substrate concentration along column.
Replacing coordinate r with its dimensionless equivalent R simplify
numerical solution of reactor model.
Before analysis of the work of the reactor it will be convenient to
analyze the behavior catalyst particle in stationary state conditions.
2.2 Dimensionless model of reaction in catalyst
particle at steady-state.
For stationary state conditions Eqs. (11) can be rewritten as follows:
\(\left(\frac{\alpha(R_{p}-R_{d})}{R_{d}+(R_{p}-R_{d})R}\frac{\partial y_{p}}{\partial R}+\frac{\partial^{2}y_{p}}{\partial R^{2}}\right)\frac{R_{p}^{2}}{{(R_{p}-R_{d})}^{2}}=(1-\varepsilon_{p})\cdot\Phi^{2}{(y_{p})}^{n}\)(18)
The boundary conditions are expressed by Eqs. (16 – 17b).