ReactionDiffusion in a Drop Using Orthogonal Collocation
ReactionDiffusion in a Drop Using Orthogonal Collocation
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 timedependent reactiondiffusion 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 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:
C
(1) firstorder reaction: ,
r=kc
(2) secondorder reaction: ,
r=k
2
c
(3) Langmuirlike kinetics 1: ,
r=c+1
kc
k
1
(4) Langmuirlike kinetics 2: ,
r=
kc
2
(c+1)
k
1
where and are appropriate rate constants.
k
k
1
The concentration of species is taken to be a function of only. (In the following IC, BC1, and BC2 stand for initial and boundary conditions.) The number balance for implies:
C
r
C
∂c
∂t
1
2
r
∂
∂r
2
r
∂c
∂r
0⩽r⩽R
IC: ,
c(0,r)=0
BC1: =0,
∂c
∂r
r=0
BC2: at ,
c=
c
eq
r=R
where is the radius of the drop, taken to be constant during the reaction, and is the concentration of species absorbed at the surface of the drop.
R
c
eq
C
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 firstorder reaction kinetics becomes:
∂C
∂τ
1
2
ξ
∂
∂ξ
2
ξ
∂C
∂ξ
D
a
0⩽ξ⩽1
IC: ,
C(0,ξ)=0
BC1: =0,
∂C
∂ξ
ξ=0
BC2: at .
C=1
ξ=1
For convenience, we set =1 so that the problem depends on a single parameter =k/D, the Damköhler number (where for kinetic models 1, 3, and 4 and for model 2). When =0, there is no reaction.
k
1
c
eq
D
a
n1
c
eq
2
R
n=1
n=2
D
a
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 for the range of Damköhler numbers studied.
τ>0.5
To quantify the accuracy of the collocation method we compute the Euclidean norm of the error for collocation points given by: =C(τ,)(τ,) where and 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 and , .
p
2
ϵ
p
∑
i=1
ξ
i
C
e
ξ
i
2

C
C
e
p=10
τ>0.01
ϵ∼O()
4
10