Continuously Stirred Tank Reactor Using Arc-Length Parameter

​
real part of eigenvalues vs.
D
a
β
2.4
B
20
x
10
0.9
x
20
7
D
am
0.034
Consider a first-order, exothermic, irreversible reaction
A→B
carried out in a well-mixed, continuously stirred tank reactor (CSTR). Assume that the fresh feed of pure
A
is mixed with a recycle stream of unreacted
A
; then the material and energy balances are given by:
V

c
A
t'
=λF
c
A,f
+F(1-λ)
c
A
-F
c
A
-V
k
0
exp
-E
RT
c
A
, (1)
Vρ
C
p
T
t'
=ρ
C
p
F(λ
T
f
+(1-λ)T-T)+V(-ΔH)
k
0
exp
-E
RT
c
A
-hA(T-
T
c
)
, (2)
where
T
and
c
A
are the reactor temperature and the exit concentration of species
A
,
ρ
and
C
p
are the density and heat capacity,
E
and
ΔH
are the activation energy and the heat of reaction,
k
0
and
h
are the reaction rate constant and the heat transfer coefficient,
(1-λ)F
is the recycle flow rate,
c
A,f
and
T
f
are the feed concentration and temperature,
T
c
is the temperature of the cooling medium, and finally
V
is the reactor's volume and
t'
is the time.
Let us define the following nondimensional variables:
t=
t'
τ
(dimensionless time) with
τ=
V
Fλ
,
x
1
=
c
A,f
-
c
A
c
A,f
(dimensionless concentration),
x
2
=
T-
T
f
T
f
(dimensionless temperature),
D
a
=
k
0
-γ
e
V
Fλ
(Damköhler number),
γ=
E
R
T
f
(dimensionless activation energy),
β=
hAτ
Vρ
C
p
(dimensionless heat transfer coefficient),
B=
(-ΔH)
c
a,f
ρ
C
p
T
f
E
R
T
f
(dimensionless adiabatic temperature rise), and finally
x
2c
=
T
c-
T
f
T
f
E
R
T
f
(dimensionless cooling medium temperature).
The governing equations, that is, equations (1) and (2) become:

x
1
t
=-
x
1
+
D
a
(1-
x
1
)exp
x
2
1+
x
2
γ
(3)

x
2
t
=-
x
2
+B
D
a
(1-
x
1
)exp
x
2
1+
x
2
γ
-β(
x
2
-
x
2c
)
(4)
This Demonstration plots the locus of reactor steady states (concentration and temperature) versus the Damköhler number when there is multiplicity for user set values of parameters
B
,
β
, and
x
2c
=7
. The unstable steady states are denoted by the dashed region of the locus. The blue dots denote the values of
D
a
at the turning points, a red dot denotes the value of
D
a
at a Hopf bifurcation, and the green dots denote the
D
a
values where two real eigenvalues change into a complex conjugate pair.
Time series plots for normalized concentration (in blue) and normalized temperature (in magenta) for a given set of initial conditions (
x
10
and
x
20
) and parameters are given. One can investigate the basin of attraction for the steady states. The approach to the steady state may be monotonic or oscillatory depending on the associated eigenvalues for that steady state. Initial conditions near the Hopf point are attracted to a stable periodic solution.

Details

Homotopy method using arc length
Consider a set of
N
nonlinear algebraic equations that depend on a parameter
α
. In vector notation, these equations are represented as
f(u,α)=0
, (5)
where the components of the vectors
f
and
u
are
f={
f
1
,…,
f
N
}
and
u={
u
1
,…,
u
N
}
.
Suppose we have found a solution to equation (5), say by Newton's method for
α=
α
0
. We represent the solution as
u
0
, so that
f(
u
0
,
α
0
)=0
.
We assume that the solution
u
is an analytic function of
α
, that is,
u(α)
is continuous and differentiable. Then, given a solution
u
0
, we can find a neighboring solution
u
on the solution curve by constructing a Taylor series expansion about
u
0
. If the solution curve is a multivalued function of
α
, it is convenient to parameterize the solution in terms of a suitable parameter along the solution curve. The preferred parameter in this study is the arc length
s
along the solution curve so that we can write equation (5) as
f(u(s),α(s))=0
.
The total derivative of
f
along the solution curve is then:
f
s
=
f
u
(u,α)
u
s
+
f
α
(u,α)
α
s
=0
, (6)
which is a set of
N
equations in
N+1
unknowns, where the unknowns are

d
u
1
ds
,
d
u
2
ds
,…,
d
u
N
ds
,
dα
ds

.
Recall that
f
u
(u,α)≡J(u,α)
is the familiar Jacobian used in Newton's method. It is convenient to write equation (6) in matrix notation
[J⋮
f
α
]
u'
α'
=0
, (7)
where
u'≡du/ds
and
α'≡dα/ds
. Thus the vector
t=
u'
α'

is tangent to the solution curve
u(α)
.
The homogeneous system of equations given by equation (7) is an underdetermined system of equations and will have a nontrivial solution if the rank
r
of
[J⋮
f
α
]
has full row rank, namely
r=N
. The nontrivial solution is not unique, however. In order to obtain a unique solution, we need to add an additional equation to equation (7). The appropriate equation is given by the definition of arc length:
u
s
.
u
s
+
2
α
s
=1
. (8)
Equations (7) and (8) represent
N+1
equations in terms of the
N+1
unknowns. Although equation (7) is linear in the unknowns, equation (8) is not. Moreover, the coefficient matrix in equation (7) depends on
u
and
α
, usually in a nonlinear way, which means that we need to solve equations (7) and (8) using an appropriate ODE solver. Matters are complicated further by the fact that
J(u,α)
is singular at bifurcation and turning points. However, the augmented matrix
[J⋮
f
α
]
has full rank at turning points. This feature allows us to devise a numerical algorithm that can integrate through turning points.
Generally one is not really interested in how
u
and
α
depend on
s
, and thus equation (8) is not really needed. The homogeneous system of equations defined by equation (7) can be solved using Cramer's rule. Thus if we define
z={
u
1
,…,
u
N
,α}
, then the vector
z
obeys the system of equations:

z
i
ξ
=
i
(-1)
det[J⋮
f
α
]
-1
,i=1,2,…,N+1
, (9)
where
det[J⋮
f
α
]
-1
is the matrix with the
th
i
​
column removed. In the homotopy literature, Garcia and Zangwill[3] call equation (9) the basic differential equation. The advantage of using it instead of the full system of equations (7) and (8) is that equation (9) is linear in the derivatives. Note that the independent variable
ξ
is not the arc length along the solution curve, but is related to arc length by
s
ξ
=
2

z
1
ξ
+
2

z
2
ξ
+…+
2

z
3
ξ
.
Hence once the solution
z(ξ)
has been found, a simple integration gives the dependence of
s
on
ξ
. The right-hand side of equation (9) can be determined readily when
N
is not too large using Mathematica as shown in the recent paper by Higgins and Binous[2].

References

[1] A. Uppal, W. H. Ray, and A. B. Poore, "On the Dynamic Behavior of Continuous Strirred Tank Reactors," Chemical Engineering Science, 29(4), 1974 pp. 967–985.
[2] B. G. Higgins and H. Binous, "A Simple Method for Tracking Turning Points in Parameter Space," Journal of Chemical Engineering of Japan, 43(12), 2010 pp. 1035–1042.
[3] C. B. Garcia and W. I. Zangwill, Pathways to Solutions, Fixed Points, and Equilibria, Englewood Cliffs, NJ: Prentice–Hall, 1981.

Permanent Citation

Housam Binous, Brian G. Higgins
​
​"Continuously Stirred Tank Reactor Using Arc-Length Parameter"​
​http://demonstrations.wolfram.com/ContinuouslyStirredTankReactorUsingArcLengthParameter/​
​Wolfram Demonstrations Project​
​Published: April 8, 2011