Question 7.PS.1: 2D Reaction–Diffusion In this last case study, we will consi...
2D Reaction–Diffusion
In this last case study, we will consider a basic question about reaction–diffusion problems. We will look at the solution to the steady state reaction–diffusion equation
D∇^2c = R(c) (7.5.1)
in a box of size L × L, where R(c) is the rate of consumption of the reactant. We will consider either a first-order reaction,
R = k_1c (7.5.2)
or a second-order reaction
R = k_2c^2 (7.5.3)
The boundaries are fixed at the bulk concentration c_0.
Our goal is to determine how the minimum concentration inside the box depends on the Thiele modulus for these two reaction models – which reaction type gives us the smaller value of the concentration?
To answer this question, we make the equations dimensionless with a length scale L and concentration c_0. For the linear problem, this gives
\frac{d^2c}{dx^2} +\frac{d^2c}{dy^2} − φ^2c = 0 (7.5.4)
with the Thiele modulus
φ =\sqrt{\frac{k_1L^2}{D}} (7.5.5)
For the nonlinear case, the same dimensionless form gives
\frac{d^2c}{dx^2} +\frac{d^2c}{dy^2} − φ^2c^2 = 0 (7.5.6)
with
φ =\sqrt{\frac{k_2c_0L^2}{D} } (7.5.7)
In both cases, the Thiele modulus measures the relative rate of reaction to diffusion.
We will discretize the domain into an n × n grid, using an equally spaced grid where Δx = Δy =Δ . For the linear problem, the interior nodes are of the form
\frac{c_{i+1,j} − 2c_{i,j} + c_{i−1,j}}{Δ^2} +\frac{ c_{i,j+1} − 2c_{i,j} + c_{i,j−1}}{Δ^2 }− φ^2c_{i,j} = 0 (7.5.8)
while the residual for the nonlinear problem is
\frac{c_{i+1,j} − 2c_{i,j} + c_{i−1,j}}{Δ^2} +\frac{ c_{i,j+1} − 2c_{i,j} + c_{i,j−1}}{Δ^2 }− φ^2c_{i,j}^2 = 0 (7.5.9)
In both cases, it proves convenient to multiply through by the grid spacing. If we also switch to the global variable k = i + (j − 1)n the linear problem becomes
c_{k+1} − (4 + φ^2 Δ^2)c_k + c_{k−1} + c_{k+n} + c_{k−n} = 0 (7.5.10)
while the nonlinear problem becomes
c_{k+1} − 4c_k + c_{k−1} + c_{k+n} + c_{k−n} − φ^2 Δ^2c^2_k = 0 (7.5.11)
For both the linear and nonlinear case, the boundary conditions are c_i = 1 for all of the nodes. In the linear case, we want to write these equations as
c_k = 1, k ∈ boundary (7.5.12)
since the problem we are solving is of the form
Ac = b (7.5.13)
In contrast, for the nonlinear problem, we are solving a nonlinear algebraic problem of the form
R(c) = 0 (7.5.14)
so we want to write the boundary equations as
c_k − 1 = 0, k ∈ boundary (7.5.15)
We already have enough information to write the program to solve the linear reaction case for a given value of the Thiele modulus. Let’s follow the approach we used in Example 7.4 and write a simple main file:
1 function c = pde_reaction_linear(n,phi2)
2 A = zeros(n^2);
3 b = zeros(n^2,1);
4 A = interior(A,n,phi2);
5 [A,b] = edges(A,b,n);
6 [A,b] = corners(A,b,n);
7 A = sparse(A);
8 c = A\b;
We then have subfunctions for the interior nodes:
1 function A = interior(A,n,phi2)
2 dx = 1/(n-1);
3 for i = 2:n-1
4 for j = 2:n-1
5 k = i + (j-1)*n;
6 A(k,k) = -4 - phi2*dx^2;
7 A(k,k+1) = 1;
8 A(k,k-1) = 1;
9 A(k,k+n) = 1;
10 A(k,k-n) = 1;
11 end
12 end
the edges:
1 function [A,b] = edges(A,b,n)
2 for k = 2:n-1 %bottom
3 A(k,k) = 1;
4 b(k) = 1;
5 end
6 for k=n*(n-1)+2:n^2-1 %top
7 A(k,k) = 1;
8 b(k,1) = 1;
9 end
10 for i = 1:n-2 %left
11 k=n*i+1;
12 A(k,k) = 1;
13 b(k) = 1;
14 end
15 for i = 2:n-1 %right
16 k=i*n;
17 A(k,k) = 1;
18 b(k) = 1;
19 end
and the corners:
1 function [A,b] = corners(A,b,n)
2 A(1,1) = 1; b(1) = 1; %lower left
3 A(n,n) = 1; b(n) = 1; %lower right
4 A(n*(n-1)+1,n*(n-1)+1) = 1; b(n*(n-1)+1) = 1; %upper left
5 A(n^2,n^2) = 1; b(n^2) = 1; %upper right
Notice that our boundary conditions are inhomogeneous, so we have to modify both the vector b and the matrix A in the boundary condition functions. The Thiele modulus appears only in the interior node equations, so we only send it to interior.
Since this is a linear problem, there is nothing interesting that we need to do – just use a direct solver and get the result. The problem can get larger very quickly. For example, if we use Δx = 0.025, we get the nice smooth profile in Fig. 7.14. This is already a reasonably large problem with 41^2 = 1681 unknowns. It is already time to start using a banded solver! However, if we wanted to go to a very large system with an even more fine discretization, it might make sense to switch to an iterative solver and use a continuation method for the initial guess, where we use the converged solution for the previous value of the Thiele modulus as the guess for the next value. For the nonlinear problem, we first write down the residual equations:
1 function R = getR(c,paramet)
2 n = paramet(1);
3 phi2 = paramet(2);
4 R = zeros(n^2,1);
5 dx = 1/(n-1);
6 d = phi2*dx^2;
7 for i = 2:n-1 %interior nodes
8 for j = 2:n-1
9 k = (j-1)*n+i;
10 R(k) = c(k+1)+c(k-1)+c(k+n)+c(k-n)-4*c(k)-d*c(k)^2;
11 end
12 end
13 for k = 1:n
14 R(k) = c(k)-1; %bottom
15 end
16 for k = (n-1)*n+1:n^2
17 R(k) = c(k)-1; %top
18 end
19 for j = 2:n-1
20 k=n*(j-1)+1;
21 R(k) = c(k)-1; %left
22 end
23 for j = 2:n-1
24 k=j*n;
25 R(k) = c(k)-1; %right
26 end
Note that we also could have used the procedure from the linear problem, breaking up the residual calculation into the interior and the boundary nodes. We chose here to keep all of the calculations in a single file. Since each row of the residual is a single equation, the file is still relatively compact.
We also need the entries in the Jacobian. For the boundary nodes, we simply get
J_{k,k} = 1, k ∈ boundary (7.5.16)
For the interior nodes, there are five non-zero entries from the differentiation of Eq. (7.5.11),
J_{k,k−n} = 1 (7.5.17)
J_{k,k−1} = 1 (7.5.18)
J_{k,k} = −4 − 2φ^2 Δ^2c_k (7.5.19)
J_{k,k+1} = 1 (7.5.20)
J_{k,k+n} = 1 (7.5.21)
We can code the Jacobian in a manner similar to the residual:
1 function J = getJ(c,paramet)
2 n = paramet(1);
3 phi2 = paramet(2);
4 dx = 1/(n-1);
5 d=2*phi2*dx^2;
6 J = zeros(n^2);
7 for i = 2:n-1 %interior nodes
8 for j = 2:n-1
9 k = (j-1)*n+i;
10 J(k,k+1) = 1;
11 J(k,k-1) = 1;
12 J(k,k+n) = 1;
13 J(k,k-n) = 1;
14 J(k,k) = -4 -d*c(k);
15 end
16 end
17 for k = 1:n
18 J(k,k) = 1; %bottom
19 end
20 for k = (n-1)*n+1:n^2
21 J(k,k) = 1; %top
22 end
23 for j = 2:n-1
24 k=n*(j-1)+1;
25 J(k,k) = 1; %left
26 end
27 for j = 2:n-1
28 k=j*n;
29 J(k,k) = 1; %right
30 end
As was the case in the residual, we chose to write the Jacobian as a single program rather than breaking it up into interior nodes and boundary nodes. For a given value of the Thiele modulus, we can solve the problem using Newton– Raphson:
1 function [x,flag] = nonlinear_newton_raphson(x,tol,paramet)
2 % x = initial guess for x
3 % tol = criteria for convergence
4 % requires function R = getR(x)
5 % requires function J = getJ(x)
6 % flag = 0 if failed, 1 if converged
7 % n and Da are the paramet
8 R = getR(x,paramet);
9 k = 0; %counter
10 kmax = 100; %max number of iterations allowed
11 while norm(R) > tol
12 J = getJ(x,paramet); %compute Jacobian
13 J = sparse(J);
14 del = -J\R; %Newton Raphson
15 x = x + del; %update x
16 k = k + 1; %update counter
17 R = getR(x,paramet);
18 if k > kmax
19 fprintf('Did not converge.\n')
20 break
21 end
22 end
23 if k > kmax || max(x) == Inf || min(x) == -Inf
24 flag = 0;
25 else
26 flag = 1;
27 end
Since we are using an iterative method, we need to have an initial guess. The easiest approach is to use a zero-order continuation method that starts from the case φ = 0, which is the no-reaction case. Here, we know that the solution is c = 1 everywhere. We can then make a small step in φ, corresponding to a weak reaction, and use the noreaction solution as the initial guess. After converging that solution, we will use it as the initial guess for the next iteration.
If we just want to get a solution at a single value of φ^2, we can use a program like the following:
1 function [x,flag] = pde_continue(phi2_max,n,nsteps)
2 del_phi2 = phi2_max/(nsteps+1);
3 x = ones(n^2,1);
4 phi2 = 0;
5 tol = 1e-8;
6 for i = 1:nsteps+1 %continue up to the desired answer
7 [x,flag] = nonlinear_newton_raphson(x,tol,[n,phi2]);
8 phi2 = phi2 + del_phi2;
9 end
10 [x,flag] = nonlinear_newton_raphson(x,tol,[n,phi2]);
The for loop of this program does zero-order continuation up to a distance Δφ^2 away from the desired value of the square of the Thiele modulus. We then compute the concentration field at the desired value of the square of the Thiele modulus in line 10. While this program will work to compute a single solution for the concentration field, it is not particularly efficient if we want to build up a plot of the concentration versus Thiele modulus since we throw away all of the intermediate calculations.
Before we proceed to answer our original question, it is worthwhile to take a look at the form of the solution for the linear and nonlinear reaction. We can do this by using the programs listed above, along with the program unpack from Example 7.4 to convert the vector solution into a matrix for easy plotting. Figure 7.14 shows the concentration profiles for φ = 1 in both cases. If we did not tell you which profile was the linear or nonlinear reaction, would you be able to figure it out from the plots? This seems rather unlikely.
To determine which type of reaction leads to a lower minimum concentration, Figure 7.15 plots the concentration at the center of the domain as a function of the Thiele modulus. From the symmetry of the problem, the minimum value of the concentration will be at the center of the domain, so this value is a good metric for the extent of reaction.
We already saw how to find the mid-point of the domain in Example 7.4, so we just took advantage of the algorithm in hotspot from that example to output the minimum concentration. For small values of the Thiele modulus, the difference between the linear and nonlinear reactions is negligible. As we increase the rate of reaction, it is clear that the linear reaction leads to the lowest concentration in the center of the domain. The reason for this behavior is that we have selected equal values of the Thiele modulus. Thus, we set the reaction rate in the linear case, k_1, equal to the maximum possible reaction rate k_2c_0 in the nonlinear case. Since the concentration inside the domain is always smaller than that at the boundary, the second-order reaction rate k_2c^2 = k_1c(c/c_0) is always smaller than the first-order reaction rate k_1c. As a result, the concentration is lower in the first-order reaction case.
