Simulation of Few Bifurcation Phase Diagrams of Belousov-Zhabotinsky Reaction with Eleven Variable Chaotic Model in CSTR

Simulation of the Gyorgyi, Rempe and Field eleven variable chaotic model in CSTR [Continuously Stirred Tank Reactor] is performed with respect to the concentrations of malonic acid and [Ce(III)]. These simulation studies show steady state, periodic and non-periodic regions. These studies have been presented as two variable bifurcation phase diagrams. We also have observed the bursting phenomenon under different set of constraints. We have given much importance on computer simulation work but not included the experimental methods in this paper.


Introduction
The classical Belousov-Zhabotinsky (BZ) reaction 1,2 is basically a chemical reaction which shows oscillations in the concentrations of few intermediates under certain conditions. It is a prototype example showing non-linear chemical dynamics. It is the simultaneous oxidation and bromination of a suitable organic substrate by acidic bromate in the presence of a metal catalyst. Studies of BZ reaction in well stirred flow reactors have revealed much more than simple periodic oscillations; a rich variety of more complex dynamical behaviors such as complex periodicity, multistability hysteresis and even chaos has been found both in experiments and in simulations. In this reaction, the concentrations of catalyst and / or intermediate species undergo oscillations in time and space. The behavior is driven by the Gibbs free energy decrease of an overall chemical reaction occurring far from equilibrium. A rigorous mechanism for the oscillatory behavior in the BZ reaction was given by Field, Koros and Noyes 3 (FKN) in 1972. Later several authors have reported the behavior experimentally under closed and open system conditions. Schmitz et al 4 have reported these a periodic oscillations experimentally under CSTR conditions in 1977. Since then, several experimental studies have been made at low as well as at high flow rates. For the better understanding of BZ-CSTR chaos, many models came in to existence. An eleven variable model was extracted from the FKN mechanism. Further 9, 7, 4 and 3 variable models were obtained 5 . All models with four or more variables reproduce the major features of the experimentally observed at low flow rate conditions. Hence in our present simulation work, we focus our attention on an 11 variable model proposed by Gyorgyi et al 6 to account both for periodic & non-periodic behaviors. Their study was aimed to model behavior with respect to flow rate only. In this context a detailed simulation of this model has been done to include various input concentrations.

Reaction mechanism and computation method
In a CSTR study, one or more fluid reagents are introduced in to a tank reactor equipped with an impeller while the reactor effluent is removed. The impeller stirs the reagents to ensure proper mixing. Simply dividing the volume of the tank by average volumetric flow rate through the tank gives the residence time since the approach is to report the data in terms of inverse of residence time. Study of the dynamics in CSTR could be made by initially having a batch system in the reactor and this batch system could be converted to the flow system by starting the flows. Another way of studying the dynamics is starting with an empty reactor vessel and studying the dynamics after filling the CSTR. During simulations these practical considerations are not necessary. The volume of the reactor is irrelevant as the flow rate is considered in terms of inverse of residence time. Hence the simulation is done by considering the flow rate, initial concentrations of various species and inflow concentration of major reactants.
The overall BZ reaction is given by 3CH 2 (COOH) 2 + 2BrO 3 -+ 2H + → 2BrCH(COOH) 2 + 4H 2 O + 3CO 2 (1) Table 1 describes the eleven variable GRF chaotic mechanisms consisting of 11 chemical species and 15 reaction steps of which 1, 3, 5 and 6 are reversible. The relevant rate constants are given in Table 2. We have taken this 11-variable model and studied the behaviors of various chemical species like MA & [Ce(III)]. The kinetics of this system in batch can be described by a system of non-linear first order differential equations obtained by applying law of mass action to each step. Such a system of differential equations 7 can be denote by Note: * indicates the radical species.
The initial concentrations of all the variables at time zero are considered as 1X10 - (2) Where, C i is the concentration of i-th species in the reactor at time t, α ij is the stoichiometric coefficient of the i-th species in the reaction j, which proceeds at the rate is calculated from mass action kinetics. The activity of H 2 O is taken as unity. This set of equations gets modified in a CSTR by the flow terms as follows.
. The dynamics of the system could be obtained by integration of this set of differential equations. There are various soft wares are available to solve such problems. We have software called Simulate with us, which incorporates the powerful Rosenbrock's integrator of stiff differential systems. Simulate sets up its own differential equations and result is provided as integrated values, when appropriate mechanism, initial concentrations and input concentrations are given. The error tolerance of integrator used is 10 -7 , initial step size is 10 -8 and maximum step size is 2.0. In the CSTR mode, there are two important initial conditions, namely, the initial concentrations of various species in the reactor and inflow concentration, namely the concentrations of various species which would be obtained in the reactor due to the flows. The initial concentration of all the species in the reactor at time zero is considered as 1X10 -

Results and Discussion
Gyorgyi et al simulated the 11 variable models at the flow rates ranging from 1.00x10 -4 s -1 to 3.25x10 -4 s -1 . They obtained small amplitude oscillations at lower flow rates (1.00x10 -4 s -1 ) and large amplitude oscillations at higher flow rates. (3.25x10 -4 s -1 ) In between these two flow rates they observed period-2, period-4 & even chaos. They also have found bistability in the neighborhood of some flow rates. (1.25-1.50x10 -4 s -1 and 2.70x10 -4 s -1 ) Thus at these flow rates, they found periodic oscillations while the flow is increasing and on the other hand they found nonperiodic oscillations when the flow is decreasing during the reverse path. We have also reproduced their results. In their paper they have reported only the effect of flow rate on chaotic eleven variable models. In addition to their work, we have done simulations with respect to various input concentrations of other chemical species and also we have given the behavior of all the chemical species present in the model. We have chosen a flow rate 3.071x10 -4 s -1 for our work. The eleven variable model presented in the Table 1

Effect of [Ce(III)]
[Ce(III)] 0 was varied in the range 0.0005-0.2 M. Figure 3 represents the traces of [Ce(IV)] at the indicated inflow concentrations and flow rate.

Bursting and Composed oscillations
Bursting and Composed oscillations have been observed by almost all experimentalists working on oscillating reactions in CSTR 8 . In this phenomenon, bursts of large amplitude oscillations are separated by periods of quiescence or small amplitude oscillations, which appear periodically over a long period. As the control parameter is varied the number of oscillations per burst increases or decreases until a bifurcation to the steady state or some other form of oscillatory behavior occurs. In our present study we have observed these bursts at the flow rate of 3.

Bifurcation studies
A bifurcation phase diagram is useful in describing the dynamical behavior of chemical systems. The results of many time series can be shown in single plot whose axes are the two or some times three constraint variables of interest. Each point on the plot indicates what kind of state or states (steady, oscillatory or chaotic) is observed by the lines separating different regions indicating bifurcation where the nature of the solution changes at that particular set of parameter. Two parameter phase plots were constructed by increasing and decreasing one of the parameters under considerations (all other parameters being fixed). The results obtained have been presented as two parameter phase diagrams 9,10 . However, presently we have performed simulations of the system at various (high/low) fixed values of the test constraints without involving the hysteresis process. The observations thus obtained are presented in Figure (