WOLFRAM|DEMONSTRATIONS PROJECT

Reaction-Diffusion in a Drop Using Orthogonal Collocation

​
first-order reaction
second-order reaction
Langmuir kinetics 1
Langmuir kinetics 2
concentration profile
error profile
Damköhler number
1
number of interior points
5
Nonlinear reaction diffusion problems expressed mathematically as boundary value problems (ODEs or elliptic PDEs) or initial value problems (parabolic PDEs) are often encountered in chemical engineering applications. A popular numerical method for solving such problems is orthogonal collocation, which reduces the problem to either a system of nonlinear algebraic equations or a system of nonlinear ODEs in time [1, 2]. The accuracy of the method depends on the number of collocation points selected.
In this Demonstration we evaluate the accuracy of the collocation method for studying a family of time-dependent reaction-diffusion problems. We assess the accuracy by comparing the collocation method with the method of lines using NDSolve.
The reaction takes place in a spherical liquid drop. Species
C
is absorbed, then diffuses and undergoes a decomposition reaction in the interior of the drop. We consider the following possible rates of reaction for the decomposition reaction. The four models are:
(1) first-order reaction:
r=-kc
,
(2) second-order reaction:
r=-k
2
c
,
(3) Langmuir-like kinetics 1:
r=-
kc
k
1
c+1
,
(4) Langmuir-like kinetics 2:
r=-
kc
2
(
k
1
c+1)
,
where
k
and
k
1
are appropriate rate constants.
The concentration of species
C
is taken to be a function of
r
only. (In the following IC, BC1, and BC2 stand for initial and boundary conditions.) The number balance for
C
implies:
∂c
∂t
=D
1
2
r
∂
∂r
2
r
∂c
∂r
+r
,
0⩽r⩽R
,
IC:
c(0,r)=0
,
BC1:
∂c
∂r
r=0
=0
,
BC2:
c=
c
eq
at
r=R
,
where
R
is the radius of the drop, taken to be constant during the reaction, and
c
eq
is the concentration of species
C
absorbed at the surface of the drop.
Dimensionless Form
For the analysis we define the following dimensionless variables:
C=
c
c
eq
,
ξ=
r
R
,
τ=
tD
2
R
.
In dimensionless form, the PDE for the first-order reaction kinetics becomes:
∂C
∂τ
=
1
2
ξ
∂
∂ξ
2
ξ
∂C
∂ξ
-
D
a
c
,
0⩽ξ⩽1
,
IC:
C(0,ξ)=0
,
BC1:
∂C
∂ξ
ξ=0
=0
,
BC2:
C=1
at
ξ=1
.
For convenience, we set
k
1
c
eq
=1
so that the problem depends on a single parameter
D
a
=k
n-1
c
eq
2
R
/D
, the Damköhler number (where
n=1
for kinetic models 1, 3, and 4 and
n=2
for model 2). When
D
a
=0
, there is no reaction.
The orthogonal collocation method is based on Jacobi polynomials that reduce the PDE to a system of nonlinear ODEs that are then solved at the collocation points using NDSolve. You can vary the number of interior collocation points and the value of the Damköhler number. This Demonstration plots the dimensionless concentration versus the dimensionless radial position for various values of the dimensionless time using different colors. The colored dots represent the solutions at the collocation points. The colored curves are the predictions found when Mathematica's NDSolve is applied to the nonlinear PDE directly. The blue curve shows the behavior of the concentration at a steady state. This curve was obtained by solving the corresponding boundary value problem using a shooting technique. It is found that a steady state is reached for
τ>0.5
for the range of Damköhler numbers studied.
To quantify the accuracy of the collocation method we compute the Euclidean norm of the error for
p
collocation points given by:
2
ϵ
=
p
∑
i=1
C(τ,
ξ
i
)-
C
e
(τ,
ξ
i
)
2
|
where
C
and
C
e
are the solutions obtained using orthognal collocation and NDSolve. This error is plotted as a function of time on the error plot. The red dots denote the times used in concentration versus radial position plots. As the number of collocation points is increased, the Euclidian norm in the error plot is reduced. For
p=10
and
τ>0.01
,
ϵ∼O(
-4
10
)
.