- Email: [email protected]
Thermal stability analysis via elecidation of hazardous reaction stoichiometries
Computers Chem. EngngVol. 22. No. 6 pp. 735-745, 1998
Pergamon PII: S0098-1354(97)00256-1
© 1998ElsevierScienceLtd Printed in Great Britain. All rights reserved 0098-1354/98 $19.00+ 0.00
Thermal stability analysis via elucidation of hazardous reaction stoichiometries C. Bruneton", C. Hoff b and P.I. Barton"'* a Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, U.S.A. b SANOFI CHIMIE, 30390 Aramon, France
Abstract Many accidents that occur in the chemical processing industries such as thermal runaways are due to unstable reactive systems. A good understanding of these systems requires knowledge of the stoichiometries of the reactions that lead to instabilities, which are typically decompositions and/or polymerizations. A screening test to elucidate these stoichiometries is presented. It eliminates implausible stoichiometries from a candidate set due to their thermodynamic infeasibility. We use computational chemistry coupled with statistical thermodynamics to estimate the thermodynamic properties of the unfamiliar species typically involved in such reactions. The real life runaway due to product instability in a process to manufacture tetrahydrofurfuryl benzenesulfonate is treated as a case study. Hartree-Fock and Density Functional Theory calculations are performed to test for thermodynamic feasibility, leading to several postulated stoichiometries corresponding to hazardous thermal decomposition. © 1998 Elsevier Science Ltd. All rights reserved Keywords: Computational chemistry; property estimation; process safety; thermal runaway; statistical thermo-
dynamics; quantum mechanics
1. Introduction Whenever a new product is developed in the pharmaceutical industry, the issues of thermal stability of the reactive system and product arise. Many thermal runways occurring in this industry are due to unstable intermediate or product molecules that trigger highly exothermic reactions. As an illustration of the importance of this issue, the Health and Safety Executive reported 135 exothermic runaway incidents in the United Kingdom during the period 1986 to 1991, mainly in batch reactors (Etchells, 1996). In this respect, the situation studied in this paper is of primary concern to the pharmaceutical industry, being the primary source of thermal runaways throughout the entire chemical industry (Barton and Nolan, 1991). Identification of the stoichiometries of the reactions that drive thermal instability is a major first step in understanding the safety issues related to a reactive system (Melhem et al., 1995). Despite the importance of this step, the high energy release that is typical of these reactions hinders the experimental techniques that are currently employed to study such reactive
* Author to whom correspondence should be addressed.
systems. In particular, the heat fluxes generated by such reactions are difficult to contain and measure with current calorimetry techniques. Hence, experiments are limited to very small reaction masses, and only very limited measurements (i.e., temperature profiles only) and experimental designs can be contemplated. This leads to incomplete understanding of unstable reactive systems, and a safety analysis that is at best approximate. Several different tools have been developed to investigate the thermal stability of a reactive mixture. One of the most widely used in the industry is the CHETAH software (Frurip, 1992; Davis et al., 1985) which allows a priori prediction of the hazards associated with species of interest in a reactive system from knowledge of their molecular structures. The software is based on recognition techniques with experimental hazard data and estimated thermodynamic data using molecular group contribution methods. However, CHETAH is only a preliminary screening tool to detect sufficiently unstable molecules, and cannot detect all the instabilities that can occur in a reactive chemical system. A more complete study requires experimental work for both the detection and observation of instabilities. Experiments are typically carried out in microcalorimeters where the initial reactive mixture is heated until a source of exothermy is 735
736
C. BRUNETON et al.
detected and screened. We propose to use such experiments as the starting point for a more detailed study of the observed instabilities, in particular an attempt to elucidate the stoichiometries of the reactions that lead to thermal instabilities. Our approach is based on a screening test that selects plausible stoichiometries from a candidate set of reactions by quantifying their thermodynamic feasibility. The candidate set is obtained using experimental information on the products formed and understanding of the chemistry of the reactants. Since the pharmaceutical industry usually manufactures newly invented species, the stoichiometries often involve molecules whose thermodynamics have not been studied. Nowadays, statistical thermodynamics coupled with computational quantum chemistry based on molecular theory offers the possibility to estimate accurately the thermodynamic parameters of an extremely large range of molecules. In this respect, we use Hartree-Fock and Density Functional Theory calculations to estimate the missing thermodynamic properties of the unfamiliar species involved in postulated stoichiometries. The manufacture of tetrahydrofurfuryl benzenesulfonate (TFB) is an example of an actual runaway due to thermal instabilities that occurred in a pharmaceutical production process. The runaway led to an explosion that damaged part of the plant where the batch reactor was located. Experimental investigations using microcalorimeters such as the Reactive System Screening Tool (RSST, Fauske and Associates Inc.) indicate that the explosion was due to a runaway driven by reactions related to the thermal decomposition of the product TFB. We propose to elucidate the stoichiometries of the reactions responsible for the runaway using our screening procedure.
that the reactions cause within the mixture (equation (1)). A,H Cp
ATod = - -
(1)
There is no general agreement over a value for an adiabatic temperature rise under which a reactive system is not dangerous. According to Grewer (1989), a ATad less than 50 K is seldom hazardous. F o r organic systems where Ce are usually of the magnitude of 0.1 Kcal tool- 1 K - 1, this corresponds to enthalpies of reaction of - 30 KJ/mol. The Gibbs free energy change on reaction is given by equation (2): A,G = A,H - TArS
(2)
The enthalpic term A,H can be viewed as the sum of an ideal mixture contribution, a contribution including mixture effects such as solvent interactions, and the pressure contribution. The ideal mixture contribution includes the enthalpy of reaction calculated in the ideal gas phase, A,H idg,and the enthalpies of vaporization or sublimation evaluated for the N condensed species. The term A,H mi~ encompasses all the mixing and solvent effects. This yields:
a,//= (ad-/~ -
~ N v,~J-/, )
+ A,/-/'~' + A,H~'~
(3)
where the species are ordered such that the first N exist in a condensed phase at the temperature of interest. According to the standard definition of the enthalpy of formation, A,H ~ag is given by (4): A,H iag= ~, viAfHi
(4)
i
2. Thermodynamic feasibility of postulated stoichiometries Thermally unstable systems usually exhibit thermal decomposition and/or polymerization reactions that drive thermal runaways. These reactions are usually complex with stoichiometries that are difficult to elucidate. Theoretical studies of the underlying chemistry complemented by analytical experiments on the resulting reaction mixture are required to derive preliminary information that yields a set of postulated stoichiometries. We perform enthalpy and Gibbs free energy calculations to select plausible stoichiometries from this set which are both thermodynamically feasible and sufficiently exothermic to be hazardous. A stoichiometry is thermodynamically feasible if the corresponding Gibbs free energy change on reaction (A,G) is negative. This is a necessary condition that a stoichiometry which represents the reactions of a runaway, must satisfy. The calculation of the enthalpy change on reaction (A,H) for a stoichiometry allows us to estimate the adiabatic temperature rise
The entropic contribution to the Gibbs free energy of reaction is calculated in a similar way according to equation (5): A,S = A,S id° -
~
vi
+ A,S"~ + A,Sp'~
ie[ 1 ... N]
(5) with
a, S'd' = E ~,s~~g
(6)
i
s~aa is the absolute entropy of the species i in its ideal gas phase, known as the ideal absolute entropy. The pressure effects on the enthalpy and entropy of reaction are insignificant compared to the other contributions. Thus it is reasonable to neglect them. We also note that the term for phase transitions cancels out in calculation of the ideal mixture contribution to the Gibbs free energy change on reaction. Therefore, the Gibbs free energy change on reaction can be defined as the sum of a contribution representing the Gibbs
737
Elucidation of hazardous reaction stoichiometries free energy change considering every species to be in the ideal gas phase A,G ~ag, and a contribution which encompasses the mixing and solvent effects A , G ~ :
entropy s~dg are given by: H~dg
The solvent and mixture interaction contributions Ad-/"~ and A,S"~ are usually small when compared to respectively the enthalpy and entropy of reaction. It is only when the solvent has particularly strong affinities with the solutes, such as strong polar interactions or formation of hydrogen bonds, that these terms can be significant. Molecules manufactured in the pharmaceutical industry are often unfamiliar. Thermodynamic properties must be estimated since no information is available in the literature concerning such molecules. We use computational quantum chemistry coupled with statistical thermodynamics to perform these estimates. In a classical sense, the enthalpy of formation of a species AyH is defined as the increment in enthalpy associated with the reaction of formation of the species in the ideal gas phase from its constituent elements in their stable states at the temperature and pressure of interest. These increments cannot be calculated directly by classical thermodynamics, so they are chosen as reference enthalpies at 298.15 K where they are called the standard enthalpy of formation. Statistical thermodynamics defines a reference state for enthalpies corresponding to a system of atoms that do not interact with each other. This allows us to compute the enthalpy of formation H~~g of species i in the ideal gas state. The enthalpy of formation of a molecule m according to the classical definition is then obtained from (8), where the HI~g have also been computed for every i constituent element involved in a reaction of formation:
j=*
ridg + R T + T~s,
i
Z
v,A~bn,
(8)
is[ 1 ... N]
The i constituent elements are ordered such that the first N exist in a condensed phase at the temperature of interest. For example, if we consider the reaction of formation for methanol from its constituent elements at 298.15 K: 1 Co,~vh~t~ + 2H20as + i 02oo~ ~ CH~OHg,,
Thus the classical enthalpy of formation for methanol can be obtained given knowledge of statistical mechanical enthalpies of formation and the heat of sublimation of carbon as follows: {,~Hidg . 1 LlidO ) afHcrl~On = Hg~1~on - ~ rl~ -v 2 "" o~ + H ~ g
+ A,.bHc. ......
(9)
Using the development of McQuarrie (1976), the enthalpy of formation HIag and the ideal absolute
(10)
and
[
n- 6 Or, I T s~d~ = stidg ..... + sldg ,o,, + R 3E j= t Lexp Ovj/T - 1 i
ln(1 - exp - ® v J T ) + Rln coel] (11) where n is the number of fundamental vibrational frequencies. The translational contributions are respectively for the internal energy and for the ideal absolute entropy uidg trans i = ~ R T
and ,dg Rln[(2~Mky'~2/3ves/R 1 s,..... = [_\ h 2 /I N]
For linear molecules, U toil ld* = R T
and
while for non-linear molecules, uidg rot~ = ~ R T
and AsH(m) "--"HI~ - Z viHiag +
rTidg
~trans~ -~ V r o t t
3n-6 {®vj ®vj + R T Z \ 2 T + T(exp ®v,/T - 1)J
(7)
A,G = A,G idg + A,G "ix
ITidg
~1/2e3/2(
Ta
~1/2
Rln
®~, is the vibrational temperature function of the 3N - 6 molecular normal vibrational frequencies v~. The ®,o,'S are the rotational temperatures, c%, is the degeneracy corresponding to the ground state energy of the molecule. The electronic energy contribution U~,g is given by the ground state energy of the molecule, which can be computed using quantum chemical calculations based on the Hartree-Fock method (HF) (Hehre et al., 1986) or Density Functional Theory (DFT) (Parr, 1989). These methods are also called Ab-Initio methods and allow us to solve numerically the time independent Schr6dinger equation. The vfs correspond to the eigenvalues of the Hessian energy matrix as a function of the molecular degrees of freedom computed after an Ab-Initio calculation. Figure 1 depicts the thermodynamic feasibility test. A stoichiometry is retained if the corresponding Gibbs free energy is negative. Quantum computational chemistry is used to estimate the enthalpies of
C. BRUNETONet al.
738
i
stoichiometry
Availablethermodynamicdata , on the species involved I
[ No available thermodynamic data
I
on the species involved
J CalculationOfArl-lidg, Arsidg Calculationof ArGidg
ArGidg= ArHidg. TArSidg of the mixing effects
AGrnix
i Stoichiometry ,o occ ,e ll
,o>0
Stoichtometry rejected I
Fig. 1. Screening procedure for stoichiometry determination.
formation and the ideal absolute entropies of the unfamiliar species involved in a postulated stoichiometry. It allows us to calculate the enthalpy and entropy of reaction and the Gibbs free energy of reaction A,G ~agconsidering the species in the ideal gas phase. The solvent effects ArGm~x can be calculated using two general approaches: classical ensemble treatments and quantum mechanical continuum models. In the classical treatments are included Monte Carlo statistical methods (Jorgensen and Ravimoha, 1985) and free energy perturbations (Bash et al., 1987; Jorgensen and Buckner, 1987). The other approach is based on the Onsager reaction field model (Onsager, 1936) where the solvent is simply a continuous unstructured dielectric with a given dielectric constant embedded in some kind of cavity. The shape of the cavity depends on the solute geometries. This has been studied by many authors (Miertus and Tomasi, 1982); (Bonaccorsi et al., 1991); (Wong et al., 1991); (Foresman et al., 1996). Except for solvents having particularly strong interactions with the solute (e.g., interactions due to hydrogen bonds), the magnitude of the solvent contribution to the Gibbs free energy change on reaction is less than 10 Kcal/mol. In this respect the sign of A,G is known when IA,GI > e with: e = I0 + ~ IviJlel[ [Kcal/mol] i
(12)
where el are the uncertainties in Kcal/mol in the Gibbs free energy of formation of species i, that are calculated from the uncertainties in the enthalpies of formation and ideal absolute entropies of species i. In other words, the sign of A,G can be obtained from a calculation of ArG~d9and the corresponding uncertainty. 3. Use of computational quantum chemistry
Computational quantum chemistry is based on quantum molecular theory where the motion and distribution of electrons is described in term of probability distributions or molecular orbitals. These are solutions of the Schrrdinger equation which is solved numerically by applying the BornOppenheimer approximation (which ignores motion of the nuclei, roughly two thousand times heavier than the electrons). Computational quantum chemistry refers to two main numerical techniques known as Hartree-Fock (HF) and Density Functional Theory (DFT). The HF method is the simplest level of calculation since some exchange correlation effects between electrons are ignored. This is due to the mean-field approximation (Hehre et al., 1986), which reduces solution of the Schrrdinger equation for the N electrons of a molecule to the solution of N one-electron SchriSdinger equations. Although electron exchange
Elucidation of hazardous reaction stoichiometries effects due to the principle of Pauli repulsion are taken into account, the application of the mean-field approximation ignores correlations between electrons which are crucial to the prediction of accurate energies (Simka et al., 1996). D F T is based on a theory (Hohenberg and Kohn, 1964; Sham, 1965) which states that the ground state energy of a system of electrons is a function of the electron charge density. D F T attempts to describe more accurately electron exchange and correlation effects using a number of functionals which include terms for exchange and correlation energy. These functionals depend on the electron density and possibly the electron density gradient. Local D F T corresponds to the first level of approximation using the Local Density Approximation (LDA) (Parr, 1989). The LDA is based on the assumption that the exchange-correlation energy per particle at any location is the same as that obtained for an homogeneous electron gas with the same density. In other words, an inhomogeneous electron density distribution is considered to be locally homogeneous. Gradient corrections are used in Non local D F T (Perdrew, 1995) to represent the variation of the electron density. Several authors have proposed a number of functionals depending on the electron density and its gradients (Perdrew, 1995; Becke, 1988; Lee et al., 1988). H F and D F T solve the Schrrdinger equation using the SelfConsistent-Field (SCF) method which is an iterative process (Hehre et al., 1986). Molecular orbitals are restricted to linear combinations of a set of known one-electron Gaussian basis functions. Both H F and D F T coupled with statistical thermodynamics allow us to compute thermodynamic properties for any kind of molecule. In this respect, they
739
offer an advantage compared to molecular group contribution methods which cannot be applied to the unconventional molecules commonly used in the pharmaceutical industry. As a matter of fact, many groups are not implemented in these latter methods, such as groups containing sulphur, boron, phosphorus or silicon atoms, and organometallic groups. In addition, charged molecules such as carbocations or carbanions cannot be treated. The accuracy of the estimation performed using H F and D F T with statistical thermodynamics depends on the method and the basis set used. As an illustration, the standard enthalpy of formation AfH~298.1SK ) of silane (Sill4) is estimated using Non Local D F T based on the Becke functional (Becke, 1988) for the exchange energy and the LPY functional (Lee et al., 1988) for the correlation energy. We also use the Becke3LPY method (Becke, 1993) where the exchange-correlation energy is represented by a functional using the Becke and LPY functionals and the Hartree-Fock exchange term. The calculations are performed with different sized basis sets, providing different degrees of accuracy. We also propose hybrid calculations where the ground state energy and the vibrational frequencies are calculated using D F T and H F respectively. Experience shows that H F vibrational frequencies scaled by a factor of 0.89 can lead to better results than D F T vibrational frequencies (Hehre et al., 1986). The calculations presented in Table 1 are carried out using the software GAUSSIAN 94 (Frisch et al., 1994) on a H P 9000/735 workstation. The use of a more sophisticated exchange-correlation energy functional and an increase in the size of the basis set lead to better estimations. We demonstrate in Table 2 the efficiency of the Ab-Initio methods coupled with statistical
Table 1. Standard heats of formation of Sill4 (Kcal/mol) (experimental value: 8.2 + 1 Kcal/mol) and corresponding CPU times Method
ASH(298.15 K)
Absolute error
CPU time
15.3
7.1
40 min
14.6
6.4
69 rain
13.5
5.3
84 min
11.1
2.9
760 min
Energy: BLPY/6-31G** a Frequencies: HF x 0.89/6-31G** Energy: B3LPY/6-31G** Frequencies: HF x 0.89/6-31G** Energy: B3LPY/6-311G** Frequencies: HF x 0.89/6-31G** Energy: B3LPY/6-311G(3df,2pd) Frequencies: HF x 0.89/6-31G** a the 6-31G** is also written 6-31G(d,p).
Table 2. Standard heats of formation (Kcal/mol) Molecule
Ab-initio methods
Group contribution methods
Experiment
C,,H6 CH3SCN
30.8 38.1
23.6 (Constantinou) 29.2 (Joback)
26.3 38.3
740
C. BRUNETONet al.
thermodynamics by estimating the standard enthalpy of formation for 1,3-butadiene and methyl thiocyanate using BLPY/6-31G(d,p) DFT. Since group contribution methods such as the Joback (Joback, 1990), the Benson (Benson, 1968) and the Constantinou (Constantinou and Gani, 1994) methods are applicable for both molecules, the results obtained are compared. The results for 1,3-butadiene contain a similar amount of error. However computational quantum chemistry leads to a more accurate estimation for methyl thioeyanate than the Joback group contribution method, which is the only group contribution method that is applicable for this molecule. As a matter of fact, the properties of methyl thiocyanate result from the divalent sulphur atom bonded to a cyanate group. This association of atoms is not implemented in any of the above group contribution methods. Further, distinguishing between isomers presents no problem for Ab-Initio methods, unlike most group contribution methods. Since molecular group contribution methods are based on correlations obtained from large databanks of experimental values of thermodynamic properties for common molecules, they are efficient when applied to familiar molecules containing groups of atoms which are included in the correlations (e.g., 1,3-butadiene). However, these methods become highly questionable when the molecule is sufficiently unconventional.
4. Example The manufacture of tetrahydrofurfuryl benzenesulfonate (TFB) offers a typical example of thermal runaway due to thermal instabilities which actually occurred in the pharmaceutical industry. The batch process involves an esterification type reaction in an organic solvent as the main reaction. The reaction is run at 20-30°C under atmospheric pressure (Fig. 2). A reactor where this reaction was carried out exploded spreading its contents throughout the plant and destroying all the surrounding equipment. This runaway has been extensively studied in the laboratory. Experiments were carried out with the product molecule TFB as the only reactant in the nonpolar organic solvent under 3 bars in a RSST. Temperature rates of 14,000°C/min and pressure rates of 1,500 bar/min were recorded revealing that the runaway was due to reactions starting from the product TFB as the only reactant. Consequently the runaway is due to secondary reactions following the main reaction of the process. This is a typical class of thermal
o
II
O Fig. 2. M a i n reaction.
A:G = AyH - TAyS
+
HO--CH2
(13)
The entropy of formation of a molecule m is calculated according to equation (14): AzS(,,) = s(~ ia - ~ v~sl
(14)
i
v~and s~are respectively the stoichiometric coefficients and the absolute entropies for the i constituent elements involved in the reaction of formation, in their stable state at the temperature and pressure of interest. Reference molecules are simple molecules containing atoms and bonds that also appear in TFB and benzenesulfonic acid. No experimental results on
O "
II
--S--CI
runaway that is initiated by a loss of temperature control of the main reactions, also called primary reactions. The reactor temperature increases and hits the onset temperature of undesired secondary reactions characterized by high heats of reaction, large activation energies and high reaction rates. These secondary reactions dictate the magnitude and time scales of heat release during a runaway. The gas phase given off by the RSST has been analyzed by gas chromatography-mass spectrometry giving some information on the products formed. However no attempts have been made to analyze the liquid phase which appears as a brown-black residue. According to the reactant, the products measured in the gas phase, and previous work done on thermal decomposition of molecules similar to TFB (Baumgartner, 1953; Sites, 1965), a series of stoichiometries for hazardous reactions are postulated (Fig. 3). The screening test to identify thermodynamically feasible and hazardous stoichiometries is applied to this set. Most of the molecules involved are common, so their enthalpies of formation, ideal absolute entropies, and enthalpies of vaporization or sublimation are available in the literature. However, no thermodynamic data exist for TFB. Therefore its ideal gas enthalpy of formation and ideal absolute entropy is estimated using Ab-Initio calculations. They are also estimated for benzenesulfonic acid since few thermodynamic data are available for this component in the literature, and the figures reported contain large experimental error. Attention is focused to choose an efficient AbInitio method providing good accuracy. In this respect we tested different methods by carrying out a series of reference calculations of standard Gibbs free energy of formation AIG on reference molecules. AIG is defined according to equation (13):
• --S--O--CH2 II
O
+ HC!
Elucidation of hazardous reaction stoichiometries
o
/-~ (a) o" ---"
.-,s.-o0
o,, ,,,-s-oH,, +
741
~o
0 0 ~
CO + H2 +
+ SO 2 + CO + H 2 f 8 O
O
It
• --S-OHII ~
l ~ H + SO2 +-'~ 02
O Fig. 3. Plausible stoiehiometries of hazardous reactions.
Table 3. CPU time for the reference calculations, HP 9000/735 Molecule Benzene Ethane Dimethyl ether Sulfur dioxide Methane thiol
Energy & Frequency B88-LPY 6--31G** 5 hours 45 min 1 hour 3 hours 21 min 77 min
molecules containing a sulfate group were found in the literature. Therefore the accuracy of a method on sulphur compounds can only be tested with respect to reference molecules containing divalent and tetravalent sulphur groups. Three different methods are used with Non Local D F T and hybrid D F T - H F , where the molecular ground state energy is calculated with Non Local D F T and the vibrational normal frequencies are calculated with HF. Differences between the calculated standard Gibbs free energies of formation and the corresponding experimental values are compared. A method is retained if these differences do not exceed 10 Kcal/mol for every reference molecule. Experimental standard Gibbs free energy of formation are extracted from (Daubert and Danner, 1989). These values are given with the experimental error below ___1 Kcal/mol. B88-LPY D F T (Becke, 1988; Lee et al., 1988) with 6-31G** basis set and
Energy B88-LPY 6-31G**
Frequency HF 6-31G*
56 min 18 min 27 min 13 min 31 rain
hybrid methods mixing B88-LPY D F T for the ground state energy calculation and H F for the vibrational frequencies calculation are tested. B88-LPY D F T gives acceptable results and is used to estimate the enthalpy of formation and ideal absolute entropy of benzenesulfonic acid. However the method appears to be rather expensive in terms of CPU time, especially for the normal frequency computation. The vibrational contribution to the Gibbs free energy of formation of a species is a rather insensitive function of the vibrational frequencies. In this respect, we suggest the use of a less sophisticated method for the computation of frequencies such as H F with a smaller basis set such as 6-31G*. This method leads to shorter overall C P U times when compared with Non Local D F T with 6-31G** (Table 3). Table 4 shows that acceptable results are obtained with B88-LPY D F T and hybrid D F T - H F
742
C. BRUNETON et al.
Table 4. Error in calculated standard Gibbs free energy of formation (Kcal/mol) Molecule
Methane Ethane Dimethyl ether Benzene Sulfur dioxide Methane thiol
Energy & Frequency Energy B88-LPY 6-31G** B88-LPY 6-31G** - 0.7 - 4.3 9.4 - 5.5 7.2 - 7.1
- 2.3 - 7.4 5.4 - 10.7 6.7 - 9.4
Table 5. Estimated thermodynamic data for TFB and benzenesulfonicacid Molecule
TFB benzenesulfonic acid
Frequency HF 6-31G*
AIH Kcalmol- 1
s~dg calmol- 1 K - 1
- 122.5
131.9
- 87.9
96.6
when multiplying the vibrational frequencies by a scaling factor of 0.89. Hybrid D F T - H F gives less accurate results. However, it has been noticed that harmonic vibrational frequencies computed at the H F level are systematically larger than the corresponding experimental values. This deficiency is corrected by scaling the vibrational frequencies by a factor of 0.89 as recommended by (Hehre et al., 1986). This leads to results similar to the B88-LPY D F T results. This hybrid method is retained for the calculation of the standard enthalpy of formation and ideal absolute entropy of TFB because of the shorter CPU times observed and accuracy comparable with B88-LPY DFT. The results obtained are displayed in Table 5. Molecular group contribution methods are not applicable for TFB and benzenesulfonic acid since groups containing hexavalent sulphur are not implemented (Benson, 1968; Joback, 1990; Constantinou and Gani, 1994). Bearing the chemical species involved and the nonpolar organic solvent used, it is reasonable to assume that solvent effects do not contribute significantly to the Gibbs free energies of reaction. Therefore, according to the above discussion, the contribution of A,Gm~x is at most 10 Kcal/mol, and the sign of A,G can be obtained from a calculation of A,G~ag with the associated uncertainty as described above. In this respect, all data are now available to perform the thermodynamic feasibility test on the postulated stoichiometries depicted in Fig. 3. We calculate the Gibbs free energy at 298.15 K and the corresponding uncertainty e for every postulated stoichiometry. This requires knowledge of the uncertainties in the enthalpies of formation and the ideal absolute entropies for each species involved in the stoichiometry. These values are
Energy Frequency B88-LPY 6-31G** HF ×0.89 6--31G* - 0.4 - 4.0 9.1 - 5.3 7.0 - 7.4
Table 6. Calculated Gibbs free energy of reaction with the corresponding uncertainties for plausible stoichiometries in Kcal/mol Stoichiometry (a) (b) (c) (d) (e) (f) (g) (h)
A,G(298.15 K) - 0.6 - 14.5 22.2 - 24.9 - 20.3 - 22.8 - 9.1 22
829s.15 32 14 15 13 12.5 24 12 21
provided in the thermodynamic database of (Daubert and Danner, 1989), which gives values of Gibbs free energies of formation for the common molecules involved in the set of plausible stoichiometries with relative experimental errors that are lower than 5%. The uncertainty in the Gibbs free energy of formation of TFB and benzenesulfonic acid is estimated to be 10Kcal/mol according to the criterion chosen to select the Ab-Initio method from the reference calculations. This criterion probably significantly overestimates the value of e, and is somewhat conservative. For example, the calculation for stoichiometry (e) leads to A,G~g --- - 20.3 Kcal/mol. The uncertainty in the Gibbs free energy of formation are respectively for dihydrofuran, tetrahydrofuran and hydrogen, 1.4Kcal/mol, 1.1 Kcal/mol and 0.02Kcal/mol. According to equation(12), this leads to e = 12.5 Kcal/mol. We conclude that the sign of A,G is negative, and stoichiometry (e) is definitely thermodynamically feasible. Table 6 gives the calculated values of A,G and e for each postulated stoichiometry. Reaction (c) and (h) are found infeasible. However (c) can be combined with favourable reactions such as (b), (d) and (e) to obtain a favourable stoichiometry. In this respect we postulate global stoichiometries from the remaining reactions paying special attention to obtain as products all the species observed experimentally in the gas phase. This leads to the stoichiometries (a), (i) and (j) (Figs 3 and 4). The stoichiometry (i) corresponds to 2(b) + 2(c) + 2(d) + (e).
Elucidation of hazardous reaction stoichiometries 2
=•O + H2
2 TFB ~
0)
~ 2 co + N + ~ /
2 ~H + 2 SO2 + 2 CO + H2 + ] ~
+/'~
/© ~'o"
743
fonic acid. The elucidation of these reactions and estimation of their heats of reaction is a complex problem since the structure of the copolymers is unknown. The problem could be treated by considering models of multicomponent copolymerization (Odian, 1991) with calculation of heats of copolymerization using computational quantum chemistry.
Fig. 4. Stoichiometries (i) and (j). 5. Conclusions
Table 7. Calculated Gibbs free energy of reaction with uncertainty at 298.15 K (Kcal/mol) Stoichiometry
A,G(298.15 K)
a i j
e(298.15K)
- 0.6 - 56 - 54
32 17 38
Table 8. Calculated enthalpy of reaction at 298.15K (Kcal/mol) Stoichiometry a i j
A,H(298.15 K) 0,3 -
A~apH(TeB)) -17 - 12 + 2A,aff/(rvn)
(A~bH~®so~r) -
The stoichiometry (j) corresponds to 2(f)+ (g). Table 7 depicts the calculated Gibbs free energies of reaction at 298.15 K with the corresponding uncertainties for the three global stoichiometries. Stoichiometries (i) and (j) are fully accepted. However no conclusion can be drawn for (a) without a calculation of the solvent effect on the corresponding Gibbs free energy of reaction. Table 8 presents the enthalpies of reaction of the proposed global stoichiometries neglecting the solvent effect. (i) is the only stoichiometry corresponding to an exothermic reaction. The magnitude of heat capacities of organic molecules is usually 0.1 Kcal/mol-lK -1. In this respect an enthalpy of reaction of - 17 Kcal/mol leads to an adiabatic temperature rise on the order of magnitude of 100 K which is sufficient to provide heat for further exothermic reactions leading to the temperature rise experimentally observed. Consequently our study provides an explanation for the onset part of the runaway. Further propositions of stoichiometries are required to elucidate the entire phenomenon characterized by the much higher adiabatic temperature rise observed in microcalorimetry experiments. This may include multicomponent polymerizations occurring in the liquid phase and leading to the heavy liquid residue observed experimentally. In this respect we can easily think of cationic multicomponent copolymerizations with the unsaturated liquid species involved in the stoichiometries (a), (i) and (j) catalysed by benzenesul-
We have presented a screening test to determine stoichiometries responsible for thermal instabilities, based on their thermodynamic feasibility. The test has been proposed to elucidate reactions due to thermal instability of product and intermediate species in the pharmaceutical industry. In this respect we have shown a new application of computational quantum chemistry to determine the thermodynamic properties of unfamiliar species typically involved in thermal runaway accidents. The technique has been applied to an example of a real thermal runaway due to the thermal decomposition of tetrahydrofurfuryl benzenesulfonate. We have been able to propose stoichiometries for the decomposition reactions consistent with analytical experimental observations of gaseous products. These stoichiometries and the calculated heats of reaction can potentially be used in a dynamic model of the runaway for validation purposes when the corresponding kinetic parameters are obtained. In this respect, transition state theory coupled with computational quantum chemistry can be used. Our screening test finds its limits in the elucidation of polymerization reactions due to thermal instabilities. Models of multicomponent copolymers in conjunction with computational quantum chemistry are needed to address this problem. Finally, although this paper is illustrated by an example drawn from the pharmaceutical industry, it should be noted that the impact of a better understanding of the chemical reactions responsible for thermal runaways is not limited to this industry. In particular, competition in commodity chemical manufacturing is encouraging manufacturers to operate their processes closer to limits and constraints. The traditional origin of such constraints is somewhat arbitrary, typically based on recommendations from studies at the bench and/or pilot scale. However, a rigorous study of the mechanisms of the reactions responsible for runaways opens the possibility for a deeper, chemically based understanding of the constraints within which a process may be operated. Those companies that can best exploit this knowledge in the operation of their processes will potentially gain a major competitive advantage. Notation
c Cp
velocity of light (ms- 1) heat capacity (Jmol- 1 K- 1)
744
ArG A,H
H~dg AsH mvapn As,bH h k M N R Si
A,S O"
T
ATa~ Ui V V vj ~rot 8o 60el
C. BRUNETON et al. exp(1) uncertainty in the Gibbs free energy of reaction (Jmol - 1) uncertainty in the Gibbs free energy of formation for a species i (Jmol-1) Gibbs free energy of reaction (Jmol- 1) enthalpy of reaction (Jmol- l) enthalpy of formation (statistical thermodynamic reference) (Jmol- 1) enthalpy of formation (classical thermodynamic reference) (Jmol-1) enthalpy of vaporization (Jmol- 1) enthalpy of sublimation (Jmol- 1) Planck constant (= 6.63.10-34 Js) Boltzmann constant (= 1.38.10-23 JK-x) molar mass (kg) number of molecules ideal gas constant (= 8.314 Jm- 3 mol - 1 K - 1) absolute entropy (Jmol- 1 K - 1) entropy of reaction (Jmol- 1 K - 1) rotational symmetry number temperature (K) adiabatic temperature rise(K) molar internal energy (Jmol- 1) volume (m 3) stoichiometric coefficient normal vibrational frequency (m- 1) rotational temperature (K) vibrational temperature (K) degeneracy of the ground state energy of a molecule.
Superscript idg mix pres
Ideal gas mixing contribution pressure contribution
Acknowledgements We are greatly indebted to Arup K. Chakraborty for the education in quantum and statistical mechanics he has given us, and to Elf Aquitaine for financial support of this project.
References Barton, J.A. and Nolan, P.F. (1991) Incidents in chemical industry due to thermal runaway of chemical reactions. In Conference on Chemical Hazards, London Press Centre. Bash, P.A., Singh, U.C., Langridge, R. and Kollman, P.A. (1987). Free energy calculations by computer simulation. Science 236, 564. Baumgartner, G.J. (1953) The Pyrolysis of some tetrahydrofurfuryl esters, PhD thesis, University of Notre Dame, USA. Becke, A.D. (1988) Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. 38(6A), 3098. Becke, A.D. (1993) Density-functional thermochemistry. III. the role of exact exchange. J. Chem. Phys. 98, 5648. Benson, S.W. (1968) Thermodynamical Kinetics, chapter 2. Wiley, New York.
Bonaccorsi, R., Cammi, R. and Tomasi, J. (1991) On the ab initio geometry optimization of molecular solutes. Journal of Computational Chemistry 12, 301. Constantinou, L. and Gani, R. (1994) Group-contribution method for estimating properties of pure compounds. AIChE Journal 40, 1697. Daubert, T.E. and Danner, R.P. (1989) Physical and Thermodynamic Properties of Pure Chemicals Data Compilation, AIChE. Davis, C.A., Kipnis, I.M., Chase, M.W. and Trewwek, D. N. (1985) The thermochemical and hazard data of chemicals, estimation using ASTM CHETAH program. In ACS Symposium Series No. 274, Chemical Process Hazard Review. Etchells, J. (1996) Runaway. The Chemical Engineer 17, September 12th. Foresman, J.B., Keith, T.A., Wiberg, K.B., Snoonian, J. and Frish, M.J. (1996) Solvent effects. 5. influence of cavity shape, truncation of electrostatics, and electron correlation on ab initio reaction field calculations. J. Phys. Chem. 100, 16098. Frisch, M.J., Frisch, LE. and Foreman, J.B. (1994) Gaussian 94, Gaussian, Inc., Carnegie Office Park, Pittsburgh, PA 15106, USA. Frurip, D.J. (1992) Using the ASTM CHETAH program in chemical process hazard evaluation. Plant~Operations Progress 11, 224. Grewer, T. (1989) Thermal Hazards of Chemical Reactions, volume 4 of Industrial Safety Series. Elsevier, Amsterdam, Hehre, W.J., Random, L., Schleyer, P. and Pople, J.A. (1986) Ab initio molecular theory, Wiley, New York. Hohenherg, P. and Kohn, W, (1964). Inhomogeneous electron gas. Phys. Rev. 136(3B), 864. Joback, K.G. (1990) Computer-aided design of molecules with desired properties, Laboratory for Intelligent Systems in Process Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA. Jorgensen, W.L and Buckner, J.K. (1987) Use of statistical perturbation theory for computing solvent effects on molecular conformation. Butane in water. J. Phys. Chem. 91, 6083. Jorgensen, W.L. and Ravimoha, C. (1985) Monte Carlo simulation of differences in free energies of hydration. J. Chem. Phys. 83, 3050. Lee, C., Yang, W. and Parr, R.G. (1988) Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B37, 785. McQuarrie. (1976) Statistical Mechanics. Harper and Row, New York. Melhem, G.A., Fisher, H.G. and Shaw, D.A. (1995) An advanced method for the estimation of reaction kinetics, scaleup, and pressure relief design, Process Safety Progress 14(1). Miertus, S. and Tomasi, J. (1982) Approximate evaluations of the electrostatic free energy and internal energy changes in solution processes. Chem. Phys. 65, 239. Odian, G.G. (1991) Principles of Polymerization, chapter 6. Wiley, New York. Onsager, L. (1936). Electric moments of molecules in liquids. J. Am. Chem. Soc. 58, 1486. Parr, R.G. (1989) Density-functional theory of atoms and molecules. Oxford University Press, New York.
Elucidation of hazardous reaction stoichiometries Perdrew, J.P. (1995) Nonlocal density functionals for exchange and correlation: theory and applications. In Density functional theory of molecules, clusters, and solids, Kluwer Academic Publishers, Dordrecht, Boston, pp, 47-66. Sham, L.J. (1965) Self-consistent equations including exchange and correlation effects. Phys. Rev. 140(4A), 1133. Simka, H., Hierlemann, M., Utz, M. and Jensen, K.F. (1996) Computational chemistry predictions of kinetics and
745
major reaction pathways for germane gas-phase reactions, Journal of the Electrochemical Society 143. Sites, R.D. (1965) A mechanism for and the identification of dimeric olefins arising from the pyrolysis of (-)-menthyl benzenesulfonate, PhD thesis, Brown University, USA. Wong, M.W., Frisch, M.J. and Wiberg, K.B. (1991) Solvent effects. 1. The mediation of electrostatic effects by solvents. J. Am. Chem. Soc. 113, 4776.
Pergamon PII: S0098-1354(97)00256-1
© 1998ElsevierScienceLtd Printed in Great Britain. All rights reserved 0098-1354/98 $19.00+ 0.00
Thermal stability analysis via elucidation of hazardous reaction stoichiometries C. Bruneton", C. Hoff b and P.I. Barton"'* a Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, U.S.A. b SANOFI CHIMIE, 30390 Aramon, France
Abstract Many accidents that occur in the chemical processing industries such as thermal runaways are due to unstable reactive systems. A good understanding of these systems requires knowledge of the stoichiometries of the reactions that lead to instabilities, which are typically decompositions and/or polymerizations. A screening test to elucidate these stoichiometries is presented. It eliminates implausible stoichiometries from a candidate set due to their thermodynamic infeasibility. We use computational chemistry coupled with statistical thermodynamics to estimate the thermodynamic properties of the unfamiliar species typically involved in such reactions. The real life runaway due to product instability in a process to manufacture tetrahydrofurfuryl benzenesulfonate is treated as a case study. Hartree-Fock and Density Functional Theory calculations are performed to test for thermodynamic feasibility, leading to several postulated stoichiometries corresponding to hazardous thermal decomposition. © 1998 Elsevier Science Ltd. All rights reserved Keywords: Computational chemistry; property estimation; process safety; thermal runaway; statistical thermo-
dynamics; quantum mechanics
1. Introduction Whenever a new product is developed in the pharmaceutical industry, the issues of thermal stability of the reactive system and product arise. Many thermal runways occurring in this industry are due to unstable intermediate or product molecules that trigger highly exothermic reactions. As an illustration of the importance of this issue, the Health and Safety Executive reported 135 exothermic runaway incidents in the United Kingdom during the period 1986 to 1991, mainly in batch reactors (Etchells, 1996). In this respect, the situation studied in this paper is of primary concern to the pharmaceutical industry, being the primary source of thermal runaways throughout the entire chemical industry (Barton and Nolan, 1991). Identification of the stoichiometries of the reactions that drive thermal instability is a major first step in understanding the safety issues related to a reactive system (Melhem et al., 1995). Despite the importance of this step, the high energy release that is typical of these reactions hinders the experimental techniques that are currently employed to study such reactive
* Author to whom correspondence should be addressed.
systems. In particular, the heat fluxes generated by such reactions are difficult to contain and measure with current calorimetry techniques. Hence, experiments are limited to very small reaction masses, and only very limited measurements (i.e., temperature profiles only) and experimental designs can be contemplated. This leads to incomplete understanding of unstable reactive systems, and a safety analysis that is at best approximate. Several different tools have been developed to investigate the thermal stability of a reactive mixture. One of the most widely used in the industry is the CHETAH software (Frurip, 1992; Davis et al., 1985) which allows a priori prediction of the hazards associated with species of interest in a reactive system from knowledge of their molecular structures. The software is based on recognition techniques with experimental hazard data and estimated thermodynamic data using molecular group contribution methods. However, CHETAH is only a preliminary screening tool to detect sufficiently unstable molecules, and cannot detect all the instabilities that can occur in a reactive chemical system. A more complete study requires experimental work for both the detection and observation of instabilities. Experiments are typically carried out in microcalorimeters where the initial reactive mixture is heated until a source of exothermy is 735
736
C. BRUNETON et al.
detected and screened. We propose to use such experiments as the starting point for a more detailed study of the observed instabilities, in particular an attempt to elucidate the stoichiometries of the reactions that lead to thermal instabilities. Our approach is based on a screening test that selects plausible stoichiometries from a candidate set of reactions by quantifying their thermodynamic feasibility. The candidate set is obtained using experimental information on the products formed and understanding of the chemistry of the reactants. Since the pharmaceutical industry usually manufactures newly invented species, the stoichiometries often involve molecules whose thermodynamics have not been studied. Nowadays, statistical thermodynamics coupled with computational quantum chemistry based on molecular theory offers the possibility to estimate accurately the thermodynamic parameters of an extremely large range of molecules. In this respect, we use Hartree-Fock and Density Functional Theory calculations to estimate the missing thermodynamic properties of the unfamiliar species involved in postulated stoichiometries. The manufacture of tetrahydrofurfuryl benzenesulfonate (TFB) is an example of an actual runaway due to thermal instabilities that occurred in a pharmaceutical production process. The runaway led to an explosion that damaged part of the plant where the batch reactor was located. Experimental investigations using microcalorimeters such as the Reactive System Screening Tool (RSST, Fauske and Associates Inc.) indicate that the explosion was due to a runaway driven by reactions related to the thermal decomposition of the product TFB. We propose to elucidate the stoichiometries of the reactions responsible for the runaway using our screening procedure.
that the reactions cause within the mixture (equation (1)). A,H Cp
ATod = - -
(1)
There is no general agreement over a value for an adiabatic temperature rise under which a reactive system is not dangerous. According to Grewer (1989), a ATad less than 50 K is seldom hazardous. F o r organic systems where Ce are usually of the magnitude of 0.1 Kcal tool- 1 K - 1, this corresponds to enthalpies of reaction of - 30 KJ/mol. The Gibbs free energy change on reaction is given by equation (2): A,G = A,H - TArS
(2)
The enthalpic term A,H can be viewed as the sum of an ideal mixture contribution, a contribution including mixture effects such as solvent interactions, and the pressure contribution. The ideal mixture contribution includes the enthalpy of reaction calculated in the ideal gas phase, A,H idg,and the enthalpies of vaporization or sublimation evaluated for the N condensed species. The term A,H mi~ encompasses all the mixing and solvent effects. This yields:
a,//= (ad-/~ -
~ N v,~J-/, )
+ A,/-/'~' + A,H~'~
(3)
where the species are ordered such that the first N exist in a condensed phase at the temperature of interest. According to the standard definition of the enthalpy of formation, A,H ~ag is given by (4): A,H iag= ~, viAfHi
(4)
i
2. Thermodynamic feasibility of postulated stoichiometries Thermally unstable systems usually exhibit thermal decomposition and/or polymerization reactions that drive thermal runaways. These reactions are usually complex with stoichiometries that are difficult to elucidate. Theoretical studies of the underlying chemistry complemented by analytical experiments on the resulting reaction mixture are required to derive preliminary information that yields a set of postulated stoichiometries. We perform enthalpy and Gibbs free energy calculations to select plausible stoichiometries from this set which are both thermodynamically feasible and sufficiently exothermic to be hazardous. A stoichiometry is thermodynamically feasible if the corresponding Gibbs free energy change on reaction (A,G) is negative. This is a necessary condition that a stoichiometry which represents the reactions of a runaway, must satisfy. The calculation of the enthalpy change on reaction (A,H) for a stoichiometry allows us to estimate the adiabatic temperature rise
The entropic contribution to the Gibbs free energy of reaction is calculated in a similar way according to equation (5): A,S = A,S id° -
~
vi
+ A,S"~ + A,Sp'~
ie[ 1 ... N]
(5) with
a, S'd' = E ~,s~~g
(6)
i
s~aa is the absolute entropy of the species i in its ideal gas phase, known as the ideal absolute entropy. The pressure effects on the enthalpy and entropy of reaction are insignificant compared to the other contributions. Thus it is reasonable to neglect them. We also note that the term for phase transitions cancels out in calculation of the ideal mixture contribution to the Gibbs free energy change on reaction. Therefore, the Gibbs free energy change on reaction can be defined as the sum of a contribution representing the Gibbs
737
Elucidation of hazardous reaction stoichiometries free energy change considering every species to be in the ideal gas phase A,G ~ag, and a contribution which encompasses the mixing and solvent effects A , G ~ :
entropy s~dg are given by: H~dg
The solvent and mixture interaction contributions Ad-/"~ and A,S"~ are usually small when compared to respectively the enthalpy and entropy of reaction. It is only when the solvent has particularly strong affinities with the solutes, such as strong polar interactions or formation of hydrogen bonds, that these terms can be significant. Molecules manufactured in the pharmaceutical industry are often unfamiliar. Thermodynamic properties must be estimated since no information is available in the literature concerning such molecules. We use computational quantum chemistry coupled with statistical thermodynamics to perform these estimates. In a classical sense, the enthalpy of formation of a species AyH is defined as the increment in enthalpy associated with the reaction of formation of the species in the ideal gas phase from its constituent elements in their stable states at the temperature and pressure of interest. These increments cannot be calculated directly by classical thermodynamics, so they are chosen as reference enthalpies at 298.15 K where they are called the standard enthalpy of formation. Statistical thermodynamics defines a reference state for enthalpies corresponding to a system of atoms that do not interact with each other. This allows us to compute the enthalpy of formation H~~g of species i in the ideal gas state. The enthalpy of formation of a molecule m according to the classical definition is then obtained from (8), where the HI~g have also been computed for every i constituent element involved in a reaction of formation:
j=*
ridg + R T + T~s,
i
Z
v,A~bn,
(8)
is[ 1 ... N]
The i constituent elements are ordered such that the first N exist in a condensed phase at the temperature of interest. For example, if we consider the reaction of formation for methanol from its constituent elements at 298.15 K: 1 Co,~vh~t~ + 2H20as + i 02oo~ ~ CH~OHg,,
Thus the classical enthalpy of formation for methanol can be obtained given knowledge of statistical mechanical enthalpies of formation and the heat of sublimation of carbon as follows: {,~Hidg . 1 LlidO ) afHcrl~On = Hg~1~on - ~ rl~ -v 2 "" o~ + H ~ g
+ A,.bHc. ......
(9)
Using the development of McQuarrie (1976), the enthalpy of formation HIag and the ideal absolute
(10)
and
[
n- 6 Or, I T s~d~ = stidg ..... + sldg ,o,, + R 3E j= t Lexp Ovj/T - 1 i
ln(1 - exp - ® v J T ) + Rln coel] (11) where n is the number of fundamental vibrational frequencies. The translational contributions are respectively for the internal energy and for the ideal absolute entropy uidg trans i = ~ R T
and ,dg Rln[(2~Mky'~2/3ves/R 1 s,..... = [_\ h 2 /I N]
For linear molecules, U toil ld* = R T
and
while for non-linear molecules, uidg rot~ = ~ R T
and AsH(m) "--"HI~ - Z viHiag +
rTidg
~trans~ -~ V r o t t
3n-6 {®vj ®vj + R T Z \ 2 T + T(exp ®v,/T - 1)J
(7)
A,G = A,G idg + A,G "ix
ITidg
~1/2e3/2(
Ta
~1/2
Rln
®~, is the vibrational temperature function of the 3N - 6 molecular normal vibrational frequencies v~. The ®,o,'S are the rotational temperatures, c%, is the degeneracy corresponding to the ground state energy of the molecule. The electronic energy contribution U~,g is given by the ground state energy of the molecule, which can be computed using quantum chemical calculations based on the Hartree-Fock method (HF) (Hehre et al., 1986) or Density Functional Theory (DFT) (Parr, 1989). These methods are also called Ab-Initio methods and allow us to solve numerically the time independent Schr6dinger equation. The vfs correspond to the eigenvalues of the Hessian energy matrix as a function of the molecular degrees of freedom computed after an Ab-Initio calculation. Figure 1 depicts the thermodynamic feasibility test. A stoichiometry is retained if the corresponding Gibbs free energy is negative. Quantum computational chemistry is used to estimate the enthalpies of
C. BRUNETONet al.
738
i
stoichiometry
Availablethermodynamicdata , on the species involved I
[ No available thermodynamic data
I
on the species involved
J CalculationOfArl-lidg, Arsidg Calculationof ArGidg
ArGidg= ArHidg. TArSidg of the mixing effects
AGrnix
i Stoichiometry ,o occ ,e ll
,o>0
Stoichtometry rejected I
Fig. 1. Screening procedure for stoichiometry determination.
formation and the ideal absolute entropies of the unfamiliar species involved in a postulated stoichiometry. It allows us to calculate the enthalpy and entropy of reaction and the Gibbs free energy of reaction A,G ~agconsidering the species in the ideal gas phase. The solvent effects ArGm~x can be calculated using two general approaches: classical ensemble treatments and quantum mechanical continuum models. In the classical treatments are included Monte Carlo statistical methods (Jorgensen and Ravimoha, 1985) and free energy perturbations (Bash et al., 1987; Jorgensen and Buckner, 1987). The other approach is based on the Onsager reaction field model (Onsager, 1936) where the solvent is simply a continuous unstructured dielectric with a given dielectric constant embedded in some kind of cavity. The shape of the cavity depends on the solute geometries. This has been studied by many authors (Miertus and Tomasi, 1982); (Bonaccorsi et al., 1991); (Wong et al., 1991); (Foresman et al., 1996). Except for solvents having particularly strong interactions with the solute (e.g., interactions due to hydrogen bonds), the magnitude of the solvent contribution to the Gibbs free energy change on reaction is less than 10 Kcal/mol. In this respect the sign of A,G is known when IA,GI > e with: e = I0 + ~ IviJlel[ [Kcal/mol] i
(12)
where el are the uncertainties in Kcal/mol in the Gibbs free energy of formation of species i, that are calculated from the uncertainties in the enthalpies of formation and ideal absolute entropies of species i. In other words, the sign of A,G can be obtained from a calculation of ArG~d9and the corresponding uncertainty. 3. Use of computational quantum chemistry
Computational quantum chemistry is based on quantum molecular theory where the motion and distribution of electrons is described in term of probability distributions or molecular orbitals. These are solutions of the Schrrdinger equation which is solved numerically by applying the BornOppenheimer approximation (which ignores motion of the nuclei, roughly two thousand times heavier than the electrons). Computational quantum chemistry refers to two main numerical techniques known as Hartree-Fock (HF) and Density Functional Theory (DFT). The HF method is the simplest level of calculation since some exchange correlation effects between electrons are ignored. This is due to the mean-field approximation (Hehre et al., 1986), which reduces solution of the Schrrdinger equation for the N electrons of a molecule to the solution of N one-electron SchriSdinger equations. Although electron exchange
Elucidation of hazardous reaction stoichiometries effects due to the principle of Pauli repulsion are taken into account, the application of the mean-field approximation ignores correlations between electrons which are crucial to the prediction of accurate energies (Simka et al., 1996). D F T is based on a theory (Hohenberg and Kohn, 1964; Sham, 1965) which states that the ground state energy of a system of electrons is a function of the electron charge density. D F T attempts to describe more accurately electron exchange and correlation effects using a number of functionals which include terms for exchange and correlation energy. These functionals depend on the electron density and possibly the electron density gradient. Local D F T corresponds to the first level of approximation using the Local Density Approximation (LDA) (Parr, 1989). The LDA is based on the assumption that the exchange-correlation energy per particle at any location is the same as that obtained for an homogeneous electron gas with the same density. In other words, an inhomogeneous electron density distribution is considered to be locally homogeneous. Gradient corrections are used in Non local D F T (Perdrew, 1995) to represent the variation of the electron density. Several authors have proposed a number of functionals depending on the electron density and its gradients (Perdrew, 1995; Becke, 1988; Lee et al., 1988). H F and D F T solve the Schrrdinger equation using the SelfConsistent-Field (SCF) method which is an iterative process (Hehre et al., 1986). Molecular orbitals are restricted to linear combinations of a set of known one-electron Gaussian basis functions. Both H F and D F T coupled with statistical thermodynamics allow us to compute thermodynamic properties for any kind of molecule. In this respect, they
739
offer an advantage compared to molecular group contribution methods which cannot be applied to the unconventional molecules commonly used in the pharmaceutical industry. As a matter of fact, many groups are not implemented in these latter methods, such as groups containing sulphur, boron, phosphorus or silicon atoms, and organometallic groups. In addition, charged molecules such as carbocations or carbanions cannot be treated. The accuracy of the estimation performed using H F and D F T with statistical thermodynamics depends on the method and the basis set used. As an illustration, the standard enthalpy of formation AfH~298.1SK ) of silane (Sill4) is estimated using Non Local D F T based on the Becke functional (Becke, 1988) for the exchange energy and the LPY functional (Lee et al., 1988) for the correlation energy. We also use the Becke3LPY method (Becke, 1993) where the exchange-correlation energy is represented by a functional using the Becke and LPY functionals and the Hartree-Fock exchange term. The calculations are performed with different sized basis sets, providing different degrees of accuracy. We also propose hybrid calculations where the ground state energy and the vibrational frequencies are calculated using D F T and H F respectively. Experience shows that H F vibrational frequencies scaled by a factor of 0.89 can lead to better results than D F T vibrational frequencies (Hehre et al., 1986). The calculations presented in Table 1 are carried out using the software GAUSSIAN 94 (Frisch et al., 1994) on a H P 9000/735 workstation. The use of a more sophisticated exchange-correlation energy functional and an increase in the size of the basis set lead to better estimations. We demonstrate in Table 2 the efficiency of the Ab-Initio methods coupled with statistical
Table 1. Standard heats of formation of Sill4 (Kcal/mol) (experimental value: 8.2 + 1 Kcal/mol) and corresponding CPU times Method
ASH(298.15 K)
Absolute error
CPU time
15.3
7.1
40 min
14.6
6.4
69 rain
13.5
5.3
84 min
11.1
2.9
760 min
Energy: BLPY/6-31G** a Frequencies: HF x 0.89/6-31G** Energy: B3LPY/6-31G** Frequencies: HF x 0.89/6-31G** Energy: B3LPY/6-311G** Frequencies: HF x 0.89/6-31G** Energy: B3LPY/6-311G(3df,2pd) Frequencies: HF x 0.89/6-31G** a the 6-31G** is also written 6-31G(d,p).
Table 2. Standard heats of formation (Kcal/mol) Molecule
Ab-initio methods
Group contribution methods
Experiment
C,,H6 CH3SCN
30.8 38.1
23.6 (Constantinou) 29.2 (Joback)
26.3 38.3
740
C. BRUNETONet al.
thermodynamics by estimating the standard enthalpy of formation for 1,3-butadiene and methyl thiocyanate using BLPY/6-31G(d,p) DFT. Since group contribution methods such as the Joback (Joback, 1990), the Benson (Benson, 1968) and the Constantinou (Constantinou and Gani, 1994) methods are applicable for both molecules, the results obtained are compared. The results for 1,3-butadiene contain a similar amount of error. However computational quantum chemistry leads to a more accurate estimation for methyl thioeyanate than the Joback group contribution method, which is the only group contribution method that is applicable for this molecule. As a matter of fact, the properties of methyl thiocyanate result from the divalent sulphur atom bonded to a cyanate group. This association of atoms is not implemented in any of the above group contribution methods. Further, distinguishing between isomers presents no problem for Ab-Initio methods, unlike most group contribution methods. Since molecular group contribution methods are based on correlations obtained from large databanks of experimental values of thermodynamic properties for common molecules, they are efficient when applied to familiar molecules containing groups of atoms which are included in the correlations (e.g., 1,3-butadiene). However, these methods become highly questionable when the molecule is sufficiently unconventional.
4. Example The manufacture of tetrahydrofurfuryl benzenesulfonate (TFB) offers a typical example of thermal runaway due to thermal instabilities which actually occurred in the pharmaceutical industry. The batch process involves an esterification type reaction in an organic solvent as the main reaction. The reaction is run at 20-30°C under atmospheric pressure (Fig. 2). A reactor where this reaction was carried out exploded spreading its contents throughout the plant and destroying all the surrounding equipment. This runaway has been extensively studied in the laboratory. Experiments were carried out with the product molecule TFB as the only reactant in the nonpolar organic solvent under 3 bars in a RSST. Temperature rates of 14,000°C/min and pressure rates of 1,500 bar/min were recorded revealing that the runaway was due to reactions starting from the product TFB as the only reactant. Consequently the runaway is due to secondary reactions following the main reaction of the process. This is a typical class of thermal
o
II
O Fig. 2. M a i n reaction.
A:G = AyH - TAyS
+
HO--CH2
(13)
The entropy of formation of a molecule m is calculated according to equation (14): AzS(,,) = s(~ ia - ~ v~sl
(14)
i
v~and s~are respectively the stoichiometric coefficients and the absolute entropies for the i constituent elements involved in the reaction of formation, in their stable state at the temperature and pressure of interest. Reference molecules are simple molecules containing atoms and bonds that also appear in TFB and benzenesulfonic acid. No experimental results on
O "
II
--S--CI
runaway that is initiated by a loss of temperature control of the main reactions, also called primary reactions. The reactor temperature increases and hits the onset temperature of undesired secondary reactions characterized by high heats of reaction, large activation energies and high reaction rates. These secondary reactions dictate the magnitude and time scales of heat release during a runaway. The gas phase given off by the RSST has been analyzed by gas chromatography-mass spectrometry giving some information on the products formed. However no attempts have been made to analyze the liquid phase which appears as a brown-black residue. According to the reactant, the products measured in the gas phase, and previous work done on thermal decomposition of molecules similar to TFB (Baumgartner, 1953; Sites, 1965), a series of stoichiometries for hazardous reactions are postulated (Fig. 3). The screening test to identify thermodynamically feasible and hazardous stoichiometries is applied to this set. Most of the molecules involved are common, so their enthalpies of formation, ideal absolute entropies, and enthalpies of vaporization or sublimation are available in the literature. However, no thermodynamic data exist for TFB. Therefore its ideal gas enthalpy of formation and ideal absolute entropy is estimated using Ab-Initio calculations. They are also estimated for benzenesulfonic acid since few thermodynamic data are available for this component in the literature, and the figures reported contain large experimental error. Attention is focused to choose an efficient AbInitio method providing good accuracy. In this respect we tested different methods by carrying out a series of reference calculations of standard Gibbs free energy of formation AIG on reference molecules. AIG is defined according to equation (13):
• --S--O--CH2 II
O
+ HC!
Elucidation of hazardous reaction stoichiometries
o
/-~ (a) o" ---"
.-,s.-o0
o,, ,,,-s-oH,, +
741
~o
0 0 ~
CO + H2 +
+ SO 2 + CO + H 2 f 8 O
O
It
• --S-OHII ~
l ~ H + SO2 +-'~ 02
O Fig. 3. Plausible stoiehiometries of hazardous reactions.
Table 3. CPU time for the reference calculations, HP 9000/735 Molecule Benzene Ethane Dimethyl ether Sulfur dioxide Methane thiol
Energy & Frequency B88-LPY 6--31G** 5 hours 45 min 1 hour 3 hours 21 min 77 min
molecules containing a sulfate group were found in the literature. Therefore the accuracy of a method on sulphur compounds can only be tested with respect to reference molecules containing divalent and tetravalent sulphur groups. Three different methods are used with Non Local D F T and hybrid D F T - H F , where the molecular ground state energy is calculated with Non Local D F T and the vibrational normal frequencies are calculated with HF. Differences between the calculated standard Gibbs free energies of formation and the corresponding experimental values are compared. A method is retained if these differences do not exceed 10 Kcal/mol for every reference molecule. Experimental standard Gibbs free energy of formation are extracted from (Daubert and Danner, 1989). These values are given with the experimental error below ___1 Kcal/mol. B88-LPY D F T (Becke, 1988; Lee et al., 1988) with 6-31G** basis set and
Energy B88-LPY 6-31G**
Frequency HF 6-31G*
56 min 18 min 27 min 13 min 31 rain
hybrid methods mixing B88-LPY D F T for the ground state energy calculation and H F for the vibrational frequencies calculation are tested. B88-LPY D F T gives acceptable results and is used to estimate the enthalpy of formation and ideal absolute entropy of benzenesulfonic acid. However the method appears to be rather expensive in terms of CPU time, especially for the normal frequency computation. The vibrational contribution to the Gibbs free energy of formation of a species is a rather insensitive function of the vibrational frequencies. In this respect, we suggest the use of a less sophisticated method for the computation of frequencies such as H F with a smaller basis set such as 6-31G*. This method leads to shorter overall C P U times when compared with Non Local D F T with 6-31G** (Table 3). Table 4 shows that acceptable results are obtained with B88-LPY D F T and hybrid D F T - H F
742
C. BRUNETON et al.
Table 4. Error in calculated standard Gibbs free energy of formation (Kcal/mol) Molecule
Methane Ethane Dimethyl ether Benzene Sulfur dioxide Methane thiol
Energy & Frequency Energy B88-LPY 6-31G** B88-LPY 6-31G** - 0.7 - 4.3 9.4 - 5.5 7.2 - 7.1
- 2.3 - 7.4 5.4 - 10.7 6.7 - 9.4
Table 5. Estimated thermodynamic data for TFB and benzenesulfonicacid Molecule
TFB benzenesulfonic acid
Frequency HF 6-31G*
AIH Kcalmol- 1
s~dg calmol- 1 K - 1
- 122.5
131.9
- 87.9
96.6
when multiplying the vibrational frequencies by a scaling factor of 0.89. Hybrid D F T - H F gives less accurate results. However, it has been noticed that harmonic vibrational frequencies computed at the H F level are systematically larger than the corresponding experimental values. This deficiency is corrected by scaling the vibrational frequencies by a factor of 0.89 as recommended by (Hehre et al., 1986). This leads to results similar to the B88-LPY D F T results. This hybrid method is retained for the calculation of the standard enthalpy of formation and ideal absolute entropy of TFB because of the shorter CPU times observed and accuracy comparable with B88-LPY DFT. The results obtained are displayed in Table 5. Molecular group contribution methods are not applicable for TFB and benzenesulfonic acid since groups containing hexavalent sulphur are not implemented (Benson, 1968; Joback, 1990; Constantinou and Gani, 1994). Bearing the chemical species involved and the nonpolar organic solvent used, it is reasonable to assume that solvent effects do not contribute significantly to the Gibbs free energies of reaction. Therefore, according to the above discussion, the contribution of A,Gm~x is at most 10 Kcal/mol, and the sign of A,G can be obtained from a calculation of A,G~ag with the associated uncertainty as described above. In this respect, all data are now available to perform the thermodynamic feasibility test on the postulated stoichiometries depicted in Fig. 3. We calculate the Gibbs free energy at 298.15 K and the corresponding uncertainty e for every postulated stoichiometry. This requires knowledge of the uncertainties in the enthalpies of formation and the ideal absolute entropies for each species involved in the stoichiometry. These values are
Energy Frequency B88-LPY 6-31G** HF ×0.89 6--31G* - 0.4 - 4.0 9.1 - 5.3 7.0 - 7.4
Table 6. Calculated Gibbs free energy of reaction with the corresponding uncertainties for plausible stoichiometries in Kcal/mol Stoichiometry (a) (b) (c) (d) (e) (f) (g) (h)
A,G(298.15 K) - 0.6 - 14.5 22.2 - 24.9 - 20.3 - 22.8 - 9.1 22
829s.15 32 14 15 13 12.5 24 12 21
provided in the thermodynamic database of (Daubert and Danner, 1989), which gives values of Gibbs free energies of formation for the common molecules involved in the set of plausible stoichiometries with relative experimental errors that are lower than 5%. The uncertainty in the Gibbs free energy of formation of TFB and benzenesulfonic acid is estimated to be 10Kcal/mol according to the criterion chosen to select the Ab-Initio method from the reference calculations. This criterion probably significantly overestimates the value of e, and is somewhat conservative. For example, the calculation for stoichiometry (e) leads to A,G~g --- - 20.3 Kcal/mol. The uncertainty in the Gibbs free energy of formation are respectively for dihydrofuran, tetrahydrofuran and hydrogen, 1.4Kcal/mol, 1.1 Kcal/mol and 0.02Kcal/mol. According to equation(12), this leads to e = 12.5 Kcal/mol. We conclude that the sign of A,G is negative, and stoichiometry (e) is definitely thermodynamically feasible. Table 6 gives the calculated values of A,G and e for each postulated stoichiometry. Reaction (c) and (h) are found infeasible. However (c) can be combined with favourable reactions such as (b), (d) and (e) to obtain a favourable stoichiometry. In this respect we postulate global stoichiometries from the remaining reactions paying special attention to obtain as products all the species observed experimentally in the gas phase. This leads to the stoichiometries (a), (i) and (j) (Figs 3 and 4). The stoichiometry (i) corresponds to 2(b) + 2(c) + 2(d) + (e).
Elucidation of hazardous reaction stoichiometries 2
=•O + H2
2 TFB ~
0)
~ 2 co + N + ~ /
2 ~H + 2 SO2 + 2 CO + H2 + ] ~
+/'~
/© ~'o"
743
fonic acid. The elucidation of these reactions and estimation of their heats of reaction is a complex problem since the structure of the copolymers is unknown. The problem could be treated by considering models of multicomponent copolymerization (Odian, 1991) with calculation of heats of copolymerization using computational quantum chemistry.
Fig. 4. Stoichiometries (i) and (j). 5. Conclusions
Table 7. Calculated Gibbs free energy of reaction with uncertainty at 298.15 K (Kcal/mol) Stoichiometry
A,G(298.15 K)
a i j
e(298.15K)
- 0.6 - 56 - 54
32 17 38
Table 8. Calculated enthalpy of reaction at 298.15K (Kcal/mol) Stoichiometry a i j
A,H(298.15 K) 0,3 -
A~apH(TeB)) -17 - 12 + 2A,aff/(rvn)
(A~bH~®so~r) -
The stoichiometry (j) corresponds to 2(f)+ (g). Table 7 depicts the calculated Gibbs free energies of reaction at 298.15 K with the corresponding uncertainties for the three global stoichiometries. Stoichiometries (i) and (j) are fully accepted. However no conclusion can be drawn for (a) without a calculation of the solvent effect on the corresponding Gibbs free energy of reaction. Table 8 presents the enthalpies of reaction of the proposed global stoichiometries neglecting the solvent effect. (i) is the only stoichiometry corresponding to an exothermic reaction. The magnitude of heat capacities of organic molecules is usually 0.1 Kcal/mol-lK -1. In this respect an enthalpy of reaction of - 17 Kcal/mol leads to an adiabatic temperature rise on the order of magnitude of 100 K which is sufficient to provide heat for further exothermic reactions leading to the temperature rise experimentally observed. Consequently our study provides an explanation for the onset part of the runaway. Further propositions of stoichiometries are required to elucidate the entire phenomenon characterized by the much higher adiabatic temperature rise observed in microcalorimetry experiments. This may include multicomponent polymerizations occurring in the liquid phase and leading to the heavy liquid residue observed experimentally. In this respect we can easily think of cationic multicomponent copolymerizations with the unsaturated liquid species involved in the stoichiometries (a), (i) and (j) catalysed by benzenesul-
We have presented a screening test to determine stoichiometries responsible for thermal instabilities, based on their thermodynamic feasibility. The test has been proposed to elucidate reactions due to thermal instability of product and intermediate species in the pharmaceutical industry. In this respect we have shown a new application of computational quantum chemistry to determine the thermodynamic properties of unfamiliar species typically involved in thermal runaway accidents. The technique has been applied to an example of a real thermal runaway due to the thermal decomposition of tetrahydrofurfuryl benzenesulfonate. We have been able to propose stoichiometries for the decomposition reactions consistent with analytical experimental observations of gaseous products. These stoichiometries and the calculated heats of reaction can potentially be used in a dynamic model of the runaway for validation purposes when the corresponding kinetic parameters are obtained. In this respect, transition state theory coupled with computational quantum chemistry can be used. Our screening test finds its limits in the elucidation of polymerization reactions due to thermal instabilities. Models of multicomponent copolymers in conjunction with computational quantum chemistry are needed to address this problem. Finally, although this paper is illustrated by an example drawn from the pharmaceutical industry, it should be noted that the impact of a better understanding of the chemical reactions responsible for thermal runaways is not limited to this industry. In particular, competition in commodity chemical manufacturing is encouraging manufacturers to operate their processes closer to limits and constraints. The traditional origin of such constraints is somewhat arbitrary, typically based on recommendations from studies at the bench and/or pilot scale. However, a rigorous study of the mechanisms of the reactions responsible for runaways opens the possibility for a deeper, chemically based understanding of the constraints within which a process may be operated. Those companies that can best exploit this knowledge in the operation of their processes will potentially gain a major competitive advantage. Notation
c Cp
velocity of light (ms- 1) heat capacity (Jmol- 1 K- 1)
744
ArG A,H
H~dg AsH mvapn As,bH h k M N R Si
A,S O"
T
ATa~ Ui V V vj ~rot 8o 60el
C. BRUNETON et al. exp(1) uncertainty in the Gibbs free energy of reaction (Jmol - 1) uncertainty in the Gibbs free energy of formation for a species i (Jmol-1) Gibbs free energy of reaction (Jmol- 1) enthalpy of reaction (Jmol- l) enthalpy of formation (statistical thermodynamic reference) (Jmol- 1) enthalpy of formation (classical thermodynamic reference) (Jmol-1) enthalpy of vaporization (Jmol- 1) enthalpy of sublimation (Jmol- 1) Planck constant (= 6.63.10-34 Js) Boltzmann constant (= 1.38.10-23 JK-x) molar mass (kg) number of molecules ideal gas constant (= 8.314 Jm- 3 mol - 1 K - 1) absolute entropy (Jmol- 1 K - 1) entropy of reaction (Jmol- 1 K - 1) rotational symmetry number temperature (K) adiabatic temperature rise(K) molar internal energy (Jmol- 1) volume (m 3) stoichiometric coefficient normal vibrational frequency (m- 1) rotational temperature (K) vibrational temperature (K) degeneracy of the ground state energy of a molecule.
Superscript idg mix pres
Ideal gas mixing contribution pressure contribution
Acknowledgements We are greatly indebted to Arup K. Chakraborty for the education in quantum and statistical mechanics he has given us, and to Elf Aquitaine for financial support of this project.
References Barton, J.A. and Nolan, P.F. (1991) Incidents in chemical industry due to thermal runaway of chemical reactions. In Conference on Chemical Hazards, London Press Centre. Bash, P.A., Singh, U.C., Langridge, R. and Kollman, P.A. (1987). Free energy calculations by computer simulation. Science 236, 564. Baumgartner, G.J. (1953) The Pyrolysis of some tetrahydrofurfuryl esters, PhD thesis, University of Notre Dame, USA. Becke, A.D. (1988) Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. 38(6A), 3098. Becke, A.D. (1993) Density-functional thermochemistry. III. the role of exact exchange. J. Chem. Phys. 98, 5648. Benson, S.W. (1968) Thermodynamical Kinetics, chapter 2. Wiley, New York.
Bonaccorsi, R., Cammi, R. and Tomasi, J. (1991) On the ab initio geometry optimization of molecular solutes. Journal of Computational Chemistry 12, 301. Constantinou, L. and Gani, R. (1994) Group-contribution method for estimating properties of pure compounds. AIChE Journal 40, 1697. Daubert, T.E. and Danner, R.P. (1989) Physical and Thermodynamic Properties of Pure Chemicals Data Compilation, AIChE. Davis, C.A., Kipnis, I.M., Chase, M.W. and Trewwek, D. N. (1985) The thermochemical and hazard data of chemicals, estimation using ASTM CHETAH program. In ACS Symposium Series No. 274, Chemical Process Hazard Review. Etchells, J. (1996) Runaway. The Chemical Engineer 17, September 12th. Foresman, J.B., Keith, T.A., Wiberg, K.B., Snoonian, J. and Frish, M.J. (1996) Solvent effects. 5. influence of cavity shape, truncation of electrostatics, and electron correlation on ab initio reaction field calculations. J. Phys. Chem. 100, 16098. Frisch, M.J., Frisch, LE. and Foreman, J.B. (1994) Gaussian 94, Gaussian, Inc., Carnegie Office Park, Pittsburgh, PA 15106, USA. Frurip, D.J. (1992) Using the ASTM CHETAH program in chemical process hazard evaluation. Plant~Operations Progress 11, 224. Grewer, T. (1989) Thermal Hazards of Chemical Reactions, volume 4 of Industrial Safety Series. Elsevier, Amsterdam, Hehre, W.J., Random, L., Schleyer, P. and Pople, J.A. (1986) Ab initio molecular theory, Wiley, New York. Hohenherg, P. and Kohn, W, (1964). Inhomogeneous electron gas. Phys. Rev. 136(3B), 864. Joback, K.G. (1990) Computer-aided design of molecules with desired properties, Laboratory for Intelligent Systems in Process Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA. Jorgensen, W.L and Buckner, J.K. (1987) Use of statistical perturbation theory for computing solvent effects on molecular conformation. Butane in water. J. Phys. Chem. 91, 6083. Jorgensen, W.L. and Ravimoha, C. (1985) Monte Carlo simulation of differences in free energies of hydration. J. Chem. Phys. 83, 3050. Lee, C., Yang, W. and Parr, R.G. (1988) Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B37, 785. McQuarrie. (1976) Statistical Mechanics. Harper and Row, New York. Melhem, G.A., Fisher, H.G. and Shaw, D.A. (1995) An advanced method for the estimation of reaction kinetics, scaleup, and pressure relief design, Process Safety Progress 14(1). Miertus, S. and Tomasi, J. (1982) Approximate evaluations of the electrostatic free energy and internal energy changes in solution processes. Chem. Phys. 65, 239. Odian, G.G. (1991) Principles of Polymerization, chapter 6. Wiley, New York. Onsager, L. (1936). Electric moments of molecules in liquids. J. Am. Chem. Soc. 58, 1486. Parr, R.G. (1989) Density-functional theory of atoms and molecules. Oxford University Press, New York.
Elucidation of hazardous reaction stoichiometries Perdrew, J.P. (1995) Nonlocal density functionals for exchange and correlation: theory and applications. In Density functional theory of molecules, clusters, and solids, Kluwer Academic Publishers, Dordrecht, Boston, pp, 47-66. Sham, L.J. (1965) Self-consistent equations including exchange and correlation effects. Phys. Rev. 140(4A), 1133. Simka, H., Hierlemann, M., Utz, M. and Jensen, K.F. (1996) Computational chemistry predictions of kinetics and
745
major reaction pathways for germane gas-phase reactions, Journal of the Electrochemical Society 143. Sites, R.D. (1965) A mechanism for and the identification of dimeric olefins arising from the pyrolysis of (-)-menthyl benzenesulfonate, PhD thesis, Brown University, USA. Wong, M.W., Frisch, M.J. and Wiberg, K.B. (1991) Solvent effects. 1. The mediation of electrostatic effects by solvents. J. Am. Chem. Soc. 113, 4776.