Page 319 - Numerical Methods for Chemical Engineering

308 6 Boundary value problems

We substitute into (6.218) the expansion (6.200) that we now write as

[q]
ϕ(r) = ϕ q χ q (r) + ϕ B r χ q (r) (6.219)
q /∈∂ [d] q∈∂ [d]
[d]
so that for each p /∈ ∂ , (6.218) becomes
  
' '
 

[q]
0 = (∇χ p ) ·  ϕ q ∇χ q + ϕ B r ∇χ q  dr − χ p fdr (6.220)
 q /∈∂ [d] q∈∂ [d] 

Upon rearrangement, this becomes
'

[q]
0 = A pq ϕ q + A pq ϕ B r − χ p (r) f (r,ϕ)dr (6.221)
q /∈∂ [d] q∈∂ [d]

where
'
A pq = [(∇χ p ) · (∇χ q )]dr (6.222)

The elements of A are easy to compute once we know the node positions and mesh topology.
Let [k] denote the region occupied by triangular element k, and let the volume of this
element be V [k] . Then, since the elements are nonoverlapping, we can partition the integral
in (6.222) into contributions from each element,
'

A pq = [(∇χ p ) · (∇χ q )]dr (6.223)
k
[k]
Now, each element integral is zero unless both p and q form a vertex of the element;
therefore, A is a very sparse matrix. For those elements k that do have p and q as vertices,
the gradients of the shape functions are locally constant since we use linear interpolation
(Figure 6.29). Thus, the elements of A are simply

[k]

A pq = ∇χ p [k] · ∇χ q [k] V (6.224)

[p,q]
k∈S E
[p,q]
S E is the set of elements that have both p and q as vertices.
To evaluate the final integral of (6.221), we interpolate the source term from its values
at the nodes,

[q] [q] [q]
f (r,ϕ) = f r ,ϕ q χ q (r) + f r ,ϕ B r χ q (r) (6.225)
q /∈∂ [d] q∈∂ [d]
to obtain
'

[q] [q] [q]
χ p (r) f (r,ϕ)dr = S pq f r ,ϕ q + S pq f r ,ϕ B r (6.226)
q /∈∂ [d] q∈∂ [d]

where we have another sparse matrix S, with elements
' '

S pq = χ p (r)χ q (r)dr = χ p (r)χ q (r)dr (6.227)
[p,q] [k]
k∈S E