TMA681 Tillämpad matematik K, 2001

STUDIO 6–7. THE CATALYST PELLET

Theory: S. Larsson, “Kompletterande föreläsningsanteckningar och övningar” (K). R. Aris, “The Mathematical Theory of Diffusion and Reaction in Permeable Catalysts. Vol. 1. The Theory of the Steady State”, Oxford University Press, 1975, Chapter 3.

1. Introduction

A substance in gas form is diffusing into a porous catalyst pellet, where it reacts in a first order reaction under isothermal and stationary conditions. This leads to the reaction-diffusion equation,

see (K 2.3),

(1)  \begin{aligned}−∇ · (D∇c) & = −kc, in V , \\ c & = c_{f}, on S, \end{aligned}

where c [mol/m^{3}] is the concentration of the interesting substance in the catalyst pellet, c_{f} is the (constant) ambient concentration, D [m^{2}/s] is the diffusion coefficient, k [s^{−1}] is the rate coefficient of the reaction, L [m] is the diameter of the pellet, and ∇c = (\frac{∂c}{∂x}, \frac{∂c}{∂y} , \frac{∂c}{∂z}).

PSfrag replacements

c
c_{f}
V
S
L

Figure 1. The catalyst pellet.

Note that we make the simplifying assumptions that the system has reached a stationary state and that the temperature is constant. The latter assumption is the same as saying that we ignore the temperature dependence in the rate coefficient k = k_{0} exp(−E/(RT)) given by Arrhenius’ law. However, we want to allow the difffusion and catalyzation properties to be different in different parts of the pellet. This means that D, k_{0}, and E may vary with position, so that D = D(x, y, z), k = k(x, y, z). Exercise 1. Introduce dimensionless variables. Let c_{ref} = c_{f}, D_{ref}, k_{ref}, L_{ref} = L be reference constants and set

(x^{∗}, y^{∗}, z^{∗}) = (x/L, y/L, z/L), u = c/c_{f}, k^{∗} = k/k_{ref}, D^{∗} = D/D_{ref}.

Show that (1) becomes

(2)  \begin{aligned}−∇^{∗} · (D^{∗}∇^{∗}u) & = −ϕ^{2}k^{∗}u, in V^{∗}, \\ u & = 1, on S^{∗}, \end{aligned}

where ∇^{∗} = (\frac{∂}{∂x^{∗}} , \frac{∂}{∂y^{∗}} , \frac{∂}{∂z^{∗}}) and the dimensionless number ϕ = L\sqrt{k_{ref}/D_{ref}} is the Thiele modulus. Hint: the chain rule gives \frac{∂}{∂x} = \frac{1}{L} \frac{∂}{∂x^{∗}}.

If the diffusion and rate coefficients are constant, then, of course, we take D_{ref} = D, k_{ref} = k, so that D^{∗}(x^{∗}, y^{∗}, z^{∗}) ≡ 1, k^{∗}(x^{∗}, y^{∗}, z^{∗}) ≡ 1, and (2) becomes

(3)  \begin{aligned}−Δ^{∗}u & = −ϕ^{2}u, in V^{∗}, \\ u & = 1, on S^{∗}, \end{aligned}

Date: May 7, 2001, Niklas Ericsson, Stig Larsson, Mathematics, Chalmers University of Technology.

2 STUDIO 6–7. THE CATALYST PELLET

where Δ^{∗} = \frac{∂^{2}}{∂x^{∗2}} + \frac{∂^{2}}{∂y^{∗2}} + \frac{∂^{2}}{∂z^{∗2}} is the Laplace operator^{1}.

PSfrag replacements

c c_{f} V S L u 1 V^{∗} S^{∗} 1

Figure 2. The pellet in dimensionless variables.

The concentration, and hence the reaction rate, will be lower in the middle of the pellet. This means that the catalyzing power of the pellet is not fully used. The efficiency of the catalyst pellet is measured by the quotient of the actual total reaction rate and the ideal reaction rate that would be achieved if the concentration and the rate coefficient were everywhere equal to the their reference values, i.e., if c = c_{f}, k = k_{ref}. This quotient is called the effectiveness factor and is given by

(4)  η = \frac{∫_{V} kc dV}{k_{ref}c_{f}V} = \frac{1}{V^{∗}} ∫_{V^{∗}} k^{∗}u dV^{∗}.

Show this equality! Homework 1. Show that η can also be expressed in terms of the flux through the boundary of the pellet: (n^{∗} is the exterior unit normal vector)

(5)  η = \frac{1}{ϕ^{2}} \frac{1}{V^{∗}} ∫_{S^{∗}} n^{∗} · (D^{∗}∇^{∗}u) dS^{∗}.

Hint: use (2) and Gauss’ divergence theorem.

PSfrag replacements

c c_{f} V S L u 1 V^{∗} S^{∗}

1 −\frac{1}{2} \frac{1}{2}

0 r

Figure 3. Slab and cylinder.

2. An exact solution

If the coefficients are constant and the geometry of the pellet is simple, namely, a slab, a circular cylinder, or a sphere, then we can solve (3) analytically. We study one of these examples here. Exercise 2. If the pellet is a slab (or a cylinder with sealed mantle surface), Figure 3, then u depends only on the axial coordinate r = x^{∗} = x/L and (3) becomes

(6)  \begin{matrix}− \frac{d^{2}u}{dr^{2}} = −ϕ^{2}u, −\frac{1}{2} < r < \frac{1}{2}, \\ u(−\frac{1}{2}) = u(\frac{1}{2}) = 1. \end{matrix}

1In Studio 2 we mentioned several reasons for writing the mathematical model in dimensionless form. Now we can add one more reason: the teacher does not have to bother to find realistic values for the physical parameters D, k_{0}, E. Since they are conveniently grouped into the dimensionless number ϕ, we can just say that we want to solve the equation for small, medium and large values of ϕ.

Solve this boundary value problem and compute η (analytically by hand). Plot η as a function of ϕ (with Matlab).

Hint: recall the hyperbolic sine and cosine from Studio 2. Answer: u(r) = cosh(ϕr)/ cosh(ϕ/2),

c(x) = c_{f} cosh(\sqrt{\frac{k}{D}} x)/ cosh(\sqrt{\frac{k}{D}} \frac{L}{2}), η = ∫_{−1/2}^{1/2} u(r) dr = tanh(ϕ/2)/(ϕ/2), u(0) = 1/ cosh(ϕ/2). Note that u(r) ≡ 1 and hence η = 1, if ϕ = 0.

3. Numerical solution

In general we cannot solve (2) analytically. If we assume that the pellet is a long straight cylinder with an arbitrary cross-section, so that u depends only on two variables x^{∗}, y^{∗}, then we can use Matlab’s PDE Toolbox to compute approximate solutions. Exercise 3. Solve (2) with D^{∗}(x^{∗}, y^{∗}, z^{∗}) ≡ 1, k^{∗}(x^{∗}, y^{∗}, z^{∗}) ≡ 1, and ϕ = 1 for several crosssections: circle, rectangle, triangle, or whatever. Compute η for each geometry. Make sure that the diameter is the same = 1 in each case. Which one has the highest efficiency?

Hint: Start the GUI (Graphical User Interface) ’pdetool’ by typing

>> pdetool

on the Matlab command line^{2}. When you start ’pdetool’ you will enter the draw mode. In draw mode, you define Ω, i.e., the geometry of the problem. (In our case Ω is the cross-section of V^{∗}.) To simplify the drawing, select Grid (to turn on the grid) and Snap (to enable the “snap-to-grid” feature) from the Options menu. The “snap-to-grid” feature helps aligning the objects you draw. In this example^{3}, we take Ω to be the square: Ω = \{(x, y) : 0 < x, y < 1/\sqrt{2}\}. From the same, i.e., the Options, menu, select Grid Spacing..., to open the Grid Spacing dialog box. Turn off both the Auto checkboxes. Change the default linear grid spacing to −0.2 : 0.2 : 1, and add an extra tick at 1/\sqrt{2}, for both the x-axis and the y-axis. Press the Apply push button, and end the Grid Spacing dialog by pressing the Done button. Also, change the axes limits by selecting Axes Limits..., from the Options menu. Set both the x-axis range and the y-axis range to [−0.2 1], press the Apply button and finally the Close button. To draw Ω, just click on the “rectangle icon” (below the File menu), put the cursor at the origin and click-and-drag using the left mouse button to create the square.

Next, enter boundary mode, by selecting it from the top of the Boundary menu. In this mode, we define the boundary conditions of the problem. Different boundary conditions may be specified for different boundary segments. You select a boundary segment by clicking on it. You can select more than one segment by, repeatedly, shift-clicking on them. In this example, we want to apply the same boundary condition, u = 1, to all (four) boundary segments, so, select them all. Then, double-click on any one of them, to open the Boundary Condition dialog box. To set the desired boundary condition, change the (default) value of r from 0 to 1, and press the OK button.

To complete the definition of the problem, enter pde mode, by selecting it from the top of the PDE menu. In this mode, we specify the coefficients of the equation. Select the square, Ω, by clicking on it, and then double-click on it, to open the PDE Specification dialog box. To specify the equation, change the value of a from 0.0 to 1 (note that a corresponds to ϕ^{2}k^{∗}) and the value of f from 10.0 to 0. Finally press the OK button. We have now defined the problem.

To compute an approximate solution to the problem, the PDE Toolbox utilizes the Finite Element Method (FEM). The idea is to divide Ω into a finite number of elements (triangles, in this case), and look for an approximate solution that is a (linear) polynomial on each triangle.

To divide Ω into triangles, enter mesh mode, by selecting it from the top of the Mesh menu (alternatively, just click on the “triangle icon” (below the Mesh menu)). If you want, you can refine the triangulation (mesh), by clicking on the icon with four triangles. This will give a more accurate solution, but, it will also require more computation and, accordingly, take a longer time.

2It is also possible to call toolbox functions from the command line (or M-files), but we will only use the GUI.

3To learn more, see e.g. the “Getting Started” and “Using the Graphical User Interface” sections in Chapter 1:

“Tutorial”, of the “Partial Differential Equations (PDE) Toolbox User’s Guide”. You find this guide by first typing » helpdesk on the Matlab command line, then click on Online Manuals (in PDF), and finally click on the link Partial Differential Equations (PDE) Toolbox User’s Guide.

To solve the problem, select Solve PDE, from the Solve menu (or, click on the “= icon”). When the computation has finished, the solution is automatically plotted. This can also be

done by selecting Plot Solution, from the Plot menu. Try different visualization options. You can specify these by selecting Parameters..., also from the Plot menu, which opens the Plot Selection dialog box. Experiment!

Finally, to compute the effectiveness factor, η, you may use the function effectiveness.m, which you can obtain at the address,

http://www.math.chalmers.se/~stig/underv/K/effectiveness.m

This function utilizes numerical integration, quadrature, to compute η. The integration is performed over one triangle at a time, and the contributions from all triangles are added together. Since Ω is the union of all triangles, this gives the integral over Ω. As input data, it needs (recall the definition of η) the solution u, and the coefficient a. To do the element-wise integration, it also needs some information about the triangulation, which is contained in two matrices called p and t. Before calling effectiveness.m, you must therefore export these data from ’pdetool’ to the “usual” Matlab workspace. To export u, select Export Solution..., from the Solve menu, and press the OK button. Analogously, export a from the PDE menu, and p, t from the Mesh menu. (Some additional data, that we don’t need, will also be exported.) To compute η, now type:

>> phi = 1; % The Thiele modulus
>> eta = effectiveness(u, a, p, t, phi)

Exercise 4. Try some other variants:

(a) Dead core: a region in the middle of the pellet has very small diffusion coefficient.

Hint: Consider the domain Ω = \{(x, y) : \sqrt{x^{2} + y^{2}} < \frac{1}{2}\}, with core Ω_{core} = \{(x, y) : \sqrt{x^{2} + y^{2}} < \frac{1}{5}\}. In draw mode, first draw a circle with centre at the origin and radius \frac{1}{2}, then draw another circle, also with centre at the origin, and radius \frac{1}{5}. In boundary mode, as before, apply the Dirichlet boundary condition, u = 1, on the boundary of Ω: \{(x, y) : \sqrt{x^{2} + y^{2}} = \frac{1}{2}\}. In PDE mode, select Show Subdomain Labels, from the PDE menu. Note that Ω_{core} has one number, say 1, and Ω − Ω_{core} has another number, say 2. We call Ω_{core} and Ω − Ω_{core} subdomains of Ω. The union of the subdomains is Ω itself. Now, different expressions for the PDE coefficients may be specified on the different subdomains. Analogously to defining boundary conditions, you select a subdomain by clicking on it. Then, double-click on it to open the PDE Specification dialog box. Apply, e.g., c = 0.1 on Ω_{core}, and c = 1 on Ω − Ω_{core} (note that c corresponds to D^{∗}). Now, you can continue as above. It might be a good idea to make one extra refinement of the triangulation in this case.

(b) Poisoning: a region near the boundary of the pellet has very small rate coefficient, k^{∗} << 1 or k^{∗} = 0.

(c) Hole: there is a hole in the middle of the pellet, where the diffusion coefficient is very large, D^{∗} >> 1, and the reaction rate is zero, k^{∗} = 0.

(d) Robin’s boundary condition: replace the boundary condition in (1) by Robin’s boundary condition (K (36), (38)). Show that in dimensionless variables this becomes

n^{∗} · (D^{∗}∇^{∗}u) + ν(u − 1) = 0, on S^{∗}.

Express the dimensionless number ν (the Biot number for mass transfer) in the original variables. Answer: ν = \frac{Ld}{D_{ref}} .

Compute solutions with the PDE Toolbox.

Hint: In the Boundary Condition dialog box, choose Neumann as condition type (“Neumann”, here, means a “generalized Neumann” boundary condition, which we call a Robin boundary condition). Now, just note that, in our setting, both g and q correspond to ν.