Consider the solid bar in Figure 6.7. The temperature of the bar is maintained at the constant value T0 at the boundaries x = −1 and x = 1. In the middle of the bar, from a region [latex]−\epsilon < x < \epsilon [/latex], there is a heating element that delivers heat at a rate q per unit volume. The steady state energy balance for this system is [latex]0 = k \frac{d^2T}{dx^2} + q(x) [/latex] (6.3.38) where k is the thermal conductivity and q(x) is the position-dependent heat generation rate term [latex]q(x)= \begin{cases} 0 &for & −1 ≤ x < −\epsilon \\ q_0 & for &−\epsilon ≤ x ≤ \epsilon \\ 0 & for &\epsilon < x ≤ 1 \end{cases}[/latex] (6.3.39) Compute the temperature profile in this bar as a function of position using finite differences.

Question 6.5: Consider the solid bar in Figure 6.7. The temperature of the...

Numerical Methods with Chemical Engineering Applications

Consider the solid bar in Figure 6.7. The temperature of the bar is maintained at the constant value T0 at the boundaries x = −1 and x = 1. In the middle of the bar, from a region $−\epsilon < x < \epsilon$ , there is a heating element that delivers heat at a rate q per unit volume. The steady state energy balance for this system is

$0 = k \frac{d^2T}{dx^2} + q(x)$ (6.3.38)

where k is the thermal conductivity and q(x) is the position-dependent heat generation rate term

$q(x)= \begin{cases} 0 &for & −1 ≤ x < −\epsilon \\ q_0 & for &−\epsilon ≤ x ≤ \epsilon \\ 0 & for &\epsilon < x ≤ 1 \end{cases}$ (6.3.39)

Compute the temperature profile in this bar as a function of position using finite differences.

We first discretize the spatial domain into n nodes. The corresponding boundary conditions are

$T_1 = T_0 and T_n = T_0$ (6.3.40)

In the interior, we apply centered finite differences to get the discretized form of the ODE

$k\left[\frac{T_{i+1} − 2T_{i} + T_{i−1}}{\Delta x^2} \right]=- q(x_i)$ (6.3.41)

We can rewrite this equation in the form

$T_{i+1} − 2T_{i} + T_{i−1}=Q_i$ (6.3.42)

where the entries in the vector Q have the value

$Q_i=\begin{cases} -\frac{q_0\Delta x^2}{k} &for &i ∈ [i_{min}, i_{max}] \\0& otherwise\end{cases}$ (6.3.43)

The lower bound for the index i where we have non-zero $Q_i$ values is

$-\epsilon =x_l+(i_{min}-1)\Delta x$ (6.3.44)

which reduces to

$i_{min}=\frac{(1-\epsilon)(n-1) }{2}+1$ (6.3.45)

given the definition of Δx and the boundary $x_l = -1$ . Likewise, the maximum index for the heater is determined from

$\epsilon =x_l+(i_{max}-1)\Delta x$ (6.3.46)

which gives us

$i_{max}=\frac{(1+\epsilon)(n-1) }{2}+1$ (6.3.47)

The problem we run into in the computer is that Eq. (6.3.45) generally will not produce an integer value. A simple approximation is to round off $i_{min}$ . The question is whether to round up or down. Consider a small problem with $\epsilon= 1/3$ and n = 5. Equation (6.3.45) gives $i_{min} = 2.33$ . The value of $x_2 = -1/2$ and the value of $x_3 = 0$ . So we should round up,

$i_{min}=\left\lceil \frac{(1-\epsilon)(n-1) }{2}+1\right\rceil$ (6.3.48)
where ⌈·⌉ is the “ceiling” operator that rounds up to the nearest integer. For the same reason that we rounded up for the start of the heater, we now need to round down to get an integer,

$i_{max}=\left\lfloor \frac{(1+\epsilon)(n-1) }{2}+1\right\rfloor$ (6.3.49)

where ⌊·⌋is the “floor” operator that rounds down to the nearest integer.
The original differential equation is linear, so we end up with a linear discretized system of the form

AT = Q (6.3.50)

where A is a tri-diagonal matrix and T is the vector of unknown temperatures. Recall that the effort to solve this equation scales with $np^2$ , where n is equal to the number of nodes and p = 1 is the half bandwidth of the matrix. Recall also that the truncation error for the centered finite differences approximation scales with $Δx^2$ , where Δx is inversely proportional to the number of nodes n. It is therefore clear that there is a trade-off between the computational effort and the error in the solution.

To solve this problem, we can write a MATLAB file to generate the matrix A and forcing vector b:

1 function T = bvp_heater(n,epsilon,q,k)
2 A = zeros(n);
3 b = zeros(n,1);
4 dx = 2/(n-1);
5 imin = ceil((1-epsilon)*(n-1)/2 + 1);
6 imax = floor((1+epsilon)*(n-1)/2+1);

7 Q = -q*dx^2/k;
8 A(1,1) = 1;
9 for i = 2:n-1
10 A(i,i-1) = 1;
11 A(i,i) = -2;
12 A(i,i+1) = 1;
13 end
14 A(n,n) = 1;
15 A = sparse(A);
16 b(1) = 1;
17 for i = imin:imax
18 b(i) = Q;
19 end
20 b(n) = 1;
21 T = A\b;

In this program, lines 8–14 write the coefficient matrix. These terms are independent of the location of the heater, but we do need to be sure to handle the boundaries (Dirichlet conditions) differently than the interior nodes. As was the case in Example 6.2, line 15 declares the matrix to be sparse so that we can take advantage of the tri-diagonal structure. We then use lines 16–20 to write the forcing vector. In line 3, we previously set all the values to be equal to zero, so we just need to write the non-zero entries. Finally, line 21 solves the BVP.

We can see that this program is very flexible, allowing us to handle a wide range of different scenarios just by changing the inputs of $\epsilon and q_0$ . (Changing k is unnecessary since it always appears with $q_0$ in the single dimensionless parameter $q_0 Δx^2/k.$ ) Figure 6.9a shows how the temperature profile depends on the heating rate q, while Fig. 6.9b shows how the temperature profile depends on the size of the heater. Note that we can run into some issues in (b) when we make $\epsilon$ very small or close to unity. In both cases, either the spatial domain with the heater (small $\epsilon$ ) or the one without the heater ( $\epsilon≈ 1$ ) become small. In order to resolve the temperature gradients in these regions, we need to have a sufficient number of nodes. Since we are using evenly spaced nodes, this can increase the problem size. Fortunately, we know how to use a banded solver!