2.2.1 Numerical method of solution of the model for catalyst particle.
For several simple models of reaction mechanism, the Eq. (18) with boundary conditions Eqs. (16-17b) were solved analytically[1,9,25]. The analytical solutions will be in the following used to control the accuracy of numerical solutions. The typical method of numerical solution is based on discretization of spatial derivatives with one of the “classical” methods, e.g. finite difference, finite element, orthogonal collocation. Next obtained set of algebraic equations is solved typically with Newton-Rapson method. Alternative is “shoot method”[25] in which the concentration y in boundary condition (16) should be selected by trial-and-error method to fulfill condition (17a).
In this work instead of solving Eq. (18) we solved equivalent model consisted of two first order differential equations.
Let’s \(Y_{1}=y_{p}\) and \(Y_{2}\) is a derivative of yp along radius of the pellet.
Then:
\(Y_{2}=\frac{dY_{1}}{\text{dR}}\) (19)
and according to Eq. (18) we obtain:
\(\left(\frac{\alpha\left(R_{p}-R_{d}\right)}{R_{d}+\left(R_{p}-R_{d}\right)R}Y_{2}+\frac{\partial Y_{2}}{\partial R}\right)\frac{R_{p}^{2}}{(R_{p}-R_{d})^{2}}=(1-\varepsilon_{p})\cdot\Phi^{2}(Y_{1})^{n}\)(20)
Set of Equations (19) and (20) have to be solved with initial conditions equivalent to boundary conditions of model Eq. (18).
We assumed that:
for R=0
\(Y_{1}=\epsilon\mathrm{\text{\ \ \ \ and\ \ }}\mathrm{Y}_{2}=0\)(21)
where \(\epsilon\) is unit roundoff.
The approximation of concentration Y1 with very low value (computer’s “zero”) is indispensable to start integration what would be impossible assuming Y1=Y2=0.
More over the condition
\(\frac{R_{p}}{(R_{p}-R_{d})}Y_{2}(R=1)=\text{Bi}\cdot(Y_{b}-Y_{1}(R=1))\)(22)
where Yb is a dimensionless concentration in the bulk phase, have to be fulfilled.
In summary, for example for calculation of position Rdof dead zone for given value of Yb and Φ we should:
Hardware and software platform as well as additional algorithms applied are as follows.
As the ODE solver the CVODE (https://computation.llnl.gov/projects/sundials) can be recommended.
The position of Rd can be easily found using well known bisection method, remembering that 0 <= Rd < 1.
Time of finding the Rd position with proposed method is equal about 20 [ms] on PC computer with Intel® i7-8700K CPU, 3.7GHz. The calculation was done declaring double precision for real variables.
The double or extended precision is a highest precision which offers commercial C, C++, C# and others compilers. This precision is enough for obtaining solution of above discussed problem with accuracy greater than 6 significant digits (comparing analytical and numerical solution), when Thiele modulus is not to large.
Precision of the method we will verify on the basis of the mentioned above analytical solutions for power type kinetic equation. If exponent in kinetic equation was not close to unity, precision of calculation was very high. For example, critical value of Thiele modulus,\(\Phi_{\text{cr}}\) (position Rd=0) for n = - 0.6 in Eq. (20), Bi = 1 and spherical geometry is equal 0.87659746 according to numerical solution and 0.87659754 from analytical solution (here and in the following it was assumed that kinetic rate is referenced to volume of the pellet). However, if parameter n is sufficiently close to 1 the problems with precision appears. When n = 0.9 then Φ = 17.169654 and 17.600041 for numerical and analytical solutions, respectively. When n = 0.99 then Φ = 58.08 and 195.25 for numerical and analytical solutions respectively.
To explain the source of this problem let’s analyze the analytical expression for concentration inside slab pellet[26]
\begin{equation} {y=\left[1-\frac{\Phi}{\Phi_{\text{cr}}}\left(1-z\right)\right]}^{2/(1-n)}\nonumber \\ \end{equation}
where
\begin{equation} \Phi_{\text{cr}}=\sqrt{\frac{2\left(1+n\right)}{\left(1-n\right)^{2}}}\nonumber \\ \end{equation}
It is clear that when parameter n is increasing to 1 the concentration drops from concentration on the surface of the slab to almost zero at very short distance. When this distance become comparable with unit roundoff the numerical calculation has to fails. The remedy is to increase the number of significant digits in representation of real number. This is possible applying another software platform e.g. Maple. The Maple enables solution of the above discussed problem with any number of significant digits. If a number 64 is accepted as a number of significant digits we obtain much better approximation the critical value of Thiele modulus, namely Φcr= 194.32 for parameter n = 0.99.