Page 317 - Numerical Methods for Chemical Engineering

306 6 Boundary value problems

2
on a closed domain ⊂ with boundary ∂ . We partition ∂ into two sections, ∂ [d]
[n]
[d]
and ∂ .On ∂ , we apply a Dirichlet boundary condition
ϕ(r) = ϕ B (r) on ∂ [d] (6.204)
[n]
On ∂ , we use a von Neumann “no-flux” boundary condition
∇ϕ · n = 0 on ∂ [n] (6.205)
We generate an unstructured mesh for the domain, and wish to find the nodal field values
that make the residual zero, subject to the imposed boundary conditions. We enforce the
boundary conditions by augmenting (6.202) with two additional integrals involving weight
[n]
[n]
functions w [d] (r) ≥ 0, w (r) ≥ 0, defined respectively on ∂ [d] and ∂ ,
' ' '
[n]

0 = χ p (r)R(r)dr + w [d] ϕ − ϕ B dS + w [∇ϕ · n]dS (6.206)
∂ [d] ∂ [n]
We evaluate the integrals in this equation, starting with the first:
' '
2
χ p (r)R(r)dr = χ p (r)[−∇ ϕ − f (r,ϕ)]dr (6.207)

As the global shape functions interpolate the field using the node values,

N n N n
2 2 2
−∇ ϕ =−∇ ϕ p χ p (r) = ϕ p (−∇ χ p ) (6.208)
p=1 p=1
2
But, since the global shape functions vary only linearly with position, −∇ χ p = 0. There-
fore, even though the exact solution has nonzero second derivatives where f = 0, the global
shapefunctions,andourapproximatesolution,nevercan.Where f = 0,theresidualisnever
exactly zero.
Although we could try using interpolation functions of higher order to allow our approx-
imate solution to have a uniformly zero residual, we instead modify the integrals to remove
any reference to second derivatives. We then obtain what is called a weak solution, in which
we do not enforce that the residual itself actually be zero, but merely that it “integrates out”
to zero locally on the scale of each element.
The integrals (6.207) that we need to compute are
' ' '
2
χ p (r)R(r)dr = χ p (r)[−∇ ϕ]dr − χ p (r) f (r,ϕ)dr (6.209)

2
We rewrite χ p [−∇ ϕ] using the chain rule identity
2
∇ · [χ p ∇ϕ] = (∇χ p ) · (∇ϕ) + χ p (∇ · ∇ϕ) = (∇χ p ) · (∇ϕ) + χ p ∇ ϕ (6.210)
to obtain for the first integral on the left of (6.209)
' ' '
2
χ p [−∇ ϕ]dr = [(∇χ p ) · (∇ϕ)]dr − {∇ · [χ p ∇ϕ]}dr (6.211)