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.
  1. Mathematical models and their solutions. Verification of a method precision.
  2. Packed Bed Reactor Model (PBRM)
The model is formulated with following assumptions:
  1. Process is isothermal.
  2. Flow rate of mobile phase is constant.
  3. Column is packed with catalyst with the same radius Rp.
  4. The radius of dead zone in catalyst is Rd and can be function of time and position inside column.
  5. Changes of physico – chemical parameters in radial direction a negligible.
  6. 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).