76
Chapter 12. Application 3 – Chemical Kinetics Simulations 12.1 Overview The simulation of chemical kinetics equations using a computer will be carried out in Assignment 3. In this chapter, the numerical methods for solving the rate equations on a computer will be discussed. In addition, the particular reaction schemes to be studied will be presented. 12.2 Introduction to Ordinary Differential Equations The rate equations of chemical kinetics are first-order ordinary differential equations. A typical ordinary differential equation can be written in the form d y (t ) = f ( y ( t ), t ) , dt
(12.1)
where y ( t ) is the unknown function to be determined from the initial condition y ( t = 0) . The functional form f ( y ( t ), t ) is known € from the differential equation. In most cases in chemical kinetics, there is no explicit time dependence on the right-hand side of Equation (12.1) and the € general form of the equation becomes € € d y (t ) = f ( y ( t )) . (12.2) dt Note that only a single ordinary differential equation is shown in Equations (12.1) and (12.2); however, the results can be generalized for any number of coupled ordinary differential equations. € 12.3 Discretization of Ordinary Differential Equations The computer solution of one set of ordinary differential equations, Newton’s equations, already has been discussed in Chapter 5. In most cases of dealing with differential equations on the computer, the continuous variable must be discretized in order for a solution to be implemented. For chemical kinetics rate equations, the independent variable that will be discretized is the time. The time variable is discretized into equally spaced grid points, t 0 , t1, t 2 , K, t n ,
(12.3)
where the discretized time variable is defined as €
t n +1 = t n + h,
(12.4)
and the time increment between discrete time steps is h. The solution y ( t ) of the ordinary differential equation will be obtained at the discrete time steps, € y 0 , y1, y 2 , K, y n ,
€
€
(12.5)
77
where the notation y n is defined as y n = y (tn ) .
(12.6)
€ 12.4 Solution by Taylor Series Expansion The most common method employed for the numerical solution of ordinary differential € equations involves Taylor Series expansion. The idea is to advance the solution y ( t ) from the present value of t = t n to the future value of t = t n +1. The solution at the future value of time, y n +1 , is expanded in a Taylor Series around the present value of time, t = t n , €
y n +1 =€ y ( t n + h )
€2 h = y n + y ′n ⋅ h + y ′n′ ⋅ + K 2! €
€
(12.7)
The derivative y ′n that appears on the right-hand side of Equation (12.7) is provided by the differential equation, € dy = y ′ = f ( y ( t )) . dt
€
(12.8)
Substituting, Equation (12.7) becomes € h2 y n +1 = y ( t n + h ) = y n + f ( y n ) ⋅ h + y ′n′ ⋅ + K 2!
(12.9)
The truncation of Equation (12.9) after the first two terms on the right leads to the Euler Method, € y n +1 = y n + f ( y n ) ⋅ h .
(12.10)
The Euler Method also is known as the first-order Taylor series method. The solution can be carried out in iterative fashion. €Starting with the known initial value, y 0 = y ( t 0 ) , the function f ( y 0 ) is computed from the differential equation. This allows evaluation of the right-hand side of Equation (12.10), from which the left-hand side, y1 , is determined. This value then can be inserted back into the right-hand side and f ( y1 ) can be computed to yield y 2 . The iterative scheme can be € € continued indefinitely. € 12.5 Solution by Runge-Kutta Methods € provides an iterative solution € The Euler Method to a first-order ordinary differential equation (or set of equations). However, because of the truncation of the Taylor Series expansion after the first derivative, the error in the method is fairly large. In addition, the method scales linearly in the time step. That is, if the time step is decreased by a factor of 10, the error in the numerical solution drops by a factor of 10. Another class of methods that have improved error scaling over the Euler Method includes the Runge-Kutta methods. These methods are based on integration of the differential equation
78
rather than a Taylor Series expansion. The derivation of the Runge-Kutta methods is quite complex and will not be presented here. Instead, the equations for the most popular of the Runge-Kutta methods, the Runge-Kutta fourth order method (RK4), is presented here for ordinary differential equations of the form given by Equation (12.2). y ( t + h) = y (t ) +
1 ( k1 + 2k2 + 2k3 + k4 ) 6
(12.11)
k1 = h ⋅ f ( y ( t )) €
€
€
1 k2 = h ⋅ f y ( t ) + k1 2 1 k3 = h ⋅ f y (t ) + k 2 2 k4 = h ⋅ f ( y (t ) + k 3 )
€ The RK4 equations are much more accurate than the Euler method and therefore they have become one of the most widely used algorithms for solving ordinary differential equations. Rather than € linear error scaling with the stepsize h, the RK4 algorithm scales as h4 . 12.6 Oscillating Chemical Reactions: The Brusselator Model The Brusselator model for an oscillating chemical reaction is described by the following set of steps, k1 A → X k2 B + X → Y + D k
(12.12)
X + X + Y 3 → 3X k
X 4 → E . The rate equations for the intermediates X and Y in the Brusselator model are given by dX = − k1 A − k 2 B X + k3 X 2 Y − k 4 X dt dY = k2 B X − k 3 X 2 Y . dt
(12.13)
The rate equations may be scaled so that dimensionless variables are employed. When this is done, the rate equations have the following form, €
79
dx = a − bx + x 2 y − x dt dy = bx − x 2 y . dt
(12.14)
In Equation (12.14), x represents the scaled concentration of intermediate species X, y represents the scaled concentration of intermediate species Y, and a and b are parameters. € for oscillating chemical reactions is similar to the Lotka-Volterra The Brusselator model model discussed in Chapter 10 in that it has a fixed point. However, the Brusselator model also is more realistic in that it exhibits limit cycle behavior. A limit cycle is a closed curve in phase space to which a set of trajectories is attracted. 12.7 Programming the Brusselator Model To program the Brusselator model for oscillating chemical reactions, the program demonstrated in class for the Lotka-Volterra model may be employed. For the Brusselator model, the lotka.f program must be modified to read in the parameters a and b rather than the rate constants k1 , k2 , and k3 of the Lotka-Volterra model. In addition, the subroutine derives must be modified to evaluate the Brusselator rate equations. € €
Program Input € The input to the Fortran program includes the Brusselator parameters a and b, the time step Δt, the total time of integration, and the initial conditions x (0) and y (0) . These values are read from the screen. The parameters a and b are stored in a COMMON block, so the COMMON statement in the main program and in the subroutines 'derivs' and 'rk4' must be modified accordingly. € € Main Body The main calculation to be carried out by the program is the advancement of the rate equations in time. This is carried out for multiple initial conditions by including a loop to prompt for each desired set of initial conditions to be entered from the screen. The solution is advanced in time using a subroutine that implements the Runge-Kutta 4th order (RK4) algorithm. The RK4 subroutine requires a subroutine called 'derivs' in which the rate equations are defined. Subroutine 'derivs' The subroutine 'derivs' provides the form of the rate equations for the Brusselator model. The variable declaration statement and COMMON block of this subroutine must be modified for the parameters of the Brusselator rate equations. In addition, the rate equations themselvesmust be entered into the main body of the subroutine. Subroutine 'rk4' The subroutine 'rk4' employs the Runge-Kutta 4th order algorithm to solve the rate equations. The only changes that must be made are to the variable declarations and the COMMON block statement. Program Output The output is written to a file. The output includes the time, x ( t ) , and y ( t ) .
€
€
80
12.8 Addendum: Fixed Points and Limit Cycles 12.8.1 Fixed Point Classification The fixed points of a dynamical system may be classified into different categories depending on the motion of trajectories in the vicinity of the fixed point. Two of the most important types of fixed points are elliptic and hyperbolic fixed points, sometimes called stable and unstable fixed points, respectively. The behavior of trajectories in the vicinity of elliptic and hyperbolic points is demonstrated in Figure 12.1.
Figure 12.1 (a) Elliptic fixed point. (b) Hyperbolic fixed point.
The motion of trajectories near an elliptic point is shown in Figure 12.1a. The trajectories orbit around the elliptic fixed point; the motion is stable in the vicinity of the fixed point. We have already seen many examples of elliptic fixed points: in the phase space portraits of the harmonic and Morse oscillators and at the center of the Henon map discussed in Chapter 6. In chemical kinetics systems, the fixed point of the Lotka-Volterra system shown in Figure 10.4 that occurs at X=1.21, Y=0.57 is an elliptic fixed point. The motion of trajectories near a hyperbolic fixed point is shown in Figure 12.2b. In this case, trajectories move away from the vicinity of a hyperbolic fixed point. Examples of hyperbolic fixed points may be found around the periodic islands in the Henon map illustrated in Figure 6.1 as well as in the van der Waals map studied in Assignment 2. Figure 12.2 shows the hyperbolic fixed points around the period 5 islands of the van der Waals map, Equation (8.9), with d=1.4. 0.3
0.2
0.1
Y
0 -0.25
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
0.25
-0.1
-0.2
-0.3
X
Figure 12.2 Hyperbolic fixed points surrounding the period 5 islands of the van der Waals map with d=1.4.
81
Some other types of fixed points are shown in Figure 12.3. These include stable node, unstable node (Figures 12.3a and b, respectively), stable focus, and unstable focus (Figures 12.3c and d, respectively).
(a)
(c)
(b)
(d)
Figure 12.3. Other fixed point types: (a) stable node; (b) unstable node; (c) stable focus; (d) unstable focus.
For a stable node, Figure 12,3a, all the trajectories in the vicinity move toward the fixed point. The opposite is true for an unstable node, Figure 12.3b: all trajectories move away from the fixed point. An example of a stable node was seen in the system of consecutive chemical reactions, Figure 10.2. All the trajectories in the vicinity of a stable focus spiral in towards the fixed point, Figure 12.3c. In contrast, all the trajectories near an unstable focus spiral out away from the fixed point. A stable node or focus is sometimes referred to as an attractor. 12.8.2 Limit Cycles All dynamical systems possess one or more types of fixed points where the motion is stationary. Some dynamical systems also possess other types of dynamical limits called limit cycles. A limit cycle is a fixed curve in phase space that represents the long time periodic behavior of a system. Consider the type of motion that might occur if a trajectory in a dynamical system approaches a stable node. This occurs in the case of the consecutive chemical reactions, Figures 10.1 and 10.2. Motion toward a stable node corresponds to the long-term dynamics coming to rest at one fixed position. A limit cycle is another type of long-term behavior of a system. In this case, however, the system oscillates between two states indefinitely. Consider a trajectory of the Brusselator chemical kinetics system that will be studied in Assignment 3 (see Chapter 13). A
82
graph of the concentration of species x as a function is shown in Figure 12.4. The initial conditions are x=1.0 and y=0.0. 2.5
2
x(t)
1.5
1
0.5
0 0
10
20
30
40
50
60
70
80
90
100
t
Figure 12.4. Plot of x as a function of time for Brusselator model with x0=1.0, y0=0.0.
The system initially oscillates with variable amplitude at short times – this is known as transient behavior. After the transients die out, the system settles into steady oscillations with fixed amplitude. Changing the initial conditions has no effect on the final state; as shown in Figure 12.5, for a different set of initial conditions, only the transients are affected. The final oscillation state is called the limit cycle. 4 3.5 3
x(t)
2.5 2 1.5 1 0.5 0 0
10
20
30
40
50
60
70
80
90
100
t
Figure 12.5. Plot of x as a function of time for Brusselator model with x0=1.0, y0=4.0.
The phase space portrait of the two trajectories of Figures 12.4 and 12.5 for the Brusselator system is shown in Figure 12.6. The closed curve to which both trajectories converge is the limit cycle.
83
4.5 4 3.5 3
Y
2.5 2 1.5 1 0.5 0 0
0.5
1
1.5
2
2.5
3
3.5
4
X
Figure 12.6. Phase space portrait of the Brusselator model showing limit cycle behavior.
