7664
J. Phys. Chem. 1992, 96, 7664-7670
because the alkyl chains stuck into the glass acts to prevent the distortion of the C-C-C bond. However the alkyl chains in n-alkane crystals do not necessarily prevent the distortion of the C-C-C bond. In crystalline n-alkanes all the carbon atoms in a molecule are in the molecular plane. The in-plane accordion motion of a flat carbon skeleton can probably easily take place because the accordion motion does not exert much strain to the surrounding. Since the accordion motion causes the simultaneous bending of all the C-C-C bonds in a molecule, the tunneling C-H hydrogen abstraction can possibly take place on the penultimate and internal carbon atoms. In the previous papers we concluded that the molecular motion of penultimate alkyl radicals in n-alkane crystals is much faster than those of internal radicalsO1'J4At a first glance, this contradicts the present conclusion that the rate of C-C-C bending motion for the internal carbon is comparable to that for the penultimate carbon. However the molecular motion causing paramagnetic relaxation is the out-of-plane twisting motion of the C-C-C bond.14 The in-plane bending motion of the C-C-C bond does not cause the relaxation because the isotropic hyperfine coupling constants of &protons are not modulated by the bending motion. It has been known that the kinds of alkyl radicals and stable products generated by y-irradiation of organic solids depend on the physical phase of the solids.1621For example, y-irradiation of crystalline and glassy 2-methylpropane (isobutane) gives 2methylpropyl and tert-butyl radicals, respectively. However no theory has been proposed for the explanation of the phase effect. We believe that a major part of that effect arises from the effect of phase on the steric hindrance to a particular molecular motion assisting the formation of radicals. The bending motion of the C-C bond is probably restricted in the crystal. However, detailed information on the structure of the crystal is necessary for obtaining a decisive conclusion.
Acknowledgment. This work was supported by a Grant-in-Aid for Developmental Scientific Research No. 04453001 from the Japanese Ministry of Education, Science and Culture. Registry No. CH3(CH2),CH3, 110-54-3; CI(CHZ)JH3, 544-10-5; CH,CH(Cl)(CH,)$H3,638-28-8; CH,CH,CH(CI)(CH2)2CH,, 234681-8; CH3(CH2)5CH3, 142-82-5; CH3(CH2)&H3, 11 1-65-9; D2, 778239-0.
References and Notes (1) Henderson, D. J.; Willard, J. E. J . Am. Chem. Soc. 1969, 91, 3014. (2) Ichikawa, T.; Ohta, N . J . Phys. Chem. 1977,81, 560. (3) Ichikawa, T.; Ohta, N. Radiat. Phys. Chem. 1987, 29, 429. (4) Ichikawa, T.; Yoshida, H. J. Phys. Chem., preceding paper in this issue. (5) Toriyama, K.; Okazaki, M.; Matzuura, K. Radiat. Phys. Chem. 1991, 37, 15, and references cited therein. (6) Iwasaki, M.; Toriyama, K.; Fukaya, M.; Muto, H.; Nunome, K. J . Phys. Chem. 1985,89, 5278. (7) Iwasaki, M.; Toriyama, K.; Muto, H.; Fukaya, M.; Nunome, K. J . Phys. Chem. 1981, 85, 61 1, and references cited therein. (8) Gillbro, T.; Lund, A. Radiat. Phys. Chem. 1976, 8, 625. (9) Claesson, 0.;Lund, A. Chem. Phys. Lett. 1977, 47, 155. (IO) Normant, J. F.;Foulon, J. P.; Deshayes, H. C. R . Seances Acad. Sei., Ser. C 1969, 269, 1325. (11) Ichikawa, T. J . Phys. Chem. 1988, 92, 1431. (12) Muto, H.; Nunome, K.; Fukaya, M. J . Phys. Chem. 1990, 94,418. (13) Tilquin, B.; Baudson, Th. J . Phys. Chem. 1988, 92, 843. (14) Ichikawa, T.; Yoshida, H. J . Phys. Chem. 1988, 92, 5684. (15) Flory, M. Faraday Discuss. 1979, 68, 14. (16) Wakayama, T.; Miyazaki, T.; Fueki, K.; Kuri, Z. J . Phys. Chem. 1970, 74, 3584. (17) Saitake, Y.; Wakayama, T.; Kimura, T.; Miyazaki, T.; Fueki, K.; Kuri, Z . Bull. Chem. Soc. Jpn. 1971, 44, 301. (18) Fukaya, M.; Wakayama, T.; Miyazaki, T.; Saitake, Y.;Kuri, Z. Bull. Chem. Soc. Jpn. 1973, 46, 1036. (19) Miyazaki, T. Radiat. Phys. Chem. 1976,8, 57; 1991, 37,635. (20) Tilquin, B.; Louveaux, M.; Bombaert, C.; Claes, P. Radiat. EfJ 1977, 32._11 _ , ._
(21) Tilquin, B.; Guyot, B.; Miyazaki, T.; Lee, K. P. Radiat. Phys. Chem.
1989, 33, 293.
Gas-Evolution Oscillators. I O . A Model Based on a Delay Equation Kedma Bar-Elit and Richard M. Noyes* Department of Chemistry, University of Oregon, Eugene, Oregon 97403 (Received: February 12, 1992; In Final Form: May 4, 1992)
We have developed a general model for a gas-evolution oscillator with the use of two differential delay equations in two unknowns denoting the population of dissolved molecules in a solution and the pressure in the gas phase. In addition to two measurable rate constants, the only disposable parameter is a delay time which relates to the time before the present at which nucleation occurred. Previous modeling efforts used over 60 ordinary differential equations with as many variables denoting populations of different sizes of bubbles. Our very simple approximateequations model a wide range of experimental situations much better than there would seem to be any reason to expect.
of other organic acids,s of benzene diazonium salts? of hydrogen peroxide,' and of ammonium nitrite.* The behavior of all of these different chemical systems results oscillator occurs when molecules of a volatile species are produced from the same overall process, which takes place in a number of continuously in solution by a homogeneous chemical reaction but parallel or successive component processes: are evolved in repetitive pulses of bubbles. The first example we (a) The reactant chemicals produce molecules of dissolved gas. know of was discovered by Morgan' in 1916 while studying the (b) The dissolved molecules may escape from the solution dehydration of formic acid in concentrated sulfuric acid. Over through its surface. 60 years later, study of the same system was resumed inde(c) If the solution becomes sufficiently supersaturated, small pendently by Showalter and Noyes2 and by Bowers and R a ~ j i . ~ nuclear bubbles will form homogeneously and spontaneously. Smith et a1.4 subsequently showed that the origin of the oscillations (d) The nuclear bubbles will grow by interaction with the was not chemical but was associated with the nucleation, growth, surrounding solution and will rise until they are destroyed at the and escape of bubbles. Oscillatory evolution of gas bubbles due surface. to such processes was also observed during decompositionreactions (e) The design of the apparatus may or may not be such that the gas in contact with solution is separated from the external 'Permanent address: Sackler Faculty of Exact Sciences, School of atmosphere by a physical leak. If such a leak does exist, its Chemistry, Tel-Aviv University, Israel.
I. Introduction
A. CbPRcteristicsof Ges-Evdution Oscilletors. A gas-evolution
0022-36S4/92/2096-1664%03.00/0 0 1992 American Chemical Society
The Journal of Physical Chemistry, Vol. 96, No. 19, 1992 7665
Gas-Evolution Oscillators magnitude may need to be considered during modeling of the system. If the solution is stirred sufficiently vigorously, supersaturation may be kept so low that component processes c and d do not take place. If stirring is less vigorous, process b and processes c d provide alternative paths to accomplish the same overall change of state. If stirring is sufficiently slow, process b can be neglected. B. Reviorro Modeling of Cas-Evolution Oscillntors. The system summarized qualitatively above will be quite complex to model quantitatively. The solution of interest contains dissolved gas molecules and also dispersions of bubbles of many different sizes. According to previous estimates," the radii of those bubbles range from about lod to 0.1 cm,while the pressures in different bubbles range by a factor of the order of 100. An approach similar to that used to model several purely chemical systems would be to divide all bubbles into classes according to size and to consider each class as a separate intermediate species. Such a treatment was first attempted by Smith et al.9J0and was developed more quantitatively by Yuan, Ruoff, and Noyes'' (YRN) to simulate the experimental observations of Kaushik et a1.I2 of the original Morgan' reaction. The model, which was called the bubbelator, described the changing composition of the solution with equations of the form of (1) and (2). dC/dt =
+
J = R , exp[-a/(C - C,,,)Z] In these equations, C is the concentration in molecules cm-3 of homogeneously dissolved gas which is being produced according to process a at a constant rate given by Rch. Process b is described by the term k,, where A and V denote the surface area and volume, respectively, of the total solution. C,,, is the concentration which would be in equilibrium with the gas above the surface. It is a consequence of the physical situation that C 1 C,, at all times. The rate of process c is denoted by J. The form of eq 2 was originally proposed by VolmerI3 and discussed further by Smith et aL9Joand Noyes.I4 Because C appears in eq 2 as a square in the denominator of an exponent, J behaves virtually as a step function as C increases. Process d is described by the M terms which are summed in eq 1, Each j designates the average of a range of bubble sizes reprmnted by a population N of bubbles per cm3of solution and having internal pressure Pj. 'khe quantity k,,in the summation is equal to kex/V,and K is the Henry's law constant for gas solubility. In order to describe process d with the use of eqs 1 and 2, it is necessary to write another M - 1 independent equations involving the rates of change of the Nj populations of bubbles including the escape of bubbles through the surface. Because the bubbelator was used to model oscillations in pressure, it was also necessary to describe process e by use of an equation which included the effects on the gas pressure of the escape of bubbles of many different sizes and also of a leak to the external atmosphere. The computations" (with M = 60) could be made to reproduce simultaneously the period and amplitude observed experimentally by using a single disposable parameter which increased k,,A above the value observed through the smooth surface. A possible justification of this parameter is that during part of each cycle, the rapid evolution of bubbles increases the area of gas-liquid interface above that at other times. No attempt was made to model that effect quantitatively. Although we are confident that the computations of YRN do indeed describe the experimental system, the computations are complicated by the need to include M different bubble sizes which range by orders of magnitude. Too small a value of M will make the model of eqs 1 and 2 too coarse, and too large an M will make the equations stiff and unwieldy. YRN" selected a compromise M 60 for most of their computations, but that compromise
seemed to be only marginally satisfactory. It would be desirable to develop an alternative model which reduced the number of independent variables even at the expense of some of the experimental details. Smith and Noyeslo first suggested that such an approximate model might be developed with the introduction of time delay into the equations. Instead of over 60 simultaneous equations, it might be possible to describe variations of the system satisfactorily by means of 2 equations which treated the concentration of dissolved molecules and presure of gas above the solution as the only variables. C. Examples of Uses of Differential Delay Equatiolrs. Many studies have already been made where the introduction of a delay was found to offer a useful way to treat the development of a system with time. In studies of population biology and ecology, delays of the order of days or even years may arise from the need for a young generation to mature.15J6 Respiration provides a physiological example. Accumulation of carbon dioxide in the lungs causes transmission of a signal to the brain which results in a delayed increase in the rate of ventilation.17J8 Instabilities and even chaotic behavior may r e s ~ l t . ' ~ , ~ In biochemistry, GoodwinZ1discussed a model for enzyme control and proposed a delay equation which used diffusion to cause repression by negative feedback. Bailey and Prairie discussed how feedback could affect the dynamic behavior of a systemZZand applied their treatment to model the hydrogenation of ethylene on a heterogeneous cataly~t.2~ Many examples already exist of applications to homogeneous chemical systems. Thus, Weiner et al.24introduced delay as a way to control the minimal bromate oscillator, and the same system was also investigated recently by Chevalier et alaZ5Ross and ~ o - w o r k e r s studied ~ ~ , ~ ~ thermochemical dimerization and isomerization reactions in which the light intensity falling on the system was dependent with a delay on the concentration of a species. Cheng et alaZ8used delay to simplify a system of ordinary differential equations much as we report here. Two very recent papers by Epstein and LuoZ9Joprovide good discussions of the application of differential delay equations to chemical problems. The extensive interest in such models has encouraged us to explore the suggestion first made by Smith and Noyeslo for the modeling of gas-evolution oscillators by similar techniques. 11. Development of Equations
The equation we develop here makes no effort to treat the range of individual bubble sizes as independent species. Instead, we note that a finite time is necessary for small bubbles to grow to large ones. Infinitesimal nuclear bubbles are formed at a rate J determined by the instantaneous concentration of dissolved molecules and presented in eq 2. Those nuclei will grow and at some time after formation will rise rapidly to the surface and discharge their contents to the surrounding atmosphere. At the time of discharge, the pressure in a grown bubble will be little more than that of the atmosphere, and the concentration in the surrounding solution will be only moderately more than Cut. On the basis of this model, we then write eq 3. dC(t)/dt = Rch - k n ~ [ C (-~c)s a t l e x p h / [ C ( t - 7 ) - c,,I2] (3) In most of the experiments which have been reported, gas pressure changed rather little during an experiment of interest so that C, could be approximated as a constant. If we then adopt the definition x = C - C,,,, we obtain eq 4. dx,/dr = Rch - knuxrexp(-a/x,-:)
(4)
In eq 4, x, refers to the concentration of dissolved molecules at time t , and x,-, refers to the concentration of such molecules which existed at time T before t. The disposable parameter T can in principle be assigned any positive value between zero and infinity.
Bar-Eli and Noyes
7666 The Journal of Physical Chemistry, Vol. 96, No. 19, 1992
The term Rchis the same as the first term in the right-hand side of eq 1. Often Rch will change sufficiently slowly that we can regard it as constant for a considerable time. However, Rch in any closed system must ultimately decay to zero because the chemicals arc consumed. Changes in the rates of chemical reaction can be accommodated if Rchis replaced by R@'. Such a replacement assumes the reaction is fmt order, and different kinetics could be used if the necessary information were available. The remaining term in eq 4 replaces all of the other terms on the right-hand side of eq 1 including the substitution based on eq 2. The rate constant k,, is not rigorously defined in terms of experimentally measurable quantities. Equation 4 describes the change of the single variable C and invokes the disposable parameter 7; it replaces M + 1 independent equations involving as many variables summarized by eqs 1 and 2 above. Of course, eq 4 is an approximation, but computations with it will be much easier than with M 1 simultaneous stiff equations. The objective of the present paper is to show that computations based on the single equation (4) can simulate a wide range of experimental observations which would otherwise require much more complexity. Equation 4 is written to describe the variation with time of the concentration of dissolved molecules. Although that concentration could in principle be measured by spectroscopicprocedures, experimental measurements in the past have been of the pressure P above the solution. The negative term for decreasing the concentration in eq 4 will cause a corresponding increase in pressure. Any decreases in pressure will occur because of escape of gas to the atmosphere through a physical leak. We therefore write eq 5.
+
dP,/dr = k2x, exp(-a/x,-:)
- kn(P,- P o )
(5)
The first term in eq 5 describes the increase in pressure in the system due to escape of gas from the solution, while the second term describes the escape of gas from the system to the external atmosphere of constant pressure Po through the leak if it exists. Here kflis the first-order rate constant for flow of gas from the space above the solution to the external atmosphere at pressure Po. Of course, kfl will be zero if the container is closed. We wish to normalize eqs 4 and 5 and to write them in terms of dimensionless variables. We do so with the substitutions of e q ~6-10. x' = x/(y1/2
(6)
t' = Rcht/crll2
(7)
period of interest. For this simplified form, x is the single variable while k and T are adjustable parameters whose relative values determine the behavior of x with time. Equation 11 has one nonzero steady state, x,, for which dxldt = 0 and x,-, = x , at all times. That steady state will satisfy eq 13, and for any k, one can calculate the appropriate x,. Note x, exp(-1 /xM2) = 1 / k
that the value of x at the steady state depends only upon k (based upon Rch and knu)but is independent of the value of 7 . We have followed the treatment of Murray16 to explore the stability of the steady state. Let 6x be given by eq 14. 6x = x - x, (14) If 6x is small enough that 6x = dx/dt can be described by terms linear in 6x, we obtain 6x = aSx,, + b6x, (15) where a = -2/x,3
6x = cehf
dx/dr = 1
- kx, exp(-l/x,-:)
dp/dt = kx, exp(-1 /x,-,2) - k l p
(11) (12)
Equations 1 1 and 12 are approximations. For instance, they treat C,, as constant over the time of interest. Such an assumption will not be valid indefinitely in a closed container. Similarly, eqs 1 1 and 12 as written were derived with the assumption that Rch is a constant. The discussion of eq 4 above concerns how that assumption might be modified. Because the behavior of eq 1 1 is not affected by changes of pressure, it can be solved independently of eq 12. Equation 12 involving the experimentally measured pressure can then be solved with the use of values of x obtained from eq 1 1 .
III. Stability of the Steady State of Eq 11 Equation 1 1 as written here treats k as a constant and thereby assumes that Rchand C ,, do not change significantly during the
(18)
(19) For any value of k, the value of x, at the steady state can be calculated with eq 13, and the real and imaginary parts of h can then be calculated by means of eqs 20 and 21. h=p+iw
p
= ae-#' cos
(UT)
w = -ae-fiT sin
+b
(UT)
(20) (21)
If CI < 0, the system will relax (with damped oscillations if p is near enough to zero) to a stable steady state x,. If p > 0, x, will be unstable and the system will relax to limit cycle oscillations. The bifurcation between the stability and instability of x, will occur when p = 0. Let a subscript B denote a value of T or o at the bifurcation point where p = 0. Equations 22 and 23 then follow from (20) and (21), and oBand T B are calculable from (24) and (25). = -b/a
(22)
= -wB/a
(23)
wB = (a2 - b2)1/2
(24)
sin
(9)
If the primes are dropped on these definitions, we then obtain normalized dimensionless equations (1 1) and (12) as equivalent to (4) and (3, respectively.
(16)
b = -l/xSs (17) Near the steady state, the solution is given by eq 18, where h is the complex number of eq 19.
COS (wBTB)
p' = (P - P o ) k , , / k 2 ~ ~ ' / ~
(13)
78
(wBTB)
= (l/We)
COS-'
(-b/a)
(25)
Because x, and b/a must both be positive, the trigonometric form of eq 25 requires that 0 Ib/a I1. The bifurcation between stable and oscillatory solutions then comes when a = b and xW2 = 2. Then k , in eq 26 is the minimum possible k at which an oscillatory solution of eq 1 1 is possible. k , = (e/2)'I2 = 1.165822 (26)
For any value of k > k,, we can use eqs 24 and 25 to calculate the minimum value of T above which sustained oscillations are possible. In logarithmic Figure 1, the dashed vertical line designates the value of k,, and the full cutye shows T B as a function of k. Combined values of k and T to the left of and below the curve in Figure 1 designate systems which will decay to a stable steady state (with damped oscillations if they are close to that curve). Combined values represented by points above and to the right of that curve will exhibit sustained oscillations. Those oscillations will be approximately sinusoidal for points near that curve but will become of relaxation type for points further from the curve. When k = k, and T is large enough that willations are possible, TW will be close to T . When k is very large, TO will be close to ~ / 2 .Therefore, the ratio of the period, 2r/w, to the delay, 7 , T ~ ,
The Journal of Physical Chemistry, Vol. 96, No. 19, 1992 7667
Gas-Evolution Oscillators
\ I
,
, ,
,
I
,
,
,
,
10
,
I
k
,
I
,
,
I
,
v
I
\ , , ,
100
180
1000
TABLE I: Parameter Values for Bifurcation Points Associated with Different Rate Constants m
50.9531 16.6550 3.3031 1.9923 1.0582 0.6014 0.3432 0.1 103 0.0563 0
m
1.4142 1.4117 1.3941 1.2216 1.1160 0.9683 0.8361 0.7134 0.5050 0.4079 0
m
104.7228 36.0451 8.6010 5.5809 3.2296 1.9599 1.1795 0.4078 0.2140 0
2.0553 2.1642 2.6039 2.8012 3.0519 3.2591 3.4370 3.6990 3.7986 4
'Rate constant. bMinimum delay for oscillations. CSteady state. "Period at bifurcation. 'Period/delay.
will vary only by a factor of 2 over the entire range of possible k values. These various points are illustrated by the entries in Table I. If the equations lead to oscillatory solutions for a particular set of k and 7 , oscillatory solutions will also be obtained for a different f such that f = T + 2 4 w and n is an integer. The value of w is unchanged in such a solution, and no new phenomena are created by the different solutions.
IV. Representative Numerical Solutions There are various methods to solve differential delay equations. Some of those methods are described by Bellman and Cooke3' and in the recent articles by E p ~ t e i n . ~In, ~our ~ case, the explicit Euler method sufficed while keeping previous data for the delay time to use in the exponential. The initial conditions used were usually p = x = 0 for times -7 < t C 0. (In order to avoid zero in the denominator of the exponent, a very small number such as 10-" was used instead.) To establish that the oscillations were stable, a few cases were checked with different initial conditions such as x = 2 during the delay time. In all cases, the same oscillations were obtained regardless of the initial conditions. A. Sdutiolrp for Changes of Concentration Only. Equation 11 describes the variation of concentration only and has been normalized with the assumption that &,,remains constant. The solid curve in Figure 2 is for k = 10, while T = 0.4 is somewhat more than the T B from Table I. Just as predicted, the oscillations are approximately sinusoidal about the predicted x, value of 0.713 with a period of 1.35. If TB denotes the period at the bifurcation, computations with different T values only a little more than 7 g showed that T TBwas proportional to 7 - T B just as expected. For such computations, the amplitude (x,,, - xmin) was proportional to ( T - T B ) ' ~ ' ,also just as expected for the neighborhood of a supercritical Hopf b i f u r c a t i ~ n . ~ ~ The dashed curve in Figure 2 shows behavior for the same value of k as the solid curve but for T = 4 or 10 times as large as for the solid curve. Although both curves in Figure 2 relate to the same unstable steady state with x, = 0.713, the form of the dashed curve is relaxational rather than sinusoidal. The total period is
-
la5
e o
495
500
time
Figure 1. (solid curve) Logarithmic plot of the delay at bifurcation, T ~ as function of k for steady-state solutions of eq 11. The vertical dashed line corresponds to k, = 1.1658.
1.1658 1.17 1.2 1.6 2 3 5 10 100 999
- _ - i
0
, , . l ,
,
Figure 2. Variation of x with t for numerical solutions of eq 11 when 6 = 0. For both curves, k = 10. For the solid curve, T = 0.4. For the dashed curve, T = 4. TABLE EI: Behaviors for Delays h n g e r than T~ P P xmax xmin 2*/wc 2 9 7.04 0.511 20.35 3 1.2 1.43 0.665 3.65 9 9.28 0.337 19.95 9.9 10.2 0.337 21.8 5 0.65 1.08 0.63 2.1 1 1.56 0.378 3. 5 5.46 0.207 11.3 9 9.45 0.202 19.5 10 0.344 0.739 1963 0.6876503 1.181 600 0.345 0.7529358 0.673 9766 1.184600 0.346 0.763 1209 0.6638538 1.188000 0.771 5304 0.6555046 0.347 1.190500 0.7788416 0.6482518 1.194000 0.348 0.349 0.785 3888 0.641 7604 1.196 500 0.791 3647 0.635 8383 0.350 1.200000 0.835 2088 0.5924270 1.230000 0.360 0.370 0.8658674 0.5620505 1.261 000 0.8904347 0.5376630 1.291 000 0.380 0.390 0.911 2967 0.5169236 1.321 000 0.9297935 0.498 5279 0.400 1.351 500 1.05 0.381 1.64 0.5 0.6 1.14 0.303 1.95 1.34 0.214 0.8 2.4 1 1.53 0.175 2.86 2.48 0.118 4.95 2 3 3.47 0.109 7.05 4 4.46 0.105 9.1 5.44 0.103 11.1 5 6 6.44 0.102 13.15 7 7.44 0.102 15.2 8.43 0.101 17.2 8 9 9.43 0.101 19.2 9.9 10.32 0.1 21.05 100 1 1.4 0.017 2.6 9 9.35 0.01 18.8 999 9 9.31 0.001 18.6 9.9 10.2 0.001 20.3 Rate constant.
Delay.
2*/7Wd 2.26 3.04 2.22 2.20 3.23 3.00 2.26 2.17 3.4349 3.4336 3.4335 3.4308 3.4310 3.4284 3.4286 3.4167 3.4081 3.3974 3.3872 3.3788 3.28 3.25 3.00 2.86 2.48 2.35 2.28 2.22 2.19 2.17 2.15 2.13 2.13 2.60 2.09 2.07 2.05
Period. dPeriod/delay.
approximately 27, and the system spends about half the time with a small x equal to about l/k. The value of dx/dt then changes abruptly from zero to unity and x rises to a sharp maximum at x,, = 7. The fall from the maximum to the minimum x is very rapid and is more so the larger the value of T . We shall relate this mathematical behavior to the experimental system in section VI below. Table I1 shows the variations in the characteristics of the behaviors obtained for computations of eq l l with the use of various selected values of k and T . Some of the trends exhibited by Table I1 are as follows: (a) When the oscillations are still approximately sinusoidal because 7 is little more than T ~the , amplitude increases only slowly with T . This behavior is discussad more quantitatively above, and some of the entries for k = 10 in Table I1 contain more siflicant figures so that readers can test the validity of the claims made there. (b) As k increases at constant T and as T increases at constant k, the period approaches 27. (c) As k and T both
7668 The Journal of Physical Chemistry, Vol. 96, No. 19, 1992
Bar-Eli and Noyes
j: 260
-
L
3 140-
;120&loo-
80 60 40 -
20-
46-
47C +
480
490
500
5
/-
0
r e
Figure 3. Variations of x (dashed curve) and of p (solid curve) for numerical solutions of eqs 11 and 12, respectively, when k = 10, T = 10, kl = 0.1, and fi = 0.
increase, x,
approaches T and x,, approaches 1/k. B. !Mutiom Related to Pnssure Changes. Equations 11 and
12 contain an identical exponential term which involves a delay. That term is negative in eq 11, which describes how dissolved molecules are removed from solution to enter first bubbles and ultimately the gas phase. The same term is positive in eq 12 which describes how molecules leave the solution and enter the gas. If kl # 0, the negative term in eq 12 relates to leakage of gas from the container to the external atmosphere. Different types of behavior are observed depending upon the magnitude of the rate of leakage, kl. (i) k, = 0. If kl = 0, the container is closed and p increases monotonically. For such a situation, d(x p)/dt = 1 and x + p will increase at a normalized constant rate equal to that of the chemical reaction producing dissolved molecules. If T is only a little more than TB, oscillations in eq 11 will be approximately sinusoidal as with the solid curve in Figure 2. Then p in eq 12 will increase monotonically with wavering slope until it reaches a limit when the chemical reactants have been consumed. If T is much greater than 78, d(x p)/dt will still be constant and equal to unity. However, p will now increase in steps with dp/dt = 0 when the dashed curve in F w e 2 is increasing steeply and Q/dt = 1 when x is small and almost constant. Such behavior is illustrated in Figure 4, which is concerned with modeling experimental behavior in section V below. (U) k l # 0. If the leak to the external atmosphere is finite, d(x p)/dr = 1 - kg. Then the sum of x p need no longer be constant. If the steady state for eq 11 is unstable and oscillatory, p will also exhibit oscillations. If kl is sufficiently small and if T - T B is also small, p may still increase initially. However, termination of the chemical reaction will ultimately cause p to attain a limit and then decrease instead of increasing indefinitely as it would if kl = fi = 0. If kl is somewhat larger and if 7 - T B is still small, p may drift to smaller values with ripples. If T - TB is large enough that oscillations in eq 11 are of the relaxational type and if kl is also sisnificantlylarge, we may obtain the more complicated behavior illustrated with Figure 3. In this fisure, the dashed m e describes the variation of the concentration by q 11 and resembles the dashed curve in Figure 2, while the solid curve describes the variation of the pressure by eq 12. Each period in Figure 3 exhibits two subperiods of significant duration. During one subperiod, dx/dt = 0 and x = l/k; then dp/dt = 1 - kg. This expression integrates to eq 27, where pmXis a maximum point in the oscillatory solid curve of Figure 3.
+
+
+
+
During the second subperiod, x is increasing with dx/dt = 1. Then dp/dt = -k,p. This expression integrates to q 28, where pb is the pressure at the break point in the descending solid curve of Figure 3.
@
20
40
60
80
00
20
time Figure 4. Variation of pressure with time in a closed system with only moderate depletion of reactants. Top plot is experimental for the dehydration of formic acid. Data are unpublished observations by Bowers. Bottom plot is for k = 10, T = 2, kl = 0 (closed system), and fi = 0.05.
Figure 3 was computed for a selected value of k l . Other behaviors may be obtained. Thus, if k, is increased while k and 7 and the dashed curve in Figure 3 remain unchanged, p will fall more rapidly from its maximum to reach a plateau where p = l/kl and will then drop almost discontinuously to near zero before it rises again. If kl is decreased from the value used in Figure 3, oscillations in p will fall to very small values before they start to rise again, but each oscillatory spike will exhibit a broadened %oulder” on the right side. We do not think it is worth illustrating the range of possibilities. C. Solutions Which Include Coaspmptioa of Reactants. The examples discussed in A and B above used eqs 11 and 12 derived with the assumption that was a term describing a constant rate of generation of solute molecules. As was pointed out in section 11, 4,must decay to zero in a real system. One way to handle this situation would be to replace unity with e+‘ in normalized equation (1 1). Such a substitution with a relatively small initial value of 0 would have the effect of causing k to increase slowly while T simultaneously decreased so that kr remained constant. A system initially in the stable region of Figure 1 would then move to the right and down to become Oscillatory. Although the nature of the oscillations would evolve in such a system, the solutions of eqs 11 and 12 during times of rather short duration would not differ much from the situations discussed in A and B above.
V. Comparison of Computations with Experimental Observations Equations 11 and 12 involve only two variables and in their normalizbd form involve the three disposable parameters k, T , and kl with fi used to describe depletion of reactants at long times. As was pointed out in the Introduction, they replace a previous effort at modeling which used 62 independent variables. In this section, we shall examine how effectively our equations can model experimental observations. Figures 4-8 each consist of two plots. The upper plot in each figure involves experimental measurements of pressure changes or proportional rates of gas flow as functions of time. Each lower plot is a computed solution of eq 12 for selected values of the parameters. The upper plots have dimensions while the lower plots are normalized and dimensionless; therefore, the plots in Figures
The Journal of Physical Chemistry, Vol. 96, No. 19, 1992 7669
Gas-Evolut ion Oscillators
0
20
60
40
80
100
time
Figure 5. Variation of pressure during repeated bursts of gas evolution in an open system. Top plot is experimental by Bowers and Rawji3 for the dehydration of formic acid. Reprinted with permission from ref 3. Bottom plot is for k = 10, T = 4,k, = 1, and /3 = 0.03.
h
0
500
1500
1000
time
Figure 7. Variation of pressure with time in a more slowly stirred solution. Top plot is experimental by Smith et a1.4 for the dehydration of formic acid a t moderately high concentration. Reprinted with permission from ref 4. Bottom plot is for k = 1, T = 5, k, = 0.1, and fl = 0.003.
"1 I 0
10
20
30
40
50
time
0
100
200
300
LO0
50C
time
Figure 6. Variation of pressure with time in a rapidly stirred solution. Top plot is experimental for the dehydration of formic acid a t somewhat low concentration from p 26 of the thesis of Smith? Bottom plot is for k = 1, T = 9, k, = 0.1, and /3 = 0.015.
4-8 cannot be compared quantitatively but only for form. Figure 4 shows some measurements by Bowerd3 in a closed
system. The experiments were continued for a time during which the rate of consumption of reactants decreased only moderately, while the computations continued until Rchhad fallen almost to zero. Both plots exhibit successive steps which decreased in height and increased in duration as the reaction continued. Figure 5 shows some observations by Bowers and Rawji3during which separated bursts decreased in amplitude because of depletion of reactants. Both plots show similar behavior with almost equally spaced pulses followed by a pressure fall almost to zero. The lower computational plot exhibits shoulders on the right side as discussed in section VB. The experimental plots do not show similar shoulders, but we believe they might be found for a slightly larger size of leak to the atmosphere. Figure 6 shows an experimental plot from Smith's thesis? The solution was rapidly stirred to reduce the total number of pulses, and the rate of chemical reaction decreased greatly during the
Figure 8. Variation of pressure with time for another gas-evolution oscillator. Top plot is experimental by Kaushik et al.34for the decomposition of ammonium nitrite. Reprinted with permission from ref 34. Copyright 1986 Journal of Chemical Education. Bottom plot is for k = 10, T = 0.5, k, = 1, and /3 = 0.03.
time shown. The computations reproduce the form of the observations very well. Figure 7 shows an experimental plot from Smith et aL4 The solution was stirred rather slowly, and the trace is somewhat noisier than the deterministic computations. However, both plots show the onset of oscillations in a reacting solution and the rise and subsequent decay of amplitude to a very small value, while the period slightly increases near the end of oscillations. Experimental plots in Figures 4-7 were all based on the Morgan' reaction involving dehydration of formic acid. The plot in Figure 8 is from the observations by Kaushik et al.34of the decomposition of ammonium nitrite,35 which was subsequently studied by Rubin et a1.* Note that just as in Figure 7, the period remains almost constant with some increase toward the end. The fact that the differential delay equation works so well at simulating Figure 8 supports our conviction that the model we have developed is concerned only with the formation and escape of gas molecules and is irrelevant to the chemical processes involved. We did not make serious efforts to adjust parameters to model experimental observations precisely but were satisfied if most of the significant features could be reproduced. We find it remarkable that such a variety of experimental behavior could be
7670 The Journal of Physical Chemistry, Vol. 96, No. 19, 1992 reproduced so well by varying 50 few parameters. We believe that Figures 4-8 justify the utility of differential delay equations to model systems such as this.
VI. Discussion The objective of this work was to develop a simplified method to model the behavior of a gas-evolution oscillator. The previous detailed effort by Yuan et a1.I' used a bubbelator model which involved over 60 independent variables in as many ordinary differential equations. Sufficient experimental information was available that only one disposable parameter was employed in the computations. The bubbelator modeled experimental observations moderately well, although the large number of independent variables sometimes made the equations uncomfortably stiff. The differential delay equations developed here need only two variables, the oversaturation x and the overpressure p in normalized equations (1 1) and (12). Those equations contain only three parameters. The two parameters k and klcould be evaluated experimentallyfor application to any real system, and the delay T is selected to best simulate observations. Our model neglects transport of molecules through the upper surface of the solution and effects of stirring on that transport. It also neglects the effects of viscosity and of surface tension and the detailed kinetics of bubble formation and growth. We might anticipate that the simplifications resulting from the use of the approximate delay equations would not reproduce observations as exactly as was possible with the much more complicated bubbelator model. However, Figures 4-8 show that by appropriate selection of the three parameters, we can reproduce very well a wider range of experimental behavior than has ever even been attempted with the bubbelator! In spite of this considerable success, we must confess to some uncertainty about the physical significance of the delay T , which is introduced as a single disposable parameter and not precisely defined. Processes in the experimental system generate four periods during a cycle. These periods are (a) the time required to raise the concentration C from a small value to the threshold for homogeneous nucleation, (b) the time required for an initial nucleus to grow enough that the resulting bubble rises to the surface and escapes, (c) the time between when the concentration rises to and starts to overshoot the threshold for nucleation and when the concentration has gone through a maximumand returned to the nucleation threshold from above, thereby stopping the nucleation process, and (d) the time from when nucleation stops to the time when the concentration of dissolved molecules has returned to a very low value while those bubbles present at the start of the period have grown and escaped to gas. These four periods do not describe a clean sequence because some are concerned with the system as a whole, while b in particular is associated with the behavior of a single bubble. These periods have been introduced to show that the development of the system is more complicated than one would expect to be able to describe with a single delay period. Equations 11 and 12 are simplified approximations, and the single delay time T probably cannot be simply correlated with any quantity in the experimental system. In spite of their approximate nature, we remain surprised that these equations model so well such a wide range of experimental behavior. In spite of the simplifications noted above, it is of interest to investigate experimentallythe effects of parameters such as viscosity, temperature, surface tension, nature of the evolving gas,
Bar-Eli and Noyes etc., on the behavior of the system. By comparison of such experimental results with calculated delays which give a best fit to the observations, we may get a better insight regarding the physical processes which occur in gas-evolution oscillators. In collaborationwith Prof. Peter Bowers of Simmons College, we hope to develop thermodynamic expressions to treat nucleation and growth of bubbles as functions of the composition of a solution. Any such results will be reported separately.
Acknowledgment. This work was supported in part by Grant No. CHE8717791 from the National Science Foundation. This paper is No. 97 in the series "Chemical Oscillations and Instabilities". No. 96 is as follows: Ruoff, P.; Fijrsterling, H.-D.; Gyiirgyi, L.; Noyes, R. M. J. Phys. Chem. 1991,95,9314-9320. We are indebted to Prof. Robert M. Mazo of the University of Oregon for helpful discussions during the development of these ideas. We are also indebted to Prof. Peter G. Bowers of Simmons College for permission to use the unpublished observations in Figure 4. References and Notes (1) Morgan, J. S. J. Chem. Soc. 1916, 109, 274-283. (2) Showalter, K.; Noyes, R. M. J . Am. Chem. Soc. 1978,100, 1042-1049. (3) Bowers, P. G.; Rawji, G. J . Phys. Chem. 1977,81, 1549-1551. (4) Smith, K. W.; Noyes, R. M.; Bowers, P. G . J. Phys. Chem. 1983, 87, 1 5 14-1 5 19. (5) Bowers, P. G.; Noyea, R. M. J. Am. Chem. Soc. 1983,105,2572-2574. (6) Bowers, P. G.; Dick, Y. M. J . Phys. Chem. 1980,84, 2498. (7) Ganapathisubramanian, N.; Noyes, R. M. J . Phys. Chem. 1981,85, 1103-1105. (8) Rubin, M. B.; Noyes, R. M.; Smith, K. W. J. Phys. Chem. 1987, 91, 1618-1 622. (9) Smith, K. W. Ph.D. Thesis, University of Oregon, Eugene, 1981. (10) Smith, K. W.; Noyes, R. W. J. Phys. Chem. 1983,87, 1520-1524. (11) Yuan, Z.; Ruoff, P.; Noyes, R. M. J . Phys. Chem. 1985, 89, 5726-5732. (12) Kaushik, S. M.; Rich, R. L.; Noyes, R. M. J. Phys. Chem. 1985.89, 5722-5725. (13) Volmer, M. Kinetik der Phasenbildung; Steinkopf: Leipzig, 1939. (14) Noyes, R. M. J. Phys. Chem. 1984,88, 2827-2833. (1 5) MacDonald, N. Time Lugs in Biological Models; Springer-Verlag: Berlin, 1978. (16) Murray, J. D. Mathematical Biology; Springer-Verlag: Berlin, 1989; p 12. (17) Banks, H. T. In Lecture Notes Appl. Math. 1981, 19, 361-383. (18) Banks, H. T. In Nonlinear Systems and Applications, an International Conference; Lakshmikantham, V., Ed.; Academic Press: New York, 1976; pp 21-38. (19) Mackey, M. C.; Glass, L. Science 1977, 197, 287-289. (20) Glass, L.; Mackey, M. C. Ann. N.Y. Acad. Sci. 1979,316,214-235. (21) G d w i n , B. C. Adu. Enz. Reg. 1965, 3, 425-438. (22) Prairie, M. R.; Bailey, J. E. Chem. Eng. Sci. 1987, 42, 2085. (23) Lyberatos, G.; Kyszta, B.; Bailey, J. E. Chem. Eng. Sci. 1984.39.947. (24) Weiner, J.; Schneider, F. C.; Bar-Eli, K. J. Phys. Chem. 1989, 93, 2704-27 11. (25) Chevalier, T.; Freund, A.; Ross, J. J . Chem. Phys. 1991,95,308-322. (26) Zimmermann, E. C.; Schell, M.; Ross, J. J. Chem. Phys. 1984,81, 1327-1 336. (27) Schell, M.; Ross, J. J . Chem. Phys. 1986,85, 6489-6503. (28) Cheng, Z.; Hamilton, I. P.; Teitelbaum, H. J. Phys. Chem. 1991,95, 4929-493 1. (29) Epstein, I. R. J . Chem. Phys. 1990, 92, 1702-1712. (30) Eptein, I. R.; Luo, Y. J. Chem. Phys. 1991, 95, 244-254. (31) Bellman, R.; Cooke, K. L. Differential Difference Equations; Academic Press: New York, 1963. (32) Hassard. B. D.; Kazarinoff, N. D.; Wan, Y. H. Theory and Applications of Hopf Bifurcation; London Mathematical Society Lecture Notes Series 0076 0552 41; Cambridge University Press: Cambridge, 1981. (33) Bowers, P. G. Private communication. (34) Kaushik, S. M.; Yuan, Z.; Noyes, R. M. J. Chem. Educ. 1986,63, 76-80.
