Question 6.PS.2: Coupled Heat and Mass Transfer Let’s now consider the proble...

In the loop for computing the interior node equations, we have broken up the operators into the linear part and the nonlinear part. We use alpha for the square of the Thiele modulus.

We also need the non-zero entries in the Jacobian. For the first row, we have

J_{1,1} = 1             (6.6.23)

from the left boundary condition on concentration (6.6.19). We get the same result for the second row,

J_{2,2} = 1             (6.6.24)

from the left boundary condition on temperature (6.6.20). For the interior nodes, using the local counter i = 2, 3, ... , n − 1, the non-zero entries in the Jacobian are

J_{2i−1,2i−3} = Δx^{−2}                (6.6.25)

J_{2i−1,2i−1} = −2Δ x^{−2 }− \phi ^2e^{−β/y_{2i}}           (6.6.26)

J_{2i−1,2i} = −\phi ^2 y_{2i−1}e^{−β/y_{2i}} βy^{−2}_{2i}                   (6.6.27)

J_{2i−1,2i+1} = Δx^{−2}             (6.6.28)

for the entries corresponding to the concentration equations (6.6.17) and

J_{2i,2i−2} = Δx^{−2}            (6.6.29)

J_{2i,2i−1} = −γ e^−{β/y_{2i}}           (6.6.30)

J_{2i,2i} = −2Δ x^{−2} − γ y_{2i−1}e^{−β/y_{2i}} βy^{−2}_{2i}            (6.6.31)

J_{2i,2i+2} = Δx^{−2}         (6.6.32)

for the entries corresponding to the temperature equations (6.6.18). For the last two rows, we again get

J_{2n−1,2n−1 }= 1                  (6.6.33)

from the right boundary condition on concentration (6.6.21) and

J_{2n,2n} = 1                     (6.6.34)

from the right boundary condition on temperature (6.6.22). The Jacobian also requires a MATLAB function:

Since the quantity e^{−β/y_{2i}} appears frequently in the Jacobian, we compute it once on line 7 and use it for the rest of the terms at a given interior node. Also, note that we are doing the counting using the local counter i, which simplifies the bookkeeping. Lines 8–11 implement the entries in the Jacobian for concentration, while lines 12–15 do the same for temperature. Lines 4 and 5 are the left boundary condition, and lines 17 and 18 are the right boundary condition.

Since we are passing parameters for the BVP through the program, we also need to modify the Newton–Raphson program (3.5) to take the inputs:

There are a number of parameters that get passed through the program, so for convenience in our main program they are passed as a vector and then unpacked in lines 8–11 to send to the subfunctions.

To solve the problem, we also need to initialize the vector y. A choice y_i = 1 for this problem is normally good enough to get the solution to converge in three iterations of the Newton–Raphson. When we are done with the calculation, we will normally want to unpack the unknowns into two separate vectors for plotting. We can do this with the following program:

This simple program just loops over the unknown vector y and extracts the odd entries into a new vector c and the even entries into a vector T. For convenience, the program also returns the corresponding entries of the independent variable in the vector x.

Parametric Analysis

We now have a set of MATLAB programs that will solve this nonlinear reaction–diffusion problem given an initial condition and a set of parameters. To look at how the concentration and temperature profiles evolve as a function of the parameters, let’s start with the simplest one. The parameter β is a dimensionless activation temperature relative to the temperature at the boundaries. In other words, it sets the bar for the difficulty in achieving a substantial reaction. It is useful then to consider what happens as β becomes large. In this case, there should be no reaction. Indeed, as we see in Fig. 6.16, the temperature profile becomes flatter and flatter as we increase β, while the concentration profile is approaching the linear one that we would expect in the absence of reaction. From a physical standpoint, the dependence on β is not particularly interesting. From a numerical standpoint, however, this serves as an important check on the accuracy of our solution and whether or not we have any bugs in the code. When you are writing more and more complicated programs, it becomes very useful to test them on cases where you know the answer. Indeed, our next test will bring up a problem with our work so far!

If we make the reaction endothermic and pick a relatively small activation temperature, say β = 0.01, we see in Fig. 6.16 that the temperature gets very low without much reaction. This particular example shows a very strong coupling between the transport of heat and species. There is also an analogy with the case of a stiff system of IVPs that we discussed in Section 4.8.2. For stiff IVPs, we have a wide separation of eigenvalues and we need to keep the time step small with respect to the largest eigenvalue so that we have numerical stability. In the case of coupled BVPs, it is possible to have one of the dependent variables change much more rapidly than the other dependent variables. The most rapidly changing dependent variable, in this case the temperature, sets the spatial discretization size Δx that we need to avoid large errors in the solution. Since the equations are coupled, errors in the solution for the temperature will be transmitted to the concentration through the reaction term. The difference between BVPs and IVPs is that there is no equivalent of the stability condition for a BVP – we will obtain a solution to the BVP for large Δx, but it will just be inaccurate. It is unlikely that the solution will diverge, although it is always possible that a nonlinear system will find an unphysical solution due to discretization error.

Let’s consider a more interesting case where we consider the balance between reaction and diffusion embodied in the Thiele modulus. If \phi= 0, then there should be no reaction. Thus, we would expect a small Thiele modulus to be the same as the large β case in Fig. 6.16. This is obviously not the case in Fig. 6.17!

The question we need to answer is whether this is due to a bug in our program, an unphysical solution caused by a bad guess, or some other problem. In general, this is a very difficult question to answer. We normally look at known cases to test out the code before considering other problems where we don’t know the answer. In the case of no reaction, the coupling in Eqs. (6.6.1) and (6.6.3) should disappear – the concentration should be a linear profile between the reservoirs and the temperature should be uniform. The solution we see in Fig. 6.17 gives us a clue to the source of the problem. The concentration profile looks correct, with a perfectly linear solution between the two reservoirs. So the problem must be in the temperature equations. If we look back at the definition of the parameter γ in Eq. (6.6.12), we see that it could be written in the form

γ =\left(\frac{k_0L^2}{D}\right)\left(\frac{Dc_0 H_{rxn}}{VkT_0}\right)           (6.6.35)

where the first term in brackets is the square of the Thiele modulus. Thus, we are not actually free to set the parameters φ and γ independently – if we wanted to set the Thiele modulus to zero and eliminate the reaction, we would have needed the second term to be infinitely large to lead to a non-zero limiting value. This is physically unrealistic.

We put in a deliberate trap in this problem in the definition of γ to bring up an important point about non-dimensionalization of a problem, especially a complicated one that requires numerical methods. This is a powerful approach, but you should be careful that your dimensionless parameters are indeed independent. If you look back to the non-isothermal CSTR problem in Section 5.7, you will see that we were careful to separate out the coefficient for the reactive contribution to the temperature equation into the Thiele modulus and Prater number, which avoids mixing different effects. Fortunately, we were able to use our numerical analysis to find the shortcoming in the problem formulation that we used here.