Page 305 - Numerical Methods for Chemical Engineering

294 6 Boundary value problems

We now apply the finite difference method, to obtain a linear equation for each interior point
(χ i , η j , ζ k ) with the label n = (k − 1)N x N y + (i − 1)N y + j,
θ θ
A n,n−N x N y n−N x N y + A n,n−N y n−N y + A n,n−1 θ n−1 + A nn θ n + A n,n+1 θ n+1
θ θ = σ H(χ i ,η j ,ζ k ) (6.174)
+A n,n+N y n+N y + A n,n+N x N y n+N x N y
Previously, when solving the Poisson equation with Dirichlet boundary conditions, we
obtained a matrix that was positive-definite and could be solved with the conjugate gradient
method. For this problem, however, we have a number of von Neumann boundary condi-
tions, e.g. at the grid points, ( χ i ,η 1 = 0,ζ k ), for which an approximation of the boundary
derivative,
∂θ −θ(χ i ,η 3 ,ζ k ) + 4θ(χ i ,η 2 ,ζ k ) − 3θ(χ i ,η 1 ,ζ k )
= 0 ≈ (6.175)
∂η 2( η)
yields a linear equation for row n of the system,
3θ n − 4θ n+1 + θ n+2 = 0 (6.176)

This row destroys the symmetry of the matrix (and its diagonal dominance as well). There-
fore, we use a method that does not require A to be positive-definite, e.g. bicgstab or gmres.
stove top 3D FD.m performs this calculation for the parameters
σ = 1 a = 0.67 b = 0.5
(6.177)
r 1 /L = 0.1 t 1 /L = 0.05 r 2 /L = 0.25 t 2 /L = 0.05
Options exist for solving the linear system by Gaussian elimination (impractical except
for very small grids), or by the bicgstab and gmres iterative methods. The program also
compares the total heat generation rate within the element to the net heat flux across the
upper surface. As these numbers must agree for the exact solution, comparing their values
provides a measure of accuracy. The results for a grid of 51 × 51 × 25 points are shown in
Figure 6.19 and Figure 6.20.

The MATLAB 1-D parabolic and elliptic solver pdepe

Above, we have focused on solving BVPs by implementing the finite difference method
directly. MATLAB also has a dedicated 1-D BVP solver, pdepe, for systems of equations of
parabolic and elliptic type. Its use is rather straightforward, and for further details type doc
pdepe.

Finite differences in complex geometries

We have used the finite difference method to solve BVPs in rather simple geometries. The
finite difference method can be extended to domains that are composites of these basic
shapes (Figure 6.21) or that can be “stretched” into one of them (Figure 6.22).
For the system in Figure 6.22, we define the transformed coordinates
x y
ξ = η = (6.178)
L h(x)