Page 101 - Numerical methods for chemical engineering

Example. 1-D laminar flow 87

and the linear shear-stress profile is

P B
τ yx (y) = y − (2.101)
x 2
Using the Carreau–Yasuda constitutive equation, this yields a nonlinear first-order differ-
ential equation for the velocity field

dv x P B
η(˙γ ) = y − (2.102)
dy x 2
To solve this system, we again use finite differences, and place evenly spaced points between
the plates,
B
y k = k( y) y = (2.103)
(N + 1)
The local velocity and shear-rate at y k are

dv x
v k = v x (y k ) ˙ γ k = (2.104)
dy
y k
At each grid point, we require (2.102) to be satisfied locally,

dv x P B
η(˙γ k ) = y k − (2.105)
dy x 2
y k
The finite difference approximations

dv x v k − v k−1
≈ (2.106)
dy ( y)
y k
yield for each point except the first a nonlinear algebraic equation

v k − v k−1 P B
f k = η(˙γ k ) − y k − = 0 k = 2, 3,. . . , N (2.107)
( y) x 2
where

v k − v k−1
˙ γ k = (2.108)
y
At the first grid point, we use the boundary condition v x (0) = 0 to obtain

v 1 P B v 1
f 1 = η(˙γ 1 ) − y 1 − = 0 ˙ γ 1 = (2.109)
y x 2 y
The Jacobian is sparse as each row has at most two nonzero elements,

∂ f k −1 v k − v k−1 dη d ˙γ k
J k,k−1 = = η(˙γ k ) +

∂v k−1 y ( y) d ˙γ dv k−1 v k−1 ,v k
˙ γ k
(2.110)

∂ f k 1 v k − v k−1 dη d ˙γ k
J k,k = = η(˙γ k ) +
∂v k y ( y) d ˙γ
˙ γ k dv k v k−1 ,v k
We could reduce the effort required to solve the problem by providing the Jacobian as a
second output argument of our function routine. In polymer flow 1D.m, fsolve is provided
with a sparse matrix that is nonzero only at those elements that are nonzero in the Jacobian,
through the JacobPattern option in optimset. This sparsity information helps fsolve to