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:
- assume Rd position;
- solve model Eqs. (19-22) with ODE (Ordinal Differential Equation)
solver
- check if condition Eq. (36) is fulfilled. If not the new value of
Rd should be chosen and calculation have to be
repeated.
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.