Mathematical Models for a Chemical Reactor

Prévia do material em texto

Mathematical Models for a Chemical Reactor 
 
 
Ignatios Vakalis 
 
Department of Computer Science 
 
California Polytechnic and State University 
San Luis Obispo, CA 93406 
 
(email: ivakalis@calpoly.edu) 
 
(Dated: Sept. 18, 2010) 
 
 
 
 
 
 
Abstract 
 
The module presents the development and solution methodology of mathematical models 
for the concentration of a chemical in a reactor. The steady state case is first examined by 
providing the development and the solution (both analytical and numerical) process of the 
underlying models. The modeling process is extended to include time dependency by 
presenting the formation of the mathematical model and its numeric solution process. 
 
 
 
 
 
 
 
 
 
Keywords: mass balance, steady state, initial/boundary- value problem, parabolic 
PDEs, implicit methods, finite difference, finite element 
 
 
Work supported by the National Science Foundation under CCLI DUE 061825 
Module Overview: 
 
The main focus of the module is on modeling the concentration of a chemical in a 
reactor. The module consists of two major parts: i) steady state; ii) time dependency, 
presenting the development of the relevant models as well as the interplay of their 
analytical and numerical solutions. 
 
The first part begins with the development of the mathematical model that describes the 
concentration of the chemical under steady-steady assumptions. It also presents both the 
analytic and numerical solution of the modeling differential equations. 
 
The second part of the module presents a detailed development and solution methodology 
of the underlying models that capture “time dependency”. The development of the partial 
differential equation (PDE) is presented along with the underlying assumptions. Since an 
analytical solution is not possible, the module gives a detailed presentation and 
comparison of numerical solutions. 
 
 
 
 
Math and Science Level: 
 
The module is suitable for a variety of courses, such as: Mathematical Modeling, Intro 
to Chemical Engineering, Computational Modeling, and Numerical Analysis. 
 
The module should be presented to students at least at the junior level in a science, 
computational science, or mathematics undergraduate curriculum. 
 
 
 
 
Background the Prerequisites: 
 
The minimal mathematical background consists of a Calculus II course. In addition, an 
understanding of differential equations is highly desirable. An intro course in Chemistry 
will provide the necessary science background. An undergraduate course in Numerical 
Analysis is not necessary, but it could add valuable and insightful experience to the 
student reader. 
 
The module requires familiarity with (and access to) a Computer Algebra System (i.e, 
Maple, Mathematica) as well as Matlab 
 
 
 
Modeling Concentration of a Chemical in a Reactor 
 
The law of conservation of mass, states that the mass of a closed system (in the sense of a 
completely isolated system) will remain constant over time. Thus the mass of an isolated 
system cannot be changed as a result of processes acting inside the system. In other 
words, mass cannot be created/destroyed, although it may be rearranged in space, and 
changed into different types of particles. 
 
In its generality, the above can be expressed as: 
 
 (rate of accumulation of a chemical in a system) = 
 (rate of flow into the system) 
 - (rate of flow out of the system) 
 + (rate of generation of a chemical in the system) 
 
 
To simplify the geometry in our modeling process, let’s assume a cylindrical reactor with 
a single entry and exit point. We also assume that the chemical is well-mixed vertically 
and laterally and also that the chemical is subject to first order decay. Chemical 
Engineering literature [Cussler E., 1977] provides us with the following general model: 
 
 
Vcz
x
xc
xx
xcDS
x
xcDSx
t
xcxcFxFc
t
cV γ−⎥⎦
⎤⎢⎣
⎡ Δ∂
∂
∂
∂+∂
∂+∂
∂−⎥⎦
⎤⎢⎣
⎡ Δ∂
∂+−=Δ
Δ )()()()()()( (1) 
 
where, V = volume, F= flow in, c = concentration, S = reactor’s cross sectional area, D 
= dispersion coefficient, γ = first order decay coefficient. (See conceptual question #1) 
Note that the concentration c is a function of space and time, thus c = c(x, t) 
 
 
Allowing and , the above equation can be transformed into: 0→Δx 0→Δt
 
),(),(),(),( 2
2
txc
x
txcW
x
txcD
t
txc γ−∂
∂−∂
∂=∂
∂
 (2) 
 (See problem #1) 
 
where W = F / S 
 
 
The above Eq. (2) is characterized as a Parabolic Partial Differential Equation [Allen, 
1998]. 
 
 
 
Steady State: 
 
 
 Model: 
To derive the mathematical equation that models the steady state of the concentration, we 
need to assume that the concentration is independent of time, and thus only depends on 
spatial variation. 
Thus, we set: 0),( =∂
∂
t
txc which implies that: )(xcc = 
 
Therefore, Eq. (2) is reduced to the following second order Ordinary Differential 
Equation (ODE): 
 
 02
2
=−− c
dx
dxW
dx
cdD γ (3) 
 
 
 
 
Solution Methodology: 
 
In this section we will present two types of solutions. Since Eq. (3) is a 2nd order ODE, 
an analytic solution is feasible and we will pursue to derive it. The module also presents 
a numerical solution. Since a numerical solution can only provide an approximation to 
the exact solution, its accuracy can be tested against the exact solution obtained by the 
analytical method. 
 
Our first step is to explicitly define the conditions for the problem at hand. Let’s assume 
that the cylindrically shaped segment of the reactor has length L , at t = 0 no chemical 
exists, and that the chemical is injected into the reactor’s inflow with constant level . inc
 
Thus, 
dx
dcDSFcFcin
)0()0( −= where, F= flow in, c = concentration, S = reactor’s 
cross sectional area, D = dispersion coefficient. 
 
To obtain a solution of the 2nd order ODE (Eq. 3), an extra condition is imposed that 
reflects an additional assumption. Specifically, 0
),( ==
dt
tLxdc
. In other words, the 
dispersion in the reactor does not affect the exact rate. 
 
 
 
 
Gathering the above observations, and observing that W = F / S , we can now list the 
following 2nd order ODE boundary value problem as: 
 
02
2
=−− c
dx
dcW
dx
cdD γ for 0 < x < L 
 
dx
dc
W
Dccin
)0()0( −= (4) 
 
0)( =
dt
Ldc 
 
 
 
Analytic Solution: 
From your course in ODEs, you recall that the solution takes the form: 
rxexc =)(
 
Substituting in Eq. (4), we obtain the characteristic equation: . 02 =−− γWrDr
The discriminant of the above is: which is always positive. Thus the solution 
of Eq. (4) can be represented as: 
γDW 42 +
 
xrxr eKeKxc 21 21)( += (5) 
where the constants and can be obtained using the boundary conditions. 1K 2K
 
 
Using a Computer Algebra System (i.e, Maple), we can obtain lengthy expressions for 
the coefficients and (See problem # 2) 1K 2K
 
 
Note that the mass fluxes for the steady-state solution can be computed using Fick’s fist 
law. Thus, using Eq. (5), we obtain: 
 
)()( 21 2211
xrxr erKerKD
dx
xdcDJ +−=−= (6) 
 
 
 
 
 
 
Numerical Solution: 
 
Finite Difference Method: 
This section presents an overview of the Finite Difference Method for numerically 
solving the ODE expressed in Eq. (4). 
Using centered finite differences, [Cheney, Kincaid 2003], for the first and second 
derivatives, we can obtain from Eq. (4) the followingexpression: 
 
( ) 02
2 11
2
11 =−Δ
−−Δ
+− −+−+
j
jjjjj c
x
cc
W
x
ccc
D γ (7) 
 
Our goal is to compute the discrete values for all points of the discrete space. To that 
end we can re-arrange the terms of Eq. (7) in order to form a system of linear algebraic 
equations. Rearranging the terms, we obtain for: 
jc
nj ≤≤0 
 
0
2
12
2
1
11 =⎟⎠
⎞⎜⎝
⎛
Δ+−−⎟⎠
⎞⎜⎝
⎛
Δ+
Δ+⎟⎠
⎞⎜⎝
⎛
Δ+− +− jjj cxW
Dc
xW
D
W
xc
xW
D γ
 (8) 
 
The above represents a liner system of (n+1) equations with (n+1) unknowns. 
 
 
In order to solve the above system, we will use the discrete version of the boundary 
conditions to calculate the values for (generated when j=0) as well as 
(generated when j = n). 
1−c 1+jc
Using the boundary condition (Eq. (4)) and discretizing the derivative, we obtain: 
 
 
x
cc
W
Dccin Δ
−−= −
2
11
0 which can be solved to obtain: 
 
011
22 c
D
xWc
D
xWcc in
Δ−Δ+=− (9) 
 
 
Substituting the above expression into the first equations of the linear system (for j=0) 
we obtain: 
 
inc cD
xWc
xW
Dc
xW
D
D
xW
W
x ⎟⎠
⎞⎜⎝
⎛ +Δ=⎟⎠
⎞⎜⎝
⎛
Δ−⎟⎠
⎞⎜⎝
⎛
Δ+
Δ+Δ+ 2222 1γ (10) 
 
 
 
 
Note that at the end of the segment (x = L), the following boundary condition my hold: 
dx
dc
W
Dcc nnn −= This simplifies to: 0=dx
dcn , and its discrete form: 02
11 =Δ
− −+
x
cc nn 
implies that: 11 −+ = nn cc
 
 
Substituting the above expression into the last equation of the linear system ( for j = n) 
we obtain: 
 
022 1 =⎟⎠
⎞⎜⎝
⎛ Δ+Δ+⎟⎠
⎞⎜⎝
⎛
Δ− − nn cW
x
xW
Dc
xW
D γ
 (11) 
 
Note that Eqs. (8, 10, 11) form a tridiagonal system, which could easily be solved either 
by hand, or by using a Computer Algebra System. (See problems #3, #4, #5). 
Upon developing the graph of c(x) as function of the spatial distance x, we observe (as 
expected), that the concentration decreases as the chemical flows though the cylindrical 
reactor. 
Appendix I presents an outline of the Thomas Algorithm for solving efficiently a 
tridiagonal system 
 
Note: 
The Boundary Value Problem (BVP) expressed by Eq. (4) can also be solved using the 
Numerical Solvers within Matlab. Appendix II, provides a sample Matlab code for the 
solution of a BVP. Problem #6 asks the user to solve the BVP using the parameter 
values specified in Problem #3. 
 
 
 
Finite Element Method 
 
Another very popular method for numerically solving PDEs is the Finite Element 
Method (FEM). [ Logan, 1997]. 
 
The Eq. (4) can be re-written as follows, using the substitutions:
d
Wp −= and 
D
q γ−= : 
02
2
=++ qc
dx
dxp
dx
cd
 (12) 
 
and the boundary conditions (see Eq. (4)) can be re-written as: 
0)( =
dx
Ldc
 and ( inccpdx
dc −−= 0)0( ) (13) 
 
In project #1, the reader is asked to develop the complete solution methodology using the 
finite element method. 
 
For any specific set of the problem parameters, graphing the numerical values of the 
concentration c(x) obtained by the Finite Difference and the Finite Element method 
respectively, depicts extremely similar graphs. 
 
 
 
 
Time Dependency: 
 
In the previous sections we have examined the modeling process of the concentration of a 
chemical in a reactor under steady-state. We have observed that the governing parabolic 
PDE (Eq. (2)) reduces to an ODE (Eq. (3)). 
In this part of module, we extend the modeling process by considering time-dependency. 
As we will see, an analytic solution is no longer possible. Thus, we present two numerical 
methods. 
 
 
Model: 
We know that the law of mass concentration states that: 
 
(rate of accumulation of a chemical in a system) = 
 (rate of flow into the system) 
 - (rate of flow out of the system) 
 + (rate of generation of a chemical in the system) 
 
which can be modeled by: 
 
),(),(),(),( 2
2
txc
x
txcW
x
txcD
t
txc γ−∂
∂−∂
∂=∂
∂
 (14) 
 
where F= flow in, c = concentration, S = reactor’s cross sectional area, D = dispersion 
coefficient, γ = first order decay coefficient, and W = F / S. 
 
The above mathematical model (Eq. 14), was formulated under the following 
assumptions: 
a) The reactor tank is well mixed in all directions 
b) Dispersion does not affect the exit rate 
c) Chemical is subject to first order-decay 
d) Prior to t=0 no concentration of the chemical is present in the tank 
e) Starting at t=0 the chemical is injected at a constant rate. 
 
Using the above assumptions, we can formulate boundary conditions for the governing 
model given by Eq. (14) as follows: 
 
 
Assumption (e) can be modeled as: 
),0(),0( ' tc
W
Dctc in += (15) 
Assumption (b) can be modeled as: 
 for t > 0 (16) 0),(' =tLc
 
Assumption (d) can be modeled as: 
 for 0 < x < L (17) 0)0,( =xc
 
 
 
Upon some investigation in the theory and solution methods of Differential Equations, 
it is evident that NO analytical solution exists for the PDE given by Eq. (14) in 
conjunction with the initial and boundary conditions described by Eqs. (15 – 17). Thus, 
the following sections examine two numerical techniques. 
 
 
 
Solution Methodology: 
 
Implicit Method: 
In this part we present an implicit method for solving PDE depicted by Eq. (14). 
 (See Glossary for the difference between explicit and implicit methods) 
 
 
We use a forward finite difference to approximate the time derivative, and central 
differences for the space derivatives, namely: 
 
t
cc
t
c lj
l
j
Δ
−=∂
∂ +1
 which has an error of order )( tO Δ (17a) 
 
 
 
x
cc
x
c lj
l
j
Δ
−=∂
∂ −++
2
1
1
1
 and ( )2
1
1
11
1
2
2 2
x
ccc
x
c lj
l
j
l
j
Δ
+−=∂
∂ +−+++
 (17b) 
 
 
We substitute the above into Eq. (14), and to the initial/boundary conditions in order to 
obtain a system of linear equations, which then must be solved. Specifically, Eq (14) 
becomes: 
 
( ) 1
1
1
1
1
2
1
1
11
1
1
2
2 +
+
−
+
+
+
−
++
+
+
−Δ
−−Δ
+−=Δ
− l
j
l
j
l
j
l
j
l
j
l
j
l
j
l
j c
x
cc
W
x
ccc
D
t
cc γ (18) 
 
 
Rearranging the terms of the above, we obtain: 
 
( ) ( ) ( ) ljljljlj ctcx
W
x
Dc
x
D
t
c
x
W
x
D
Δ−=⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ−Δ+⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ++Δ−⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ+Δ
+
+
++
−
1
2
21
2
1
12
1
2
1
12 γ (19) 
 
 
Observe that the initial condition: c(x,0)= 0 translates in the discrete space as: 
 for j = 0,1, 2, …, n and l=0. 0=ljc
 
 
A complication arises, when we try to explicitly write out the set of equations described 
by Eq. (19). Specifically, for j = 0 we generate the term: 011 =+−lc
 
This term can be evaluated using one of the boundary conditions. Specifically, when we 
discretize Eq. (15) , we obtain: 
 
 ⎟⎠
⎞⎜⎝
⎛
⎥⎥⎦
⎤
⎢⎢⎣
⎡
Δ
−+=
+
−
+
+
W
D
x
cc
cc
ll
j
in
l
2
1
1
1
1
0 
 
Solving the above and substituting in Eq. (19) as applied for j=0, we obtain the first eq. 
of our system that describes the concentration of the nodal points (discrete space): 
 
( ) ( ) linll ctcx
W
D
Wc
D
W
x
W
x
D
t
c
x
D
0
2
1
0
2
2
1
12
122212 Δ−⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ+−=⎟⎟⎠
⎞
⎜⎜⎝
⎛ +Δ+Δ++Δ−⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ
++ γ (20)Now, the boundary condition for t > 0 expressed by Eq. (16), can be 
discretized as follows: 
0),(' =tLc
 
 02
1
1
1
1 =Δ
− +−++
x
cc ln
l
n
 which produces: (21) 
1
1
1
1
+
−
+
+ = lnln cc
 
 
Substituting Eq. (21), into Eq. (19) as applied for j=n, we obtain the last eq. of our 
system that describes the concentration of the nodal points (discrete space): 
( ) ( ) lnlnln ctcx
D
t
c
x
D
Δ−=⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ++Δ−⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ
++
−
1212 12
1
12 γ (22) 
 
 
We are almost at the end….. 
 
Observe that Eqs. (19, 20, 22) form a system of liner equations with unknowns the 
discrete values of the concentration at the nodal points. Furthermore, Eqs. (19,20,22) 
form a tridiagonal system, which could easily be solved either by hand, or by using a 
Computer Algebra System. (See problem #8). 
Upon developing the graph of c(x) as function of the spatial distance x, we observe (as 
expected), that the concentration decreases as the chemical flows though the cylindrical 
reactor. (Appendix I presents an outline of the Thomas Algorithm for solving efficiently 
a tridiagonal system) 
 
Note: 
The modeling problem at hand governed by Eqs. (14-17), can also be solved using the 
Numerical Solvers within Matlab. Appendix IV provides a sample Matlab code for the 
solution of a parabolic PDE. Problem #7 asks the user to modify the code presented in 
Appendix IV, in order to solve Eq. (14) using an implicit method. 
 
 
 
Crank - Nicholson Method 
 
From Eqs. (17a, 17b), we observe that the time difference approximation is first –order 
accurate, and the spatial difference approximation is second-order accurate. 
 
Can we do better? 
 
Indeed yes!. The Crank –Nicholson method provides an improved implicit scheme that 
is second –order accurate in both space and time (i.e, in both independent variables of 
the 2nd order PDE). 
 
For this scheme, the first time derivative can be approximated at point (midpoint) 
by the expression: 
2/1+lt
t
cc
t
c lj
l
j
Δ
−=∂
∂ +1
 
Compare the above with Eq. (17a) that depicts a forward finite divided difference for 
point . 1+lt
 
 
 
To estimate the first and second derivatives in space at midpoint we average the 
difference approximations at points and . Thus: lt 1+lt
 
 
⎥⎥⎦
⎤
⎢⎢⎣
⎡
Δ
−+Δ
−=∂
∂ +−++−+
x
cc
x
cc
x
c lj
l
j
l
j
l
j
222
1 11
1
111 
 
and 
 
( ) ( ) ⎥⎥⎦
⎤
⎢⎢⎣
⎡
Δ
+−+Δ
+−=∂
∂ +−+++−+
2
1
1
11
1
2
11
2
2 22
2
1
x
ccc
x
ccc
x
c lj
l
j
l
j
l
j
l
j
l
j 
 
 
Substituting into Eq. (14) and then rearranging the terms, we can obtain: 
 
 
( ) ( ) ( ) 12112112
1
24242
++
+
+
− ⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ+Δ+−⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ−Δ+⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ+Δ
l
j
l
j
l
j ctx
Dc
x
W
x
Dc
x
W
x
D γ
 = 
 
= ( ) ( ) ( ) ljljlj ctx
Dc
x
W
x
Dc
x
W
x
D
⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ+Δ−−−⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ−Δ−⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ+Δ− +−
1
24242 21212
γ
 (23) 
 
 
 
Again, observe that the initial condition: c(x,0)= 0 translates in the discrete space as: 
 for j = 0,1, 2, …, n and l=0. 0=ljc
 
 
 
A complication arises, when we try to explicitly write out the set of equations described 
by Eq. (23). Specifically, for j = 0 we generate the term: 011 =+−lc
 
This term can be evaluated using one of the boundary conditions. Specifically, when we 
discretize Eq. (15) , we obtain: 
 
 ⎟⎠
⎞⎜⎝
⎛
⎥⎥⎦
⎤
⎢⎢⎣
⎡
Δ
−+=
+
−
+
+
W
D
x
cc
cc
ll
j
in
l
2
1
1
1
1
0 
 
Solving the above, and substituting in Eq. (23) as applied for j=0, we obtain the first 
equation of our system that describes the concentration of the nodal points (discrete 
space): 
 
( ) ( ) 11210
2
2 2
1
2
++ ⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ+⎟⎟⎠
⎞
⎜⎜⎝
⎛ +Δ+Δ+Δ+−
ll c
x
Dc
D
W
x
W
x
D
t
γ
 = 
 
= ( ) ( ) inll cx
DW
D
Wc
x
Dc
D
W
x
W
x
D
t
⎟⎠
⎞⎜⎝
⎛
Δ+⎟⎠
⎞⎜⎝
⎛−⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ−⎟⎟⎠
⎞
⎜⎜⎝
⎛ −Δ−Δ−Δ−−− 222
1
2 12
1
0
2
2
γ
 (24) 
 
 
 
Now, the boundary condition for t > 0 expressed by Eq. (16), can be 
discretized as follows: 
0),(' =tLc
 
 02
1
1
1
1 =Δ
− +−++
x
cc ln
l
n
 which produces: (25) 
1
1
1
1
+
−
+
+ = lnln cc
 
 
 
Substituting Eq. (25), into Eq. (23) as applied for j=n, we obtain the last equation of our 
system that describes the concentration of the nodal points (discrete space): 
 
 
( ) ( ) 12112
1
22
++
− ⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ+Δ+−⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ
l
n
l
n ctx
Dc
x
D γ
 = ( ) ( ) lnln ctx
Dc
x
D
⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ+Δ−−−⎟⎟⎠
⎞
⎜⎜⎝
⎛
Δ− −
1
22 212
γ
 (26) 
 
 
We are almost at the end….. 
 
Observe that Eqs. (23, 24, 26) form a system of liner equations with unknowns the 
discrete values of the concentration at the nodal points. Furthermore, these equations 
form a tridiagonal system, which could easily be solved either by hand, or by using a 
Computer Algebra System. (See problem #9). 
Upon developing the graph of c(x) as function of the spatial distance x, we observe (as 
expected), that the concentration decreases as the chemical flows though the cylindrical 
reactor. 
 
Conceptual Questions: 
 
 
 
 
1) List and explain the necessary assumptions for the development of the 
concentration of a chemical for the steady state case. 
 
 
 
2) Provide a detailed explanation of each term in the formation on Eq. (1). How 
Fick's first law of flux diffusion is embedded in Eq. (1)? (see Glossary) 
 
 
 
3) Investigate and provide reasons for selecting centered finite differences in solving 
Eq. (4) using the finite difference method. Are there other alternatives for 
discretizing the derivatives? How do the various methods rank in terms of 
accuracy? 
 
 
4) List and explain the necessary assumptions for the development of the 
concentration of a chemical for the time-dependency case. 
 
 
 
 
 
Problems: 
 
1) Derive Eq. (2) from Eq. (1). Show all steps. 
 
 
2) Use Maple to obtain the analytic form of the expression for the 1K and 2K in 
the solution Eq. (5). (Hint: Use the boundary conditions – Eq. (4), and solve the 
resulting algebraic system for 1K and 2K ) 
 
 
3) Solve the tridiagonal system described by Eqs. 8, 10,11 to obtain a numerical 
solution of the values: 1,1,2,3.0,10 ===Δ== DWxL γ . Then obtain the 
values of the concentration of the reactor at points x = 0, 2, 4, 6, 8, 10. Plot the 
graph of c(x) as a function of the spatial distance x 
 
4) Compare the values of the concentration obtained in problem #3 with those that 
can be derived from the analytical /exact solution described by Eq. (5). How 
accurate is the numerical solution? 
 
 
5) Solve again the tridiagonal system described by Eqs. 8, 10,11 to obtain a 
numerical solution of the values: 4,1,2,3.0,10 ===Δ== DWxL γ . Plot the 
graph of c(x) as a function of the spatial distance x. Compare the graph with 
that obtained in problem #3. How does the increase of the value of D (depicting 
the turbulence mixing) affects the shape of the graph? Why? 
 
 
6) Solve the BVP described by Eq. (4) using Matlab. Use the parameter values 
given in problem #3. (Hint: See Appendix II) 
 
 
7) Modify the Matlab code given in Appendix IV to numerically solve the PDE with 
the appropriate initial and boundary conditions given by Eqs. (14-17) 
 
 
8) Solve the tridiagonal system described by Eqs. 19, 20, 22 to generate a 
numerical solution of the values: 1,1,2,3.0,10 ===Δ==DWxL γ . Then 
obtain the values of the concentration of the reactor at points: x = 0, 2, 4, 6, 8, 10. 
Plot the graph of c(x) as a function of the spatial distance x. 
 
 
 
9) Solve the tridiagonal system described by Eqs. 23, 24, 26 to generate a 
numerical solution of the values: 1,1,2,3.0,10 ===Δ== DWxL γ . Then 
obtain the values of the concentration of the reactor at points: x = 0, 2, 4, 6, 8, 10. 
Plot the graph of c(x) as a function of the spatial distance x. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Projects: 
 
1) Upon consulting a Numerical Analysis book, develop a complete formulation of 
the Finite Element Method (FEM) in solving the Parabolic PDE in Eqs. 12. 13. 
Graph the results in obtaining the values of c(x) for the same parameter values 
given in problem #3. Compare the results. 
 
 
 
 
 
2) Consider a first order chemical reaction in a tabular flow rector. On the 
assumptions of steady state, laminar flow [Rogers,1992], and negligible axial 
diffusion, the material balance equation is given by: 
 
 0
11 2
2
2
2
0 =−⎟⎟⎠
⎞
⎜⎜⎝
⎛
∂
∂+∂
∂+∂
∂⎥⎦
⎤⎢⎣
⎡ −− kc
r
c
rr
cD
z
c
R
rν 
 
where, v0 = velocity of central stream line, R = radius tube, k = reaction constant, 
D = diffusion constant, c = concentration, z = axial distance, r = radial distance 
from center. 
 
Using the following substitutions, 
 
R
rV
kR
D
c
cCkz ==== ,,, 2
00
βνα 
the original PDE can be transformed to: 
 
 
 ( ) C
V
C
VV
CCV −⎟⎟⎠
⎞
⎜⎜⎝
⎛
∂
∂+∂
∂=∂
∂− 11 2
2
2 βα 
 
 
a) Choosing a set of appropriate boundary conditions, what class of PDE is the 
above equation? ( Hint: See Appendix III) 
b) Set up the equation for numerical solution using the finite difference method. 
c) Will the above step result in an explicit or implicit set of equations? 
d) What are the stability conditions? Explain. 
 
 
 
 
 
References 
 
 
Allen M., et.al., Numerical Analysis for Applied Science, Wiley, New York, (1998). 
 
 
 
Cheney W., Kincaid D., Numerical Mathematics and Computing, Brooks Cole (2003). 
 
 
 
Coombes K., et. al., Differential Equations with Matlab, John Wiley (2000). 
 
 
 
Cussler E., Diffusion: Mass Transfer in Fluid Systems, Cambridge University Press, 
(1997) 
 
 
 
Logan D., A first Course in the Finite Element Method, PWS Publishing Company 
(1997) 
 
 
Rogers, D.F., Laminar flow analysis, Cambridge U. Press (1992) 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Appendix I : 
 
Thomas Algorithm 
 
A tridiagonal system can be represented as: 
 
 A* X = * = 
nn
n
ab
c
ab
cab
ca
00
00
00
1
33
222
11
L
OOOM
O
MO
L
−
⎥⎥
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎢⎢
⎣
⎡
−
n
n
x
x
x
x
1
2
1
M
⎥⎥
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎢⎢
⎣
⎡
−
n
n
d
d
d
d
1
2
1
M
 
Applying an LU decomposition, the above can be written as: 
 
⎥⎥
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎢⎢
⎣
⎡
nn lb
lb
lb
l
00
0
00
0
000
33
22
1
L
OOOM
O
MO
L
* 
⎥⎥
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎢⎢
⎣
⎡
−
1000
0100
10
001
1
2
1
L
OOOO
O
MO
L
nμ
μ
μ
 
 
The solution can be found efficiently by a forward and backward substitution as: 
 
,
1
1
1 l
bw = 
i
iii
i l
wbd
w 1−
−= for i =1, 2,3, …,n 
 
nn wx = and 
i
iii
i l
wwx 1+−= μ for i =1, 2,3, …,n-1 
 
 
 
 
 
 
 
 
 
 
 
Appendix II : 
 
Solving BVPs using Matlab 
 
The following Matlab function twopoint implements the solution of a second order BVP 
given by: 
 
)()()()()()()( 2
2
xFxzxE
dx
xdzxD
dx
xzdxC =++ (a1) 
where z = z(x) 
 
The user must provide the values of the parameters C(x), D(x), E(x) and F(x) at each 
nodal point of the descritized space. The boundary conditions must also be supplied. 
 
 
function y = twopoint(x,C,D,E,F,flag1, flag2, p1,p2) 
% 
% x is a row vector of n+1 nodal points 
% C, D,E, F are vectors describing C(x), D(x), E(x) and F(x) 
% If y is specified at node 1 , flag1 must be set to 1. 
% If y’ is specified at node 1 , flag1 must be set to 0. 
% If y is specified at node n+1 , flag2 must be set to 1. 
% If y’ is specified at node n+1 , flag2 must be set to 0 
% 
n = length(x) – 1; 
 h(2:n+1) = x(2:n+1) – x(1:n); 
 h(1) = h(2); 
h(n+2) = h(n+1); 
r(1:n+1) = h(2:n+2). / h(1:n+1); 
s=1+r; 
if flag1 ==1 
 y(1) =p1; 
 else 
 slope0=p1; 
if flag2 ==1 
 y(n+1) =p2; 
 else 
 slopen=p2; 
end 
W= zeros(n+1, n+1); 
If flag1 ==1 
 c0=3; 
 W(2,2)=E(2) – 2*C(2)/(h(2)^2*r(2)); 
 W(2,3)=2*C(2)/(h(2)^2*r(2)*s(2) + D(2)/h(2)*s(2)); 
 b(2) = F(2) – y(1)*(2*C(2)/(h(2)^2*s(2)) – D(2)/(h(2)* s(2))); 
else 
 c0=2; 
 W(1,1) = E(1) – 2*C(1) / (h(1)^2*r(1)); 
 W(1,2) = 2*C(1)*(1+ 1/r(1))/(h(1)^2*s(1)); 
 b(1) = F(1) + slope0*(2*C(1)/h(1) – D(1)); 
end; 
if flag2 ==1 
 c1=n-1; 
 W(n,n)=E(n) – 2*C(n)/(h(n)^2*r(n)); 
 W(n,n-1)=2*C(n)/(h(n)^2*r(n)*s(n) - D(n)/h(n)*s(n)); 
 b(n) = F(n)– y(n+1)*(2*C(2)/(h(n)^2*s(n)) + D(n)/(h(n)* s(n))); 
else 
 c1=n; 
 W(n+1,n+1) = E(n+1) – 2*C(+1) / (h(n+1)^2*r(n+1)); 
 W(n+1,n) = 2*Cn+(1)*(1+ 1/r(n+1))/(h(n+1)^2*s(n+1)); 
 b(n+1) = F(n+1) + slope0*(2*C(n+1)/h(n+1) + D(n+1)); 
end 
for i= c0:c1 
 W(i,i)= E(i) – 2*C(i) / (h(i)^2*r(i)); 
 W(i, i-1)= 2*C(i)/(h(i)^2**s(i)) - D(i)/h(i)*s(i)); 
 W(i,i+1)= 2*C(i)/(h(i)^2*r(i)*s(i)) + D(i)/h(i)*s(i)); 
 b(i)= F(i); 
end 
 
z = W(flag1 +1:n+1- flag2, flag1 +1 : n+1 – flag2) / b(flag1+1:n+1 – flag2); 
 
if flag1 == 1 & flag2==1, y = [y(1); z; y(n+1)]; end 
if flag1 == 1 & flag2==0, y = [y(1); z]; end 
if flag1 == 0 & flag2==1, y = [z; y(n+1)]; end 
if flag1 == 0 & flag2==0, y = z; end 
 
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
 
 
 
 
 
 
 
 
 
 
 
 
Appendix III : 
 
Classification of 2nd order PDEs 
 
The general form of the 2nd order PDEs with two independent variables x, t can be 
written as: 
 
 
0),(,),(,,,),(),(),(),(),(),( 2
22
2
2
=⎟⎠
⎞⎜⎝
⎛
∂
∂
∂
∂+∂
∂+∂∂
∂+∂
∂
t
txz
x
txzztxf
t
txztxC
tx
txztxB
x
txztxA 
 
For the cases examined in this module, the coefficients A, B, C are constants. 
 
 
 
 
The PDEs are classified into three categories, based on the following criterion: 
 
• If 0 then the PDE is 42 <− ACB elliptic 
 
 
 
• If 0 then the PDE is 42 =− ACB parabolic 
 
 
 
 
• If 0 then the PDE is 42 >− ACB hyperbolic 
 
 
 
 
 
 
 
 
 
 
Appendix IV : 
 
Solving a general Parabolic PDE using Matlab 
 
The following Matlab code provides a numerical solution to the general parabolic PDE: 
 
 22
2 ),(),(
x
txuK
t
txu
∂
∂=∂
∂
 for 0 < x < L 
 
 
function u =differ-equation(nx, hx, nt, ht, init, lowb, hib, K) 
% 
Clear all; 
alpaha = K * ht/hx^2; 
A = zeros( nx – 1, nx – 1); 
u = zeros(nt+1, nx+1); 
u(:, 1)= lowb * ones(nt+1, 1); 
u(:, nx+1)= hib*ones(nt+1, 1); 
u(1,:) = init; 
A(1,1) = 1 + 2*alpha; 
A(1,2) = - alpha; 
for i = 2 : nx – 2 
 A(i, i) = 1 + 2*alpha; 
 A(i, i-1) = - alpha * 
 A( i,i+1) = - alpha; 
end 
A(nx - , nx – 2) = - alpha; 
A(nx – 1, nx – 1) = 1 + 2* alpha; 
b(1,1) = init(2) + init(1) * alpha; 
for i=2: nx – 2 
 b(i,1)= init(i+1); 
end 
b(nx -1, 1) = init(nx) + init(nx+1) *alpha; 
[L, U] = lu(A); 
for j=2: nt+1 
 y = L\b; 
 x = U\Y; 
 u(j,2:nx) = x’; 
 b=x; 
 b(1,1)=b(1,1) + lowb*alpha; 
 b(nx-1,1)= b(nx-1,1) + hib*alpha; 
endPartial Solutions-Hints (For Instructor’s Use ONLY) 
 
 
Problem #1 
 
Recall that: 
x
xcxxc
x
c
Δ
−Δ+=Δ
Δ )()(
, t
tctxc
t
c
Δ
−Δ+=Δ
Δ )()(
 
 
 
 
Also, t
tcttc
t
c
t Δ
−Δ+=∂
∂
→Δ
)()(lim
0
, x
xcxtc
x
c
x Δ
−Δ+=∂
∂
→Δ
)()(lim
0
 and 
 
 
x
x
xc
x
xxc
x
c
x Δ
∂
∂−∂
Δ+∂
=∂
∂
→Δ
)()(
lim
02
2
 
 
 
 
 
 
Project #1 
Develop a complete formulation of the Finite Element Method (FEM) is solving the 
Parabolic PDE in Eqs. 12. 13. 
 
 
Solution: 
Let’s consider a discetization . Assume a tank of length L. Then Nh = L. 
Multiply Eq. (12), with the smooth function v(x) and integrate over the interval [0, l]: 
 
0)( '
0
'' =++∫ dxqcpccL ν (b1) 
 
 
Integrating by parts, we produce: 
0)0()0()()(
00
'
0
'''' =++−− ∫∫∫ LLL dxqcdxpcdxccLLc ννννν (b2) 
 
 
 
From analysis, we know that c(x) can be written in the form: 
 
 where, )()(
0
xaxc
N
j
jj∑
=
= φ )()( xx jφν = and )( ii xcc = at node i 
 
 0 if ji ≠ 
 with =)( ji xφ 
 1, if i =j 
 
 
 
Substituting to Eq. (b2), we obtain: 
0)0()0(
00
'
0
'
0
'
00
' =++−− ∫∑∑ ∫∫∑∑
====
dxajqdxapdxaa j
L
i
N
j
j
N
j
L
o
ijj
L
i
N
j
ji
N
j
jj φφφφφφφφ (b3) 
 
Since, , we have that: )()(
0
xaxc j
N
j
jφ∑
=
= )0()0()(
0
''
0 ∑
=
==−−
N
j
jin cccp φ
 
 
 
Thus, Eq. (b3) can be solved to produce: 
0)0()(
00
'
0
'
0
'
0
0 =++−− ∫∑∑ ∫∫∑
===
dxaqdxapdxaccp j
L
i
N
j
jj
N
j
L
o
ijj
L
i
N
j
jiin φφφφφφφ (b4) 
 
 
Note that: hxx kk /)( 1−−=ϕ if ),( 1 kk xxx −∈ 
 = if hxxk /)( 1 −+ ),( 1+∈ kk xxx 
 
 
Clearly, Eqs. (b4) form a system of linear equations that can be solved to obtain the 
discrete values of the concentration on various nodal points. jc
 
 
 
 
 
 
 
 
Glossary 
 
Deterministic Systems /Models: 
A deterministic system is a system in which no randomness is involved in the 
development of future states of the system that the model describes. Deterministic models 
produce the same output for a given starting condition 
 
Fick's first law: 
It relates the diffusive flux to the concentration field, by stating that the flux goes from 
regions of high concentration to regions of low concentration, with a magnitude that is 
proportional to the concentration gradient (spatial derivative). In one spatial dimension it 
states: 
x
cDJ ∂
∂−= where c= concentration, D = diffusion coefficient 
 
Finite Difference: 
Finite-difference methods are numerical methods for approximating the solutions to 
differential equations using finite difference equations to approximate derivatives. 
 
Laminar Flow: 
Laminar flow, sometimes known as streamline flow, occurs when a fluid flows in parallel 
layers, with no disruption between the layers. It is the opposite of turbulent flow. In 
nonscientific terms laminar flow is "smooth," while turbulent flow is "rough”. 
 
 
Explicit methods 
Explicit methods calculate the state of a system at a later time from the state of the 
system at the current time. Mathematically, if Y(t) is the current system state and 
Y(t + Δt) is the state at the later time (Δt is a small time step), then, for an explicit 
method: 
 
 
 
Implicit Methods 
Mathematically, if Y(t) is the current system state and Y(t + Δt) is the state at the later 
time (Δt is a small time step), then, an implicit method one solves an equation 
 
to find Y(t + Δt). 
It is clear that implicit methods require an extra computation (solving the above 
equation), and they can be much harder to implement. Implicit methods are used because 
many problems arising in real life are stiff, for which the use of an explicit method 
requires impractically small time steps Δt to keep the error in the result bounded .