Conditions for Critical Effects in the Mass Action Kinetics Equations for

Article pubs.acs.org/JPCA

Conditions for Critical Effects in the Mass Action Kinetics Equations for Water Radiolysis Richard S. Wittman,* Edgar C. Buck, Edward J. Mausolf, Bruce K. McNamara, Frances N. Smith, and Chuck Z. Soderquist Energy and Environment Division, Pacific Northwest National Laboratory, Richland, Washington 99352, United States ABSTRACT: We report on a subtle global feature of the mass action kinetics equations for water radiolysis that results in predictions of a critical behavior in H2O2 and associated radical concentrations. While radiolysis kinetics have been studied extensively in the past, it is only in recent years that high-speed computing has allowed the rapid exploration of the solution over widely varying dose and compositional conditions. We explore the radiolytic production of H2O2 under various externally fixed conditions of molecular H2 and O2 that have been regarded as problematic in the literaturespecifically, “jumps” in predicted concentrations, and inconsistencies between predictions and experiments have been reported for α radiolysis. We computationally map-out a critical concentration behavior for α radiolysis kinetics using a comprehensive set of reactions. We then show that all features of interest are accurately reproduced with 15 reactions. An analytical solution for steady-state concentrations of the 15 reactions reveals regions in [H2] and [O2] where the H2O2 concentration is not uniqueboth stable and unstable concentrations exist. The boundary of this region can be characterized analytically as a function of G-values and rate constants independent of dose rate. Physically, the boundary can be understood as separating a region where a steady-state H2O2 concentration exists from one where it does not exist without a direct decomposition reaction. We show that this behavior is consistent with reported α radiolysis data and that no such behavior should occur for γ radiolysis. We suggest experiments that could verify or discredit a critical concentration behavior for α radiolysis and could place more restrictive ranges on G-values from derived relationships between them.

■

INTRODUCTION It is notable that the basis of the chemical kinetics equations, the Law of Mass Action, dates back to the 1860s1,2before the atomic theory of interacting molecules was established. As in thermodynamics, experiment and theory led to general principles that have mostly survived modern developments in physics and chemistry. This survivability is not surprising because the laws of mass action and thermodynamics are primarily the consequence of the statistical nature of many interacting degrees of freedom in our case, interacting molecular species. In radiolysis, many of the species are very short-lived, rapidly de-exciting or interacting with surrounding species. On the time scale of interest for radiolysis kinetics reactions, the formation of short-lived radical species and their reactions in water have been extensively studied.3−8 As in the case of many other authors,3−9 we take these reactions and generation rates as a starting point, but we emphasize that we are not competing with those works or proposing a more sophisticated radiolysis model. Our goal here is to highlight a subtle numerical behavior of the radiolysis kinetics equations in general that has not previously been recognized and then give the analytical basis for such behavior. Whether such behavior actually occurs in nature is an open question to be answered experimentally. As in previous work, we notice a particular dependence of H2O2 production on the concentrations of O2 and H2.5−7 The conditions considered are particularly applicable to the H2O2driven corrosion of spent nuclear fuel in an environment © 2014 American Chemical Society

depleted of O2 and in the presence of H2 overpressure generated from structural iron reacting with water.10 While many complex chemical and physical processes can be imagined to operate at the exposed UO2 surface,9−13 we focus on radiolysis of water from a uniform α dose. Unlike previous work, we identify a reduced reaction set exhibiting the same global features analytically and then show how [H2O2] behaves like an order parameter that characterizes distinct phases of the system. A critical point in H 2 O 2 concentration is predicted that implies a ([H2], [O2]) region with a discontinuous boundary that may be the source of disagreement between model predictions and data.7 This work attempts to explain how a ([H2], [O2]) region with a boundary, discontinuous in [H2O2], can emerge in radiolysis kinetics. The next section describes the approach taken to solve the radiolysis kinetics equations numerically and outlines the reasoning for working with a reduced reaction set analytically at steady-state. It is then shown explicitly how a critical jump in H2O2 concentration arises from stability conditions on model parameters. Finally, we summarize the conditions for global critical jumps in radiolysis kinetics and offer suggestions on how experiments could confirm or possibly correct model parameters. Received: October 1, 2014 Revised: November 24, 2014 Published: November 25, 2014 12105

dx.doi.org/10.1021/jp509856g | J. Phys. Chem. A 2014, 118, 12105−12110

The Journal of Physical Chemistry A

■

Article

state as a function of H2 and O2 concentration for a 5 MeV α dose rate of 25 rad/s. As expected, the numerical solution (Figure 1) took much longer to evaluate (out to t = 108 s) and is shown on a slightly coarser mesh. The initial H2O2 concentration was assumed to be zero for Figure 1 (left) and is compared with the smallest [H2O2] values that solve the analytical model [Figure 1 (right)]. Setting the initial H2O2 concentration to 0.5 M resulted in Figure 2 (left) solutions to be compared with the greatest values that solve the analytical model [Figure 2 (right)]. The steady-state solution, which used the reduced set of equations and the approximations and techniques of the next section, agrees with the numerical integration of the full set of equations (Figures 1 and 2). Both upper and lower solutions to the kinetics equations show a jump in H2O2 concentration at a critical boundary in H2 and O2 concentrations. Also, as the initial H2O2 concentration is continuously lowered from 0.5 M, the numerical steady-state [H2O2] jumps from Figure 2 solution down to Figure 1 solution with no intermediate steady-state values. Figures 1 and 2 therefore show that the smallest and largest [H2O2] values that solve the analytical approximation correspond to different initial conditions assumed for the kinetics equations. Additionally, a continuous surface can be mapped-out (Figure 3) for the analytical approach revealing a third intermediate solution value of [H2O2] not accessible to the time-dependent solution. The black line on the surface of solutions at an α dose rate of 25 rad/s (Figure 3) indicates a separation between a region having a single steady-state [H2O2] and one having three possible values. Because not all values are accessible to the evolving kinetics solutions the nature of the [H2O2] jump would not be obvious in models that simply solve the time-dependent equations. While the continuous surface of Figure 3 gives some insight into how critical jumps in concentration can occur in the kinetics equations, even more can be inferred by extracting analytical relations implied by the reduced reaction set. In the next section, we show how the analytical approximation was obtained to give relations that characterize the critical boundary. Analytical Model Solution and Interpretation. Assuming that the large number of calculations summarized in Figures 1 and 2 support the claim that the reduced reaction set is practically equivalent to the full set, the model solution steps described here allow us to extract useful relationships that can give insight into the behavior of radiolysis kinetics equations and can help interpret experiments. Table 1 reactions imply seven independent steady-state equationsone for each time-dependent species. Three species can be simply eliminated algebraically. Additionally, we numerically confirm that setting k13, k14, k29, and k36 to zero has negligible effect on the solution because (1) reactions 14 and 29 (see Table 1) generate products of small reactant concentrations and (2) reactions 13 and 36 (see Table 1) have small rate constants. This simplification allows a direct solution of [e−], thus reducing the number of equations to three (eqs 1−3). The remaining three equations for [H2O2], [·OH] and [·H] are

RADIOLYSIS KINETICS ANALYSIS Taking the conventional reactions describing water radiolysis, we solve the kinetics equations numerically to compare with analytical steady-state solutions derived from a reduced reaction set. After showing that the important features of the full model are replicated by the reduced reaction set, we derive analytical relations that describe the observed boundaries of critical concentrations. Model Description. There are many investigations of water radiolysis that contain reaction rate constants and G-values appropriate to the radiation environment considered here. Most of them refer to radiation chemistry databases (e.g., NIST14) or the original work of refs 3 and 4. To be specific, our full model reaction set is identical with Pastina and LaVerne7 with the addition of a single reaction for the direct thermal decomposition of H2O2.15 Even though at 25 °C such a direct decomposition reaction physically occurs at a very slow rate, its inclusion ensures the possibility of steady-state under all conditions. In fact, the critical conditions can be understood as a parametric boundary of a region where no steady-state exists without including direct H2O2 decomposition. Numerically, we know such regions exist with one indicated for α radiolysis in ref 7. We numerically confirm all reported results of ref 7 by solving the kinetics equations in FORTRAN with the stiff ordinary differential equation (ODE) solver DLSODA of ODEPACK.16,17 Additionally, solution consistency over time is verified by checking that charge balance and atom balance are preserved to within machine precision. The rate constant numbering of Table 1 corresponds to the 79 reactions of ref 7 with number 80 used for the H 2 O 2 Table 1. Subset of Reactions Sufficient To Represent Full Radiolysis Model Predictions with Reaction Rate Constants (Mn s−1) at 25°C

a

reaction no.a

reaction

kr

3 10 13 14 21 22 26 29 36 42 43 50 51 57 80

H2O2 → H+ + ·HO2− H+ + ·HO2− → H2O2 e− + H2O → ·H + OH− ·H + OH− → e− + H2O ·HO2 → O2− + H+ O2− + H+ → ·HO2 e− + H2O2 → ·OH + OH− e− + O2 → O2− ·H + H2O → H2 + ·OH ·H + H2O2 → ·OH + H2O ·H + O2 → ·HO2 ·OH + H2 → ·H + H2O ·OH + H2O2 → ·HO2 + H2O ·HO2 + O2− → ·HO2− + O2 H2O2 → ·OH + ·OH

1.1 × 10−1 5.0 × 1010 1.9 × 101 2.2 × 107 1.3 × 106 5.0 × 1010 1.1 × 1010 1.9 × 1010 1.1 × 101 9.0 × 107 2.1 × 1010 4.3 × 107 2.7 × 107 8.0 × 107 2.5 × 10−7

Reaction numbering of ref 7.

decomposition reaction. We are specifically interested in a system where H2 and O2 concentrations are set by the external environment. Experimentally, such a system could be explored with an aqueous system in continuous contact with a gas of controlled partial pressures of H2 and O2. The full model numerical steady-state solutions and an approximate analytical solution derived from the equations based on reduced reaction set of Table 1 (with derivatives set to zero) agree within a few percent and are visually identical. Figure 1 shows both numerical (left) and analytical (right) based H2O2 concentrations at steady12106

−k43[O2 ][·H] + (k51[H 2O2 ] + k50[H 2])[·H] 1 = G̃ H2O2 + (G̃ ·OH + G̃ ·OH2 − G̃ ·H − G̃ e−) 2

(1)

(k43[O2 ] + k42[H 2O2 ])[·H] − k50[H 2][·OH] = G̃ ·H

(2)



Article

Figure 1. Comparison of steady-state [H2O2] at 25 rad/s (5 MeV α) for the full radiolysis model ([H2O2]0 = 0 M) (left) and for the smallest analytical solution to Table 1 reaction set (right).

Figure 2. Comparison of steady-state [H2O2] at 25 rad/s (5 MeV α) for the full radiolysis model ([H2O2]0 = 0.5 M) (left) and for the greatest analytical solution to Table 1 reaction set (right).

Gḋ). Equations 1 and 2 result from combinations of the original seven rate equations at steady-state, but notice that eq 3 is the time derivative of [·OH] and can be used to assess the stability of the ·OH concentration around steady-state. Eliminating [·OH] and [·H] gives a single cubic equation for [H2O2] according to A 0 + A1[H 2O2 ] + A 2 [H 2O2 ]2 + A3[H 2O2 ]3 = 0

where [·OH] and [·H] can be calculated from eqs 1−3. For example, using eqs 1−3 to instead solve for [·OH], we can plot d[·OH]/dt (eq 3) as a function of [·OH] to determine how the rate of change of [·OH] depends on small perturbations from steady-state. Of course, if at steady-state the slope is negative, the change will return [·OH] to its steady-state value; and, if it is positive, the change will force [·OH] away from its steady-state value indicating an unstable solution. Figure 4 corresponds to one ([H2], [O2]) location in Figure 3 where the three solutions for [H2O2] are indicated with arrows. In analogy with potential theory, the red curve (Figure 4) is the negative integral of d[·OH]/dt with respect to [·OH], showing the high and low solutions as stable “potential” minima with the center solution as an unstable maximum. This feature is consistent with the steadystate behavior of the full time-dependent model where the center solution is not realized. The explicit forms of the coefficients of eq 4 depend on G-values, dose rate, and rate constants and are given in eqs 5−8 according to

Figure 3. Surface of [H2O2] showing a region on the near side of the black line where three possible values solve the analytical radiolysis model based on Table 1.

−2k51[H 2O2 ][·OH] − k50[H 2][·OH] + 2k 80[H 2O2 ] 1 + G̃ H2O2 + (3G̃ ·OH + G̃ ·HO2 + G̃ ·H + G̃ e−) 2 =0

(4)

(3)

The “G with tilde” notation is short-hand for the radiolytic generation G-values with a factor of the dose rate included (G̃ ≡ 12107



Article

[O2 ] ≥ S−[H 2]

(11)

where in the lowest order (zeroth order) of [O2]/[H2], the boundary slope S− is determined to be. k42k50 (G·OH + G·H + Ge−)2 k43k51

S− =

/[(2G H2O2 − G·OH + G·HO2 − 3Ge−) (2G H2O2 + 2G·OH + G·HO2) − G·H(2G H2O2 + 4G·OH + G·HO2 + 2G·H + 2Ge−)]

Equation 11 (equality) can now be understood as an analytical approximation for the critical boundary of Figure 1 (or black line of Figure 3). In fact, as k80 decreases to zero, the (Figure 1) H2O2 concentration on the high [O2] side of the boundary continues to increase, resulting in an ([H2], [O2]) region where no steadystate solution exists for k80 = 0. As mentioned, Figures 1−3 assume G-values specifically for α radiation where GT = 1.0 molecule/100 eV and S− = 7.6 × 10−4. We now ask if the same critical behavior exits for G-values of other radiation fields. Table 2 shows the water radiolysis G-values appearing in eq 10 for the radiation fields considered in ref 7. Notice that eq 10 test condition (GT > 0) is satisfied only for α radiolysis.

Figure 4. Equation 3 derivative of [·OH] as a function of [·OH] (black) showing the values of [H2O2] corresponding to Figure 3. The red curve is minus the integral with respect to [·OH].

A0 =

k43k50 (2G̃ H2O2 + G̃ ·HO2 + G̃ ·H − G̃ e− + G̃ ·OH)[H 2] 2 [O2 ]

A1 =

(12)

(5)

Table 2. Water Radiolysis G-Values of Eq 10

k43k51 (2G̃ H2O2 + G̃ ·HO2 + G̃ ·H − 3G̃ e− − G̃ ·OH)[O2 ] 2 − k42k50(G̃ ·OH + G̃ ·H + G̃ e−)[H 2] (6)

G-values (molecules/100 eV) species

5 MeV-α

γ

10 MeV-p

2 MeV-p

0.10 1.00 0.10 0.35 0.15 1.00

0.66 0.70 0.02 2.70 2.60 −11.1

0.57 0.74 0.03 1.18 0.90 −4.08

0.20 0.76 0.05 0.63 0.30 −0.56

·H H2O2 ·HO2 ·OH e−

A 2 = −2k 80(k43k51[O2 ] + k42k50[H 2]) k k + 42 51 (2G̃ H2O2 + G̃ ·HO2 − 3G̃ ·H − 3G̃ e− − G̃ ·OH) 2

GT

(7)

A3 = −2k42k51k 80

Table 2 implies that a steady-state solution exists for the other radiation fields even without including the H2O2 decomposition reaction. Furthermore, the critical behavior seen in Figures 1−3 seems to be special to α radiolysis. Table 2 shows an increasing value of GT with linear energy transfer (LET). It would be interesting to determine if a similar critical behavior would be predicted for ion beams with similar or greater LET than for 5 MeV α’s. It is possible that a more microscopic analysis could reveal that the extent of LET could determine if and when a characteristic correlation length can become finite. At the level of mass action kinetics with G-values, this analysis implies that for α radiolysis we can consider [H2O2] as the order parameter of the steady-state system. At a fixed temperature of 25 °C, the critical point of the system for a dose rate of 25 rad/s is given in Table 3. Figure 5 shows the critical point (black dot) on a phase diagram where the red dot indicates the location considered in Figure 4. Notice that for concentrations greater than the critical point, [H2O2] can change continuously from a value greater than

(8)

Equations 4−8 now allow us to derive at least approximate relationships between steady-state concentrations and model parameters. For example, in the limit of zero dose rate it is straightforward to show that lim

̇ 0 d→

[H 2O2 ]2 = k43k50(2G H2O2 + G·OH + G·H + G·HO2 ḋ

− Ge−)[H 2][O2 ]/[4k 80(k43k51[O2 ] + k42k50[H 2])]

(9)

5

where, as is known, the H2O2 concentration is seen to go as the square root of the dose rate (ḋ). Here it should also be mentioned that if we had not set k29 to zero, eq 4 would have been fourth order in [H2O2] with the higher order coefficients containing the expressions of eqs 5−8, and an additional lowest order positive term having an overall factor of k29. Elementary algebra implies that if the H2O2 decomposition rate k80 goes to zero there is the possibility that no solution exists for a test condition GT defined according to GT ≡ 2G H2O2 − 3Ge− − G·OH − 3G·H + G·HO2 ≥ 0

Table 3. Critical Concentrations at ḋ = 25 rad/s

(10)

species

concn (M)

This condition is practically true even when k29 is set to zero because it is unlikely that A0 < 0 (eq 5). Given that eq 10 condition is true, it can be shown that no solution exists when

[H2]c [O2]c [H2O2]c

3.650 × 10−3 4.206 × 10−6 7.284 × 10−3

12108



Article

Figure 6. Comparison of [H2O2] model predictions near critical values of [O2] with data of ref 7.

Figure 5. Critical point in ([H2], [O2]) shown with the black dot where red curves indicate first-order transitions in [H2O2] and the black curve indicates continuous change in [H2O2]. The red dot is the concentration point of Figure 4.

parameter uncertainties dominate. The current model parameters imply that if submicromolar concentrations of O2 can be maintained at 1-atm H2, then the H2O2 concentration should in fact agree with the black dashed curve of Figure 6. It is true that the agreement obtained for [O2] = 1 μM does not require a discontinuity, but very low H2O2 concentrations for the suggested measurements, where [O2] < 0.5 μM, would tend to support the model parameters and indicate a discontinuity. If such data are still in disagreement with predictions, then eqs 10−12 can be used to adjust the allowed ranges of G-values supporting the data.

[H2O2]c to less than [H2O2]c. Discontinuous changes in [H2O2] are implied for changes that cross both the red curves (Figure 5)i.e., approaching the red dot from the right gives the lower value of [H2O2] and approaching the red dot from the left gives the higher value of [H2O2] shown in Figure 4. As understood from eqs 10−12, the upper critical boundary is a boundary above which no steady-state exists without thermal decomposition. Therefore, the region of Figure 5 where [H2O2] < [H2O2]c can be understood as a phase where radical destruction regulates [H2O2], and the region where [H2O2] > [H2O2]c is a phase where thermal decomposition regulates [H2O2]of course both mechanisms operate in both regions. At least within eqs 1−8 approximation, a critical exponent of 1/3 can be determined for both reduced [H2] and [O2] according to [H 2O2 ][H2] = [H2]c − [H 2O2 ]c ∝ ([O2 ] − [O2 ]c )1/3

■

CONCLUSIONS Under the driven conditions of α radiolysis of water the kinetics equations predict the existence of a critical point in H2 and O2 concentrations for steady-state [H2O2] that has not previously been noted. The critical point concentrations are shown to scale with dose rate. This prediction is not inconsistent with current α and γ radiolysis measurements, but whether such discontinuous concentration boundaries occur physically for α radiolysis is an open question. If such a critical behavior is discredited by data then new restrictive ranges on G-values can be inferred from derived relationships presented here. Of course, the conventional radiolysis parameters would be supported if concentration jumps in H 2O2 can be observed in experiments at very low (micromolar) O2 concentrations. If the latter is the case, then we predict that similar effects should be observable for radiation having even greater LET than 5 MeV α’s. It would be interesting to explore the possibility of analogous critical behavior occurring in other mass action kinetics model applications such as atmospheric physics and chemistry, radiation medicine, and health physics.

(13)

and [H 2O2 ][O2] = [O2]c − [H 2O2 ]c ∝ ([H 2] − [H 2]c )1/3

(14)

In eqs 13 and 14 the exponent of 1/3 seems to be confirmed numerically with the full model and the location of the critical point is within about 5% of Table 3 values. One further implication of eqs 1−8 approximation is the universality of the critical conditions shown in Figures 1−5 and Table 3. In all figures, the concentrations of H2O2, H2, and O2 scale with dose rate−e.g., the shape of Figure 3 surface is the same for all dose rates with a simple scale of axes. In considering the measurements of ref 7, Figure 6 shows the H2O2 concentration as a function of cumulative dose where the black curves assume a zero initial O2 concentration. With no added H2, the solid curve and closed diamonds are in agreement. For the saturated H2 case the model prediction (black dashed curve) and measurement (open diamonds) are in strong disagreement. Consistent with the critical boundary implied by eq 11, a jump in [H2O2] is predicted for an O2 concentration of about 0.5μM. Figure 6 indicates that a reasonable agreement with the measurement could be established if [O2] ≈ 1 μM. While it is plausible that micromolar levels of O2 were present in the experimental measurements represented in Figure 6, further work is needed to determine if experimental or model

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank C. F. Jové-Colón, D. C. Sassani, J. L. Jerden, and W. L. Ebert for useful discussions. This work was supported by the U.S. Department of Energy, Office of Nuclear Energy, Used Fuel 12109



Article

Disposition (UFD)−Repository Science (RS) Program under Contract. Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by Battelle under Contract DEAC05-76RL01830.

■

REFERENCES

(1) Guldberg, C. M.; Waage, P. Studies Concerning Affinity. Videnskabs-Selskabet i Christiana 1864, 35. (2) Guldberg, C. M. Concerning the Laws of Chemical Affinity. Videnskabs-Selskabet i Christiana 1864, 111. (3) Elliot, A. J.; McCracken, D. R. Computer modelling of the radiolysis in an aqueous lithium salt blanket: Suppression of radiolysis by addition of hydrogen. Fusion Eng. Des. 1990, 13, 21−27. (4) Elliot, A. J. A Pulse Radiolysis Study of the Temperature Dependence of Reactions Involving H, OH and eaq− in Aqueous Solutions. Radiat. Phys. Chem. 1989, 34, 753−758. (5) Ershov, B. G.; Gordeev, A. V. A model for radiolysis of water and aqueous solutions of H2, H2O2 and O2. Radiat. Phys. Chem. 2008, 77, 928−935. (6) Pastina, B.; Isabey, J.; Hickel, B. The influence of water chemistry on the radiolysis of the primary coolant water in pressurized water reactors. J. Nucl. Mater. 1999, 264, 309−318. (7) Pastina, B.; LaVerne, J. A. Effect of Molecular Hydrogen on Hydrogen Peroxide in Water Radiolysis. J. Phys. Chem. A 2001, 105, 9316−9322. (8) Wittman, R. S.; Buck, E. C. Sensitivity of UO2 Stability in a Reducing Environment on Radiolysis Model Parameters. MRS Proc. 2012, 1444. (9) Christensen, H.; Sunder, S. An evaluation of water layer thickness effective in oxidation of UO2 fuel due to radiolysis of water. J. Nucl. Mater. 1996, 238, 70−77. (10) Carbol, P.; Fors, P.; Gouder, T.; Spahiu, K. Hydrogen suppresses UO2 corrosion. Geochim. Cosmochim. Acta 2009, 73, 4366−4375. (11) Shoesmith, D. W.; Kolar, M.; King, F. A Mixed-Potential Model to Predict Fuel (Uranium Dioxide) Corrosion within a Failed Nuclear Waste Container. Corrosion 2003, 59, No. 9, 802−816. (12) Jerden, J.; Frey, K.; Cruse, T.; Ebert, W. Waste Form Degradation Model Status Report: Electrochemical Model for Used Fuel Matrix Degradation Rate; Argonne National Laboratory, 2012; 1. (13) Trummer, M.; Roth, O.; Jonsson, M. H2 inhibition of radiation induced dissolution of spent nuclear fuel. J. Nucl. Mater. 2009, 383, 226− 230. (14) Manion, J. A. NIST Chemical Kinetics Database, NIST Standard Reference Database; National Institute of Standards and Technology: Gaithersburg, MD, 2013; 1. (15) Hiroki, A.; LaVerne, J. A. Decomposition of Hydrogen Peroxide at Water−Ceramic Oxide Interfaces. J. Phys. Chem. A 2005, 109, 3364− 3370. (16) Hindmarsh, A. C. ODEPACK, A Systematized Collection of ODE Solvers. Scientific Computing 1983, vol. 1 of IMACS Transactions on Scientific Computation, 55−64. (17) Brown, P. N.; Hindmarsh, A. C. Reduced Storage Matrix Methods in Stiff ODE Systems. Appl. Math. Comput. 1989, 31, 40−91.

12110


Conditions for Critical Effects in the Mass Action Kinetics Equations for

Recommend Documents