Question 6.PS.1: Reaction–Diffusion in a Packed Bed We consider a packed bed ...
Reaction–Diffusion in a Packed Bed
We consider a packed bed reactor, shown in Fig. 6.12. At the left edge of the bed, x = 0, the concentration of the reactant is equal to the bulk concentration, c = c_0. There is a wall at the right edge of the bed, x = l, with the corresponding no-flux boundary condition,
\frac{dc}{dx}|_{x=l}= 0 (6.4.1)
For a first-order reaction in the bed, the mass balance takes the form
0 = D\frac{∂^2c}{ ∂x^2} − k(x)c (6.4.2)
We assume that, due to the size of the particles, the reaction rate is spatially dependent. We adopt a simple model, where the reaction is at full strength in the center half of the particle and zero otherwise. For example, in the first particle of the bed, the reaction rate constant is
k(x)= \begin{cases} 0 &for& x ∈ [0, 0.05l]\\ k_0 &for &x ∈ (0.05l, 0.15l)\\ 0 &for &x ∈ [0.15l, 0.20l] \end{cases} (6.4.3)
This periodic function repeats five times over the length l of the array.
The goal of this case study is to use finite differences to see how the concentration profile inside the reactor depends on the Thiele modulus.
If we are going to work in terms of a Thiele modulus, we should convert the problem to dimensionless form,
\tilde{x} = x/l, \tilde{c} = c/c_0, \tilde{k} = k/k_0 (6.4.4)
With this choice of dimensionless variables, you can convert the unsteady equations into
0 = \frac{d^2\tilde{c}}{d\tilde{x}^2 }− \tilde{k}(\tilde{x})φ^2\tilde{c}. (6.4.5)
The dimensionless reactivity has the same spatial dependence as eq. (6.4.3) in each particle, but now the maximum value is \tilde{k} (\tilde{x} ) = 1 instead of k_0. For a first-order reaction, the Thiele modulus is
φ =\sqrt{\frac{k_0l^2}{D} } (6.4.6)
The equation for the interior nodes is
c_{i+1} + c_{i−1} −\left[2 + φ^2 Δx^2k_i\right] c_i = 0 (6.4.7)
where
Δx =\frac{1}{n − 1} (6.4.8)
and k_i is the reaction rate constant evaluated at node i. The left node equation is
c_1 = 1 (6.4.9)
The fictitious node equation at the right boundary gives
\frac{c_{n+1} − c_{n−1}}{2Δx} = 0 (6.4.10)
so c_{n+1} = c_{n−1}. Using this result in the discretized form of the differential equation gives
2c_{n−1} −\left[2 + φ^2 Δx^2k_n\right] c_n = 0 (6.4.11)
Since the reaction rate constant k_n is zero at the boundary, this result simplifies to
c_{n−1} − c_n = 0 (6.4.12)
Let’s write first a function that solves the differential equation for a given value of the Thiele modulus and the number of nodes n, in case we want to do some mesh refinement later on.
1 function c = linear_react(n,phi2)
2 dx = 1/(n-1);
3 k = zeros(n,1);
4 for i = 1:n
5 x = mod((i-1)*dx,0.2);
6 if x >= 0.05 && x <= 0.15
7 k(i) = 1;
8 end
9 end
10 A = zeros(n);
11 b = zeros(n,1);
12 A(1,1) = 1;
13 b(1,1) = 1;
14 for i = 2:n-1
15 A(i,i-1) = 1;
16 A(i,i) = -2 - dx^2*phi2*k(i);
17 A(i,i+1) = 1;
18 end
19 A(n,n-1) = 1;
20 A(n,n) = -1;
21 A = sparse(A);
22 c = A\b;
The beginning of this program sets up the reaction rate vector using the modulo operator. We have already seen the modulo operator when we set the output frequency for integrating IVPs, for example in Program 4.1. Here, we are using the modulo operator to map the position in the “global” domain x ∈[0, 1] back onto a “local” domain x ∈ [0, 0.2]. The for loop starting in line 4 loops through the n nodes. The value of x for that node is x_i = (i − 1) Δx. When we compute the modulus of x_i with 0.2 in line 5, we get the remainder after dividing by 0.2, which maps the result back to the first catalyst particle. For example, if x_i = 0.55, the operation x=mod(0.55,0.2) returns x = 0.15, since you can divide 0.2 into 0.55 twice with a remainder of 0.15. Line 6 then checks to see if this position is inside the catalyst particle. We need to have an “and” operator && to be sure that we are to the right of the start of the particle and to the left of the end of the particle. If this is the case, the reaction rate constant is set to 1. If not, then the reaction rate constant is zero because we initialized the vector k to be all zeros in line 3.
Aside from setting the reaction rate vector, this program should look remarkably similar to the one we used for Example 6.2, since the underlying differential equations are the same. Lines 12 and 13 of the present program implement the Dirichlet boundary condition on the left, while lines 19 and 20 implement the Neumann boundary condition on the right. For the interior nodes, we need the local reaction rate k_i, which is called in line 16. Finally, we solve this tridiagonal system using a banded solver in lines 21 and 22.
Figure 6.13 shows the concentration profile in this model of a packed bed for different values of the Thiele modulus. Clearly, as the Thiele modulus increases and reaction becomes dominant over diffusion, there is a sharp reduction in concentration away from the reservoir at x = 1. If you look closely, you may also be able to see the discrete nature of the reaction, since it only takes place inside the particles. You may have thought that the concentration profile would be more step-like with this reactivity profile. However, remember that diffusion takes place in both the catalyst pellets and the interstitial space. Although the reaction is only occurring inside the particles, diffusion both inside the particles and in the space between particles is smearing out the discreteness of the reaction rate.
