Previously discussed analytical methods for solving the direct problem in
chemical kinetics are not sufficient for analysis of different reaction kinetic
schemes. First, even given the mathematical model represented by an ODE system,
it is not always possible to integrate the equations analytically. The reason
for that may be just an absence of such a solution. This refers primarily to a
large number of kinetic models in which differential equations are non-linear
relative to the sought functions.
Second, if the analytical solution is obtained, it is often to lengthy
and awkward.
Finally, there is a large class of real mathematical models that are described by partial differential equation sets[1], which cannot be integrated numerically.
Thus, the series of mathematical problems that can be solved
with the previously discussed numerical integration methods is quite narrow.
That is why, in order to solve the direct problem, we have to rely on more
universal approaches. Such approaches are based on
using numerical integration of differential equations and systems.
Among the built-in
tools of the Mathcad suite an important role belongs
to those designed for the numerical solving ordinary differential equations (ODEs)
and their systems. Let’s look through these tools
starting with the solving unit GIVEN/ODESOLVE. The ODESOLVE function first appeared in Mathcad 2000 Pro. On default
this function used Runge-Kutta method of the fourth-level of precision[2].
The organization
of the GIVEN/ODESOLVE solver resembles greatly that of the GIVEN/FIND solving block: it
starts with the GIVEN keyword. An ODE or a system as well as the initial conditions should be placed in the solver body. The solving is performed with a call up of built-in function ODESOLVE using
the following format:
ODESOLVE(x,b,[steps]),
where x is the unknown, and b
is
the upper integration limit. The last parameter, steps, determines
the number of integration steps and is optional. If this parameter is absent,
the number of steps is set up automatically.
As an example, let’s calculate a kinetic curve for an intermediate in a consecutive
reaction where the second step is of the second order:
,
.
Analytical solution of this problem was discussed in the chapter
1.3.2, so we can compare numerically- and analytically-obtained results.
The way of using GIVEN/ODESOLVE solver
and the results of its work are shown in fig. 3.1. The following values were
used: initial concentration A =
0.5M, rate constants =0.05 s-1 and
=0.1 M-1.s-1. Referring to fig. 3.1,
the kinetic curve calculated numerically (dots) matches the curve calculated
using analytical formula (line). The authors want to point out that during the
assigning of the results of symbolic computation one should use only the name
of the desired function (without argument). For the construction of a graph
both function name and its argument are specified fig.
3.1).
Fig. 3.1. Kinetic curve calculation
for an intermediate in a consecutive second-order reaction using ODESOLVE function
(download Mathcad 11 file,
Mathcad 13 file,
Mathcad 14 file; on-line calculation)
The results of this example allow one
to think that the ODESOLVE function is a sufficient tool for solving the direct
kinetic problem. Now we proceed to a numerical solution of an ordinary differential equation (ODE) set
using a solver. Solving process for an ODE set does not differ much from that
for a single ODE: all equations in the set as well as the starting condition
equalities are given in the solver. The only
difference is in the number of arguments that ODESOLVE function should have. An
array of desired function names is required for solving sets of equations. The
variable, where the calculation result is stored, is an array of function names
as well. Below is an example of a numerical calculation of all component concentrations in a
hypothetical multi-step reaction.
.
(3.1)
We can see that the process goes through 6
elementary steps with 5 compounds; and most of the reactions have order higher
than first. Finding the analytical solution of the direct problem is barely
possible in this case. However, having numerical integration tools we do not
need an analytical approach.
The solution of this example is shown in fig. 3.2. One can see that the use of ODESOLVE
function is rather simple.
Fig. 3.2. Calculation of the kinetic
curves for all components in a multi-step reaction using ODESOLVE function
Working with ODESOLVE the user can
choose an algorithm for the numerical ODE solution. If the cursor is above
ODESOLVE keyword, the right-click will activate the corresponding context menu
(fig. 3.3). In Mathcad 2001−13 the user could
choose between integration with a fixed (Fixed) or adaptive (Adaptive) step as
well as a special method for stiff systems integration (Stiff). A new option,
Adams/BDF, was added in Mathcad 14.
Fig. 3.3. Choosing an ODESOLVE algorithm in different
Mathcad versions.
An obvious
advantage of the GIVEN/ODESOLVE solver is that differential equations and ODE
sets are written in the usual manner. That is why
Mathcad documents designed for kinetic problems are rather clear. Another feature of the ODESOLVE
solver is the possibility to use not only differential equations but also usual
algebraic ones.
For example, let’s
consider solving the direct problem for a parallel reaction with the following
mechanism:
,
.
Here reagent A is consumed in two parallel steps
through interaction with reagents B
and C. Each step has the second
order. To find kinetic curves for all the compounds in this reaction it is
enough to use a mathematical model consisting of two differential and three
algebraic equations (the latter represent the material balance of the reaction
system):
,
,
,
,
.
All equations can
be put in the solver without changes. The corresponding numerical
solution is shown in fig. 3.4.
Fig. 3.4. Solving a system of “differential—algebraic” equations
So far, we have considered the
problem of calculating kinetic curves for compounds. In mathematics
this is called a Cauchy problem. It is known that for the Cauchy problem entry
conditions are given, i.e. the value(s) of the desired function y(x) at point x=x0. In
kinetics, as a rule, this point corresponds to the starting point of a
reaction. Often one needs to find a partial solution of a differential equation
using known function values in several points. Such kind of problems is called a boundary-value problem (boundary problem).
The
two-point boundary problem may be
exemplified by the calculation of the kinetic curve for an intermediate
in a consecutive reaction using known values of
concentration at two moments of time (fig. 3.5). Assume the initial A concentration =1 M, rate
constants = 0.008 s-1, k2 = 0.004 s-1, and the concentration of
an intermediate B after 7.2 and 445.1
s from the beginning of the reaction Cb=0.05M. To calculate the kinetic curve for
the compound B we will use the
equation:
,
that
is converted into:
.
The solution of this boundary problem is given
in fig. 3.5 and it does not require additional explanation.
Fig. 3.5. Solving a boundary-value problem with
solver
In addition to Given/Odesolve solver,
one can use built-in Mathcad suite integrators for solving the direct kinetic
problem. These integrators are designed specifically
for the numerical solving of ODEs and ODEs sets. There are several built-in
functions of this kind. It is important to mention that each of the integrators
requires right-hand member array of the ODEs set. Recall that the mathematical
model for a complex reaction is obtained through
multiplying the stoichiometric array by the rate array. The result of this
operation is the array of the ODE set right-hand values. That is why built-in
integrators are very convenient for solving kinetic problems with
stoichiometric arrays of high dimensionality.
We will begin the learning of the
integrators with the rkfixed function. It implements the fourth-order Runge-Kutta method (rk) with
fixed step of integration (fixed). According to Mathcad syntax
this function has five required arguments:
rkfixed(v,x1,x2,npoints,f),
where
v is the intial values vector; x1 and x2 — independent variable values that
determine the interval of integration; npoints — the integration steps
quantity; f — the array function of the ODE set right-hand values.
Creation of the array function f requires
following a specific procedure. Usually this function is
declared right before rkfixed using the following syntax:
f(t,x):=[array of the ODE set right-hand values],
where t is the
independent variable and x is the desired functions array. The right-hand part
of this construction is an array of the ODE set right-hand values; and each of these
functions is represented by an index variable x, i.e. x0,
x1 etc.
The quantity of the array x elements as well as the quantity of the
right-hand values is equal to the quantity of equations in the system. Let’s clarify how rkfixed works by an example of a numerical
solution of the direct problem for a three-step consecutive chemical reaction
(fig. 3.6).
In the beginning, the right-hand
value array Model(k,x) is defined. Current
concentrations of the compounds A, B, C,
D are given as x1, x2, x3, x4 (here 1, 2, 3, 4 are vector indices, not lettered ones!) After assigning the
rate constant values, a function F(t,x)
is declared. The first argument must be an independent variable — time. This
declaration allows the rkfixed integrator to understand the function F(t,x). The numerical calculation result is stored in a
matrix S. Based on the rows and columns quantities in this matrix one may
conclude: the solution matrix has (n+1) columns if the set has n equations. It
is easy to see that the first column starts with a value tn=0 and ends with a value tk=100.
The matrix has (N+1) rows for N integration steps. The second column contains
the value of the variable x2 at each integration step. In this case these
are the values of the reagent A concentration. Correspondingly, the 3rd, 4th and 5th columns
represent concentrations of the reagents B, С and D. Thus, the matrix S give a
pictorial view of the reaction mixture changes over time. Finally, we can make
a graphical representation of the calculation results as kinetic curves for
each reagent.
Thereby, Mathcad tools for numerical
solution of differential equation sets allow one to calculate quickly kinetic
curves of all reactants in a complex chemical process. As the independent
variable in chemical problems is time t,
discussed methods can be used in modeling any time
function (having postulated a differential model for the process). One can find
many analogies of the kinetic models of chemical reactions in other fields of
knowledge (microbiology, sociology etc.) It is useful to discuss the
corresponding examples in order to form some practical skills in creating
differential models.
Let’s consider
one of such problems. Certain microorganisms propagate in proportion to the
colony size (with an aspect ratio k)
but at the same time produce some excrement, which is a poison for them.
Fig 3.6. Numerical
solution of the direct problem for a consecutive reaction with two
intermediates
The rate of
the colony disappearance is proportional to the amount of poison and current
microorganism population with an aspect ratio k1. The poison formation rate is proportional to the number of
microorganisms (with a ratio k2). Suppose
the initial colony size equals N0, and the amount of poison Z is 0 at the beginning. One is required to make the
corresponding set of differential equations and solve it numerically, and
present graphically the microorganisms population trends
along with the amount of poison in the system. Assume k=0,1; k1=0,0001; k2=0,01, N0=2000.
First, we should create a
differential equation set in accordance with the problem specifications.
Changes in the microorganisms population is determined by an increase kN as a result of
reproduction and a decrease –k1NZ due to poisoning. Therefore, the
first differential equation of the system will be of this form:
.
The rate of poison amount change will be described as
.
A differential model for the process has been defined; now we can create a corresponding Mathcad
document to solve the problem (fig. 3.7).
Fig 3.7. Microorganism
population and poison amount trends
As the figure 3.7. implies, the number of microorganisms first
increases with time, achieves the highest value at some point, and after that
the colony becomes extinct. The curve Z(t) is a
typical saturation one. At the beginning, the poison accumulation rate is
small, but it increases with the lapse of time until it reaches the maximum.
Certainly, after full disappearance of the microorganisms the amount of poison
stabilizes and becomes constant. A chemical analog for this model would be a
complex chemical reaction where compounds N
and Z participate in an intermediate
step. In this case, the compound Z is
an autocatalyst for the decomposition of N
in accordance with discussed mathematical correlations.
Further in the
book we will see that a minor (at first sight) modification of the starting
differential equations set and initial conditions may cause a significant
change in the dynamic outlook of the solution.
One may
often have a case when in order to obtain enough accurate result it is necessary
to use a variable integration step: decreasing in the area of large derivative
changes and, vice versa, increasing
when the derivative changes slowly. This approach is
implemented in an algorithm that function Rkadapt uses. During the work of this function the step of integration is adapted in accordance
with the derivative trend in the selected interval of integration. Let’s consider the following ODE system as an example:
,
.
Further, we will discuss minutely this system which is a mathematical model of the widely-known
kinetic scheme bruesselator. Now we will compare the results of its numerical
solution using the functions rkfixed
and Rkadapt.
The corresponding
plots are shown in fig. 3.8. They show that using a
fixed step can lead to an instable solution which can be interpreted
wrongly from a physical point of view (dashed line). The function Rkadapt, as we can see, allows us to eliminate the
mistakes of rkfixed, and reveals
the true behavior of the desired function in the given independent variable
range (due to an adaptable step of integration). In practice, Rkadapt is preferable in the solving of many direct
problems, especially in cases when the starting kinetic model is non-linear.
Fig 3.8. Comparison of the results for calculations with fixed step
of integration (solid line) to an adaptive one (dashed line)
We also
want to mention that the function Rkadapt requires the same five arguments specified in an rkfixed body. Even though integration utilizes a changeable
step, the result will still be represented for evenly
distributed points as specified by user.
There is
one more circumstance related to a variety of built-in integrators. It is the
existence of so-called stiff ODE
sets. The concept of stiffness may be illustrated by the
example of the kinetic equation for a multi-step reaction
.
It is considered that the mathematical model is
stiff if among the eigenvalues λi
of the rate constant matrix there exist such
eigenvalues for which . Usually this condition holds if the rate constants matrix
has elements very different in modulo (three and more orders).
Let’s consider the
following kinetic scheme:
,
that can described with only two differential
equations in accordance with the stoichiometric matrix rank. The corresponding
ODE set can be written as
.
The
eigenvalues vector for this matrix is equal to
,
which allows one to consider the mathematical model
as a stiff one. In this case the integration should be
performed using special built-in functions for stiff systems. Mathcad suite has
integrators Stiffb, Stiffr, Radau, and starting from Mathcad 14 — an integrator AdamsBDF. Fig. 3.9 shows the Mathcad document
illustrating the solution of the direct kinetic problem for a stiff model using
several integrators.
Fig 3.9. An example of a
kinetic scheme described with a stiff set of differential equations (Mathcad
11.2)
As fig. 3.9
shows, the integrator rkfixed can not solve the
problem at all. It shows a diagnostic message «Found a number with a magnitude
greater than 10^307 while trying to evaluate this expression». In this case the integrators Rkadapt, Radau, Stiffb, and Stiffr do work but give different computation
results for the chosen step of integration. One can find that a five-fold
increase in the number of steps virtually levels the results. However, the
integration of more complex stiff sets still should be done
by using the specially designed built-in functions
The
functions Stiffb and Stiffr require an additional argument J(t,x), which is a
matrix of partial time (t,
zero column) and the x
vector derivatives for a kinetic function. A way to form the vector J(t,x) is shown in
fig. 3.9. It is worth mentioning that, in the early versions, this stage was
done manually (paying attention to the fact that the symbol editor cannot differentiate expressions with index
variables). The user had to transform the index variables into lettered ones
and vice versa (see the fig. 3.9
where the index and lettered variables are formatted
with different styles). Because of that the Mathcad 2001i version included the Radau function that did not require the argument J(t,x). Although it
was very convenient, the user had to accept some loss of precision. In the Mathcad 14 version the
functionality of the Radau
was expanded. In addition, this version had tools for automatization of the
matrix J(t,x) symbolic calculation.
The command dsolve of the
Maple system was previously discussed as a method of
analytical solution of the direct problem in chemical kinetics. This command
also can be used for a numerical solving of ODEs or
their sets. In this case one should use the following
syntax:
dsolve({ode,ic},numeric,vars,options)
Here ode is the
differential equation (or ODE set) with the initial conditions ic. The option numeric is a directive for numerical computations (one may use the construction
'type=numeric' instead of the keyword numeric); vars
is the desired function (or the desired function set in case of ODE set); options are additional options given in the keyword=value
form.
The option numeric (or type=numeric) indicates that dsolve will return a numerical calculation result. The
most important additional option is the method. It determines which
numerical method will be used in the desired function
search. The Maple system gives a choice among a variety of methods (see the
list in table. 3.1).
On default, during a solution of the
Cauchy problem the dsolve function automatically utilizes the
Runge-Kutta method of the 4–5th order of accuracy. The major options
are given in the table 3.2.
Table 3.1 – Numerical methods for solving ordinary differential
equations in Maple® suite.
No. |
The ‘method’ value |
Numeric method used
by ‘dsolve’ solver |
1. |
rk45 |
Runge–Kutta–Fehlberg method of the 4–5th order of
accuracy |
2. |
classical or classical[foreuler] |
Euler method |
3. |
classical[heunform] |
Modified Euler method using Heun’s formula |
4. |
classical[impoly] |
Euler method subkind |
5. |
classical[rk2] |
Runge-Kutta method of the 2nd order of accuracy |
6. |
classical[rk3] |
Runge-Kutta method of the 3rd order of accuracy |
7. |
classical[rk4] |
Runge-Kutta method of the 4th order of accuracy |
8. |
classical[adambash] |
|
9. |
classical[abmoulton] |
|
10. |
rosenbrock |
Rosenbrock method |
11. |
bvp |
Numerical method to solve the boundary value problem |
12. |
dverk78 |
Runge–Kutta method of the 7–8th order of
accuracy |
13. |
lsode or lsode[choice], where choice can be adamsfunc, adamsfull, adamsdiag, adamsband, backfunc, backfull, backdiag,
backband |
Modification of the |
14. |
gear, gear[bstoer], gear[polyextr] |
Gear method and its modifications |
15. |
taylorseries |
Method utilizing the |
Table 3.2 – Some major options for the dsolve
command used in numerical calculations
Option |
Purpose |
'output' =
keyword or array |
Controls the output
order of the calculation results. Can
have the symbol values ‘procedurelist’,
‘listprocedure’, as well as ‘array’
or ‘Array’. By default 'output'=
procedurelist. |
'stop_cond' =
list |
Controls the process of computation finishing when the
‘list’ conditions are met. |
'range' = numeric..numeric |
Determines the
independent variable range (the length of the integration interval). |
'stiff'=boolean |
This option is given as a Boolean expression. For example, if 'stiff'=true, the equation is considered stiff and the
default ‘method’ value is set to ‘rosenbrock’ instead of ‘rk45’. |
An important component of using the
dsolve function in numerical calculations is the output order of the results.
In particular, when the following construction is being used:
>dsol:=dsolve({sys1,invalues},numeric);
only
a message about the successful procedure implementation is displayed:
dsol:=proc(rkf45_x)...end
proc
After that, the user has to decide how
the results should be visualized. It is possible to
print the answer in the form of individual values of the desired function, an
array, etc. However, the most visual output form is the graphical one. The
plotting of the results is provided by the command odeplot from
the graphical library plots. Figs. 3.10 and 3.11 show a solution of the differential equation set,
which describes the kinetics of the first-order reversible reaction with arbitrary rate constant values.
Fig 3.10. Numerical solution of the direct kinetic
problem using Mathcad tools
Fig 3.11. Kinetic curves
for reversible reaction participants calculated using numerical calculation
results
In some reactions, one can see the periodic changes of the reagent
concentrations over time. Correspondingly, the rate of the reaction has an
oscillating character. Such reactions are called oscillating
or periodic. Nowadays several dozens of homogeneous and
heterogeneous oscillating reactions have been explored.
Investigations of the kinetic models for these complex processes have allowed formulating
a series of general conditions, which are required for the stable oscillations
of the reaction rates and intermediate concentrations:
q
Stable oscillations appear usually in
open systems, where reagent concentrations can be maintained
constant.
q
The complex reaction must include
autocatalytic steps as well as product-inhibited ones;
q
The reaction mechanism must include
steps with the order higher than first.
These conditions are required but not sufficient for the oscillation to
occur in the system. An important role is played also
by the ratio between the rate constants of certain steps and starting reagent
concentrations. An investigation of the oscillating reactions is still an
important chemical kinetics problem because it is crucial in understanding
catalysis, periodic process laws for living systems, and chemical technology.
Sometimes
chemical problems can be answered using the knowledge
from other sciences that are not related to chemistry at first sight. For
example, some information about a complex reactions flow can
be gained from the mathematical models of the interspecific competition.
A classical example is the “predator–prey” model, which describes the
population trends for predators and prey in living conditions (the Lotka–Volterra model). The
analogies between this model and many reaction schemes are evident.
The main
point of the model is the following. Let’s consider
some closed ecological system that includes two populations: “predators” and herbivorous
“prey”. The population sizes are K
and B correspondingly. The prey
population is reproduced by means of nutrition. The prey consume grass only. The amount of grass is characterized by the T
values. Predators eat solely the prey. Their population trend is determined by the amount of herbivors. There are no
natural enemies for the predators. Instead, the latter experience the natural
loss (diseases, age-specific death). Thereby, the conditions can
be expressed with the following scheme:
,
,
.
Here Вх corresponds
to the number of dead predators. Using chemical terminology, one can say that
the resulting process is a conversion of the starting reagent T into the product Вх. The process is accompanied
by the formation of intermediates K
and B, which serve for different
kinetic functions.
The stages of population expansion are autocatalytic with the
reproduction factors (rate constants) k1
and k2,
but the presence of predator mortality (k3)
prevents the unlimited growth of both populations. Undoubtedly, the overall kinetics
of the process is affected by the amount of grass
necessary for the prey population increase. Assume that we have an unlimited
amount of grass, i.e. T(t)=const. Then, based on the given
conditions, one can write the following differential equation set:
,
.
If the constants for each step and initial
population sizes are given, then the numerical solution will allow one to
predict prey and predator population trends (fig. 3.12).
Fig. 3.12. Population
trends for predators (dashed line) and prey (solid line) in the Lotka–Volterra
model
As fig. 3.12 shows, the population trends for both
populations exibit a continuous wave pattern. For the given initial conditions these oscillations have a constant period and
amplitude. There is an interdependence between the two
population sizes: increase of one of them impedes the growth of other. In the case
of some chemical process descibed by the Lotka-Volterra model, the
concentrations of the intermediates K
and B would be oscillating.
Often it is more convenient to present the solutions of such systems not
as the “concentration over time” trends but rather as the dependence of one of
the concentration on the other. In the second part of fig. 3.12, we have shown the
prey population as a function of the predator population — the so called “phase
path”. The population dynamics can be represented as a
movement along the phase path. The position of the point corresponds to the
population sizes at a given moment of time. One can see that the phase path for
these initial conditions is a closed graph. In the case of continuous oscillations,
the point follows the phase path repeatedly.
Now let’s change the form of the Lotka-Volterra
model by dividing both parts of the equations by k2, and assuming that k2dt=dτ, k/k2=a, k3/k2=b:
,
.
The system has only two parameters now, а and b. Both of them are
positive according to the physical meaning. If one solves the resulting system
for a number of initial conditions, he will end up with a phase path set called
the phase
portrait of the system. For the given case we
will have the phase portrait as a set of closed concentric graphs (fig. 3.13).
Due to the a>0,
b>0 conditions, all the phase
paths are situated in the first quadrant of the coordinate plane. The
oscillations occur near some stationary point, which can be
determined by putting and equal to zero. One can
easily find that ,. The corresponding point (Xst;Yst)
is called the critical point of the system.
Fig. 3.13. Phase portrait of the Lotka—Volterra system with a critical
point
In the case when all the phase paths are closed and
“incapsulated” one into the other, the crytical point is called
a center.
If the initial grass-eater and predator populations
equal b and a correspondingly, the
simulation will not reveal oscillations in the system. Any deviation from the
stationary values will lead to oscillations.
Even though the “predator–prey” model is rather idealized, many kinetic models for real chemical
systems are based on it. For example, D.A.Frank–Kamenetsky
used the Lotka-Volterra model to explain the processes of higher hydrocarbon
oxidation.
The hydrocarbon oxidation kinetics is extremely
complex because it includes many consequent-parallel steps. Thus, the full
mechanism description is problematic. In such cases
researchers confine themselves to model descriptions. Each of the model steps
may represent a series of elementary stages, and each of the model symbols may
correspond to a whole set of compounds playing the same kinetic function.
For example, in the case of a continuous supply of a
gasoline–air mixture into the reactor (heated to certain temperature), one can
see periodical flashes of the cold flame appearing with a constant frequency.
In this case the full combustion does not occur. The oxidation products include
aldehydes, organic peroxides and other compounds. Some regularities
have been established for this process. In particular, the flash frequency
increases with the increase of oxygen concentration and temperature. In order
to explain this effect, D.A.Frank–Kamenetsky suggested the following kinetic
scheme involving two types of intermediates, X and Y:
,
,
Here A is
the starting compound, B is the
product, X is the superoxide type
molecules or radicals, and Y means
the aldehyde type molecules or radicals. One can see that the scheme pustulates
an autocatalysis with the two intermediates. Let’s
assume the reagent concentration does not depend on time (А=const), i.e. its consumption rate is compensated by its insertion
into the reactor. That will give us the following equation set:
,
After the division of both equation parts by k2 we end up with the Lotka-Volterra type
set:
,
,
where a=k1A/k2, b=k3A/k2,
k2dt=dτ. It was shown before that the solution for this system has a
behavior oscillating in time
Let’s show the numerical solution of the Lotka-Volterra
model using Maple suite tools (fig. 3.14). Here the DEplot command from the
DEtools library was used in addition to the dsolve. In
this case, in addition to the integral curves set for the phase paths, the directional
field is visualized. The directional field is
a series of arrows, each of which represents the motion direction along the
phase path at the given point. In many cases the
directional field increases the clearness of the phase portrait. The
directional field can also be plotted using the Maple commands phaseportrait and dfieldplot.
Fig. 3.14. “Predator—prey” model analysis using Maple
To construct the directional field in the Mathcad
environment we recommend the user function field, which was created by T. Gutman (fig. 3.15).
There are other types of critical points. For example,
consider the kinetic scheme consisting of the elementary steps:
,
in which the first step has the zero order. The concentration changes for
the intermediates are described by the equation set:
,
.
Assume , and . We can plot a series of phase portraits for the different
initial concentrations of X and Y based on the numerical solution of the
direct kinetic problem (fig. 3.16).
Fig. 3.15. Phase portrait of the Lotka—Volterra system using a
directional field
The computation results, presented in fig. 3.16, show
that all phase paths converge at one point. Its coordinates are determined by
the values of the intermediates stationary concentrations, and :
,
,
which gives , , or, taking into account the given rate constant values, , . Such type of critical point is called
node,
and oscillations are impossible in this system.
Fig. 3.16. Phase portrait of the system with “node”−type critical
point
Now consider the following equation set:
,
.
This system is also often used
in the differential biological models. If one analyzes the corresponding phase
portrait with the directional field (fig. 3.17), one can see that all the phase
paths approach a critical point and then move away. In this case, we have the “saddle”
critical point[3].
Fig. 3.17. System with a “saddle” critical
point
Finally, there is one more remarkable critical point
type — focus. In order to illustrate it, we will examine one of the
kinetic models of photosynthesis.
In the past there was suggested a mechanism for the
dark steps cycle of photosynthesis. Sugars with different
numbers of carbon atoms, 3 to 7 (trioses, tetroses, pentoses etc.), participare
in this cycle. Having labelled the number with a subscript, one can
create the scheme of the process:
С5
+ С1 + Х → 2 C3 ,
2C3
® C6,
C6
+C3 Û C6,
C4
+ C3 Û 2 C7,
C7 + C3 Û 2 C5.
Here X means
triphosphorpyrydinenucleotide and C1
is carbon dioxide. This kinetic scheme was analysed by D.S.Chernavsky, who assumed some concentrations remaining constant
and ended up with the following differential equation set:
,
.
Let’s solve the set with Mathcad tools (fig. 3.18), using
the numerical integration method with an adaptive step. The constants a0,
a1, a2, b1, b2,
b3 values have been chosen arbitrarily.
The results show that there are periodical concentration oscillations, which
decay over time.
Fig. 3.18. Modelling the photosynthesis kinetics
The phase path is of the spiral form in this case. The
spiral “wraps” around a critical point called focus.
The investigation of the critical point character is closely related to the question of the system stability.
Here the chemical kinetics borrows some terms from the dynamic system theory,
such as Lyapunov’s stability criteria.
Without a deep discussion of the mathematical
apparatus, we will show how the mathematical suites allow one to determine the
critical point type. Assume the mathematical model of a process described by a
set of two differential eqations. In order to find the critical point type one
has to:
q
calculate the critical point coordinates on a phase plane. For
this one has to solve the corresponding algebraic equation set, which is
obtained through the equating of desired functions derivatives to zero;
q
compute
the Jacobian matrix for the system using the critical point coordinates;
q
find the eigenvalues , of the latter matrix.
Obtained eigenvalues allow one to establish the
critical point type and the stability of the stationary state. Six cases are
possible here. They are outlined in fig. 3.19.
λ1,
λ2 – комплексные и Re(λi)=0 Центр |
λ1,
λ2 – действительные и разного знака Седло |
λ1,
λ2 – действительные и отрицательные Устойчивый
узел |
λ1,
λ2 – действительные и положительные Неустойчивый
узел |
λ1,
λ2 – комплексные и Re(λi)<0 Устойчивый
фокус |
λ1,
λ2 – комплексные и Re(λi)>0 Неустойчивый
фокус |
Fig. 3.19. Possible critical point types and phase portraits versus
different Jacobian matrix eigenvalues
As an example we will consider
the previously discussed microorganism propagation model with slight
modifications. Let’s assume that the poison produced
during the microorganism life can decompose (for example, by means of the Sun’s
radiation). The poison decomposition will represent an elementary zero-order
reaction with the rate constant k3.
The new mathematical model will look this way:
.
.
The solution of the direct kinetic problem is shown
in fig. 3.20. We can see that in the case of the assumed rate constants the
microorganism population is oscillating. The critical point type is the node, because all the Jacobian
eigenvalues are imaginary.
Fig. 3.20. Oscillation mode of the population trend in
microorganism colony
The first “chemically grounded” model of a oscillating reaction is a model, which was proposed by I.
Prigozhin and is called bruesselator.
The model is based on a hypothetical reaction with the
following mechanism:
,
,
,
.
It is assumed that the
concentrations of the reagents A and B do not change over time. The
concentrations of D and E are not included in the mass action
law. That is why one need only two equations for the formal
kinetic description of the reaction:
.
It is possible to reduce the number of the controlling
parameters in this system by substituting some variables: , . After
these
changes
the
system
takes the form:
,
,
where , .
A remarkable feature of the bruesselator is the
variety of the critical point types and, consequently, of the phase portraits
depending on the a
and b parameters ratio (table 3.3).
Table
3.3 – Critical point
types for the bruesselator
Parameters a and b ratio |
Critical point type |
b < (a-1)2 |
Stable node |
b < a2+1 |
Stable focus |
b = a2+1 |
Center |
b > a2+1 |
Unstable focus (limit
cycle) |
b > (a+1)2 |
Unstable node |
The case b>a2+1
requires additional examination. The critical point type is the unstable focus.
One can see the appearance of the so-called limit cycle in the phase
portrait. In this case, any point in the phase plane will end up following the same
closed phase path regardless of the initial conditions. This means that stable
asymptotic concentration oscillations (auto-oscillations) of the same amplitude
and frequency will appear with the course of time. Correspondingly, this case
is essentially different from the Lotka-Volterra model. In the latter, one can
find closed phase paths as well, but there is no the only path that does not
depend on the initial conditions. The point set, which
“attracts” all phase paths, was called by I. Prigozhin an attractor. Thus,
the bruesselator has the attractor, while the Lotka-Volterra system does not
have. The appearance of the bruesselator limit cycle can be
seen in fig. 3.21. The Mathcad tools were used to plot the phase
portrait assuming a=1, b=3,25.
Nowadays many real chemical systems are
known, in which processes, accompanied by the concentration oscillations,
take place. These can be both heterogeneous and homogeneous reactions. In
particular, the hydrogen peroxide reduction on the mercury drop surface can
progress periodically in specific conditions. The conjugate process of mercury
surface oxidation is accompanied by a change of the
surface tension. It leads to the drop shape changes. The oscillating mode of
the reaction can be observed through the periodical changes
of the mercury drop shape, which resembles a heartbeat (“mercury heart”).
Fig. 3.21. Brusselator phase
portrait with a limit cycle
The oscillating reactions in homogeneous aqueous media
are of a special interest. Probably, the oxidation of organic acids and their
esters by the bromate ion is investigated the most. B.P. Belousov
(1951) observed the periodic color changes during the oxidation of citric acid
by bromate ion in sulfuric acid solution in the presence of cerium ions. The detailed investigation of this process was done by A.M. Zhabotinsky.
The discovery of this reaction stimulated the investigation of periodical
processes in chemical systems. It became evident that homogeneous oscillating
reactions underlie the most important biochemical processes: generation of
biorhythms and nerve impulses, muscles contraction, etc. As of today, the
reaction of catalytic oxidation of different reducing agents by bromic acid (HBrO3I), following the
auto-oscillating mode, is called the Belousov-Zhabotinsky
reaction. This reaction goes in the acidic water solution and is accompanied by the concentration oscillations for the
oxidized and reduced catalyst forms and intermediates. As a catalyst
one can use transition metal ions, such as manganese or cerium. The reducing
agents can be different organic compounds (malonic acid, acetylacetone, etc.)
We want to mention that, in spite of many publications
dealing with Belousov-Zhabotinsky reaction, the true mechanism of this process
in still unknown. Many kinetic schemes were proposed
to explain the existence of the concentration oscillations. One of the possible
mechanisms is shown in the table 3.4.
Here one can find several important conjugated processes.
1. During the step (1) HBrO2 is formed. It acts as an autocatalyst in the
following reactions.
2.
The extensive chain reaction
of the oxidant BrO3−
with the autocatalyst provides the conditions for the Me+ ions oxidation (steps 4—7).
3.
The oxidation is inhibited due to the chain termination (step 3).
4.
The oxidized form of the
catalyst is reduced during step 14.
The way other reagents react can be
deduced from the given scheme. We have to admit that in spite of the
large number of steps, this kinetic model should be
considered as simplified. However, the solution of the direct kinetic
problem for this scheme at the given conditions (see table 3.4) shows the
presence of stable concentration oscillations. A fragment of the corresponding
Mathcad document is shown in fig. 3.22. Here the
mathematical model was developed in compliance with
the kinetic scheme given in the table 3.4. It was assumed
that the hydrogen ion concentration is constant during the reaction. One can
find the corresponding document on the book’s site.
Table 3.4 – Possible
mechanism of the Belousov-Zhabotinsky reaction
№ |
Step number |
Reaction |
Kinetic parameters values |
1 |
1-2 |
|
k1=2,1; k2=1,0.104 |
2 |
3 |
|
k3=3,0.106 |
3 |
4-5 |
|
k4=4,2 ; k5=4,2.107 |
4 |
6-7 |
|
k6=8,0.104 k7=8,9.103 |
5 |
8 |
|
k8= 3,0.103 |
6 |
9-10 |
|
k9=8,0.109 k10=1,1.102 |
7 |
11 |
|
k11= 4,6.103 |
8 |
12 |
|
k12= 1.106...1.107 |
9 |
13 |
|
k13=1,0.106 |
10 |
14 |
|
k14=2,0.10-1 |
11 |
15 |
|
k15= 3,20.109 |
Fig. 3.22.
Concentration oscillations in
the Belousov–Zhabotinsky reaction
Somewhat
different scheme for the Belousov-Zhabotinsky reaction was suggested by Field,
Korös and Noyes.
The model is called oregonator. It includes
following stages:
,
,
,
,
.
Here A
corresponds to the BrO3−
ion; B corresponds to all organic
reagents that can be oxidized; P is HOBr; X is HBrO2;
Y is the Br− ion; Z
is the reduced form of the catalyst. The mathematical model can
be written as a set of three differential equations:
,
,
.
It is assumed that the concentrations of the
compounds A and B remain constant during the reaction. By using dimensionless variables,
one can transform the set:
,
,
,
where ; ; ; . Now the model
has three controlling parameters:
, , ,
values of which influence greatly the system dynamic behavior.
The document shown in fig. 3.23 can
be used as a template in the computer modeling of the oregonator model.
By changing the controlling parameters, one can see a variety of the reaction
modes with different amplitudes and oscillation frequencies, an appearance of
limit cycles, changes in phase paths trajectories, etc. The value of the stoichiometric
factor f is also of great importance.
Compared to bruesselator, this model is more complex in analysis of possible
stationary states and plotting of the phase portraits. The
reader
can
find
more
details
in
the
specialized
literature.
Fig. 3.23. One
of the direct problem solutions for the oregonator problem
By this point, we were considering only the chemical kinetics problems
in case of constant temperature of the reaction mixture. But
the temperature can change due to ambient conditions (forced heating of
cooling) as well as due to internal factors (heat liberation or adsorption
during reaction). Previously discussed methods are not sufficient in this case.
The mathematical model of the reaction becomes complicated because the temperature
is now a function of time. While equations describing a material balance of the
system are sufficient for isothermic kinetics, in case of altering temperature
one have to consider energy balance as well.
If the reaction has a thermal activation character, changes in
temperature lead to changes the rate constant. This relationship is often described by the Arrhenius equation:
,
where is the activation energy (J/mol),
is the pre-exponential factor. The Arrhenius equation is based on the collision
theory. The theory exploits ideas of an energy barrier and effective collisions
of the reacting particles, which happen in a unit of space over a unit of time.
The k0 value is
proportional to the total collision number. The activation energy determines
the energetic conditions for an active collision — the collision when a transformation
of the reagents into the products is possible. It is a certain excess of energy
in comparison with the average reactants energy, which have to be applied for the reacting species to react. The Arrhenius
equation implies that the reaction rate increases when the temperature rises.
The smaller the activation energy is, the greater such increase will be.
So, the rate
constant depends on temperature. In the case of the altering temperature
the rate constant also becomes a function of time. Consequently, when solving
the direct kinetics problem, we have to add the corresponding equations (the
temperature over time relationships) to the reaction model.
Let’s consider one
of the non-isothermic kinetics cases. Some amount of germanium (IV) chloride is being heated. The heating is
accompanied by a consequent decomposition:
,
.
Assume the heat exchange is organized in a way
that the heating appears with a constant rate 10 K/min. One is asked to establish how the
gross mass of solids will change if 0.002 mol
of GeCl4 is being heated.
The initial temperature is 298 K. It is known that the Arrhenius’ relationships for the rate constants
have the following form:
, .
In these equations the
pre-exponential factors have dimensions of min−1,
activation energies are given in cal/mol.
A mathematical model of the process consists of three differential
equations that describe the changes in reactant amounts:
,
,
as well as of a heating rate
equation:
.
Initial
conditions for the given ODE set are: =0,002; ; =0; .
To integrate this ODE set one can use any of the previously discussed
mathematical suite built-in tools. In fig. 3.24 we
have shown how to solve the problem using the Mathcad built-in function AdamsBDF. The plots allow one to track the
trends of compound concentrations over time. As one can see, during the first
15 mins of heating the amount of the starting material virtually remains
constant. After that, decomposition occurs with a notable rate. When GeCl4 has decomposed, the
solid phase of the reaction mixture consists solely of GeCl2. This composition remains unchanged until
approximately the 30th minute of heating, when the intermediate
begins to decompose into the final product. Finally, some time later the
mixture will consist of pure germanium. One can track the change in mass of the
initial sample in the same way.
When describing the processes of non-isothermic kinetics, it is
convenient to use a unitless variable — conversion of the starting compound X. For example, if we have a single
first-order reaction under the programmed temperature changes conditions, we
can describe its kinetics with a set of two equations:
.
with initial conditions , . The numerical solution of this set for given
kinetic parameters values and linear heating mode (fig. 3.24) shows that a
conversion vs. time plot has a distinct S-shape. A slope of such a curve
changes depending on the given temperature change rate. It important to note
that such curve can be obtained experimentally with a
special device called derivatograph. Information about the system behavior
within a given temperature range is completely enough for solving the inverse
problem, i.e. for the determination of the reaction kinetic parameters using
experimental data.
Fig. 3.24. Solution of the GeCl4
decomposition problem
Fig. 3.25. Conversion
vs. time and temperature for different heating rates
We want to point out that so far we considered only the cases when the
temperature changes were controlled by external factors (i.e. was determined by
the experiment conditions).
If an exothermic reaction takes place in an isolated system, in other
words, when the heat exchange with environment is absent (adiabatic reactor), a
temperature will apparently increase over time. The rate of this increase depends
both on the kinetic parameters (rate constant) and on the thermodynamic
properties of the system (thermal conditions of the reaction, heat capacity).
For a well-mixed periodic reactor, where a single first-order reaction A→B occurs, the mathematical model is described
by this set of equations:
,
,
.
Here ρ is density, kg/m3,
and Cp is specific heat
capacity of the reaction mixture, J/(kg.K).
ΔH is the
reaction heat effect (taking the sing into account), J/mol. To be specific, these parameters depend on temperature. In
addition, heat capacity and density can change as the reaction goes. One should
account for that when performing important calculations. In order not to
overcomplicate, we assume that these values are constant. We will define an
additional parameter: . One can see that
,
An
integration of this equation with the initial conditions T(0)=T0, gives:
.
In the case of adiabatic conditions, the system
will warm up and reach the final temperature that corresponds to an exhaustion
of the reagent:
.
The
temperature Tad is called adiabatic temperature.
The modeling results for the behavior of this system are
shown in fig. 3.26. The initial parameters were:
=1.105 s-1, ET=Ea/R=5000 K, kmol/m3, J=100 K.m3/kmol. The reader can see that performing the calculations is rather
simple.
Fig. 3.26. Temperature and reagent concentration changes in a periodic adiabatic
reactor
In the real conditions, some heat from the
reaction is liberated into environment through the
reactor walls. A differential equation for the temperature changes is given in the following form in this case:
.
In
order to follow operational trends of such non-adiabatic reactor, one have to
introduce additional parameters: heat transfer coefficient h, W/(m2.K), reactor volume V, m3, wall surface S, m2.
In fig. 3.27 one possible way to compute the reagent concentration and reaction
mixture temperature as a function time is shown.
Fig. 3.27. Operation dynamics of a periodic nonadiabatic reactor
For
a well-mixed flow reactor working in an adiabatic mode:
,
,
Here is, and T0
is the reagent temperature when entering the reactor.
The solution of Cauchy problem for this reactor type allows one to conclude:
the dynamic portrait can change strikingly depending on temperature of the
initial mixture. Such situation is illustrated on fig.
3.28. One can see that for the time t=60 s many different kinetic curves as well as temperature—time
relationships are possible, even though the initial temperatures differ only
for 1 K. In both cases a stationary
state is reached. However, for T0 = 274 K
the stationary conversion is low and does not exceed 18%. If the initial
temperature equals 275 К, other
stationary state is reached. The latter corresponds to
rather high conversion (>84%).
Fig. 3.28. Temperature and concentration trends in a flow adiabatic reactor
Thus, the processes taking place in technological reactors can have a multistationarity
even for relatively simple kinetic schemes (in this case
we consider a simple non-reversible first-order reaction). In practice,
reactors work usually in the conditions close to stationary. Therefore, a
problem of optimal organization of the reaction conditions becomes of great
importance. In the discussed example the first stationary state is undesirable from the
efficiency point of view.
For example, assuming A=1.105
s-1, Ea/R=5000 K, J=100 K.m3/kmol, t=60 s, initial concentration CA0 = 1
kmol/m3, and initial temperature T0
= 270 K, we can find three stationary
states. Their quantitative properties are determined
by the solutions of an algebraic equation set, to which the differential
equation set is transformed when both and equal
zero:
,
.
Fig. 3.29. Computations of possible stationary
states and analysis of their stability
The presence of multistationarity can be illustrated
in other way. The last two equations allow one to conclude that
.
A
left part of the obtained correlation depends linearly on the temperature T. Its value is proportional to the
cooling rate caused by a hot (T>T0)
airflow out of the reactor. A right part corresponds to the heat generation
rate in the reactor due to the reaction exothermicity. It is a non-linear
function of temperature. If one plots the temperature dependences of the
equation right and left parts, one will see their interception in points
corresponding to the calculated stationary temperatures and concentrations
(fig. 3.30).
Fig. 3.30. Graphical representation of possible stationary states
When considering systems with many stationary states, it is important to
investigate the stability of the latter. Stability of a stationary state is directly connected to the thermal stability of the
reactor. It may happen that a small perturbation of the system takes it out of
the unstable state. The process will convert into the other one, now stable. In
this case calculations (fig. 3.29) show that two out
of three possible stationary states are stable: for them the Jacobian matrix
eigenvalues are real and of the same sign (stable node). The third stationary
state has real, but negative, Jacobian eigenvalues (saddle point). A comparison
of these results with the plots shown in fig. 3.30 allows one to conclude: a
stationary state is stable if a slope of the heat elimination curve is smaller
than a slope of the heat liberation.
Finally, we can prove the conclusions by plotting a phase portrait of the
system (fig. 3.31). Here dots correspond to the possible stationary states. The
phase paths for different initial temperatures are marked with bold lines.
Fig. 3.31. Phase portrait for exothermic reaction in an adiabatic flow reactor
Fig. 3.31 allows us to see that one or another stationary state is realized according to the initial conditions.
The discussed examples by no means cover all possible problems of chemical kinetics as well as other differential models for chemical-engineering processes, on which the reactor theory is based.
[1] We do not
discuss such mathematical models in this book.
[2] In this book we will not discuss the essence of such and such
numerical integration methods. This information is available in virtually all handbooks
on numerical methods.
[3] In the
document shown in fig. 3.17 (as well as some other documents in this chapter)
there were used user functions IntCurves, VField (T. Gutman). The
reader can find the corresponding documents on the book site.