22
Ind. Eng. C h e m . Res. 1990, 29, 22-29
Quantitative Use of Model Compound Information: Monte Carlo Simulation of the Reactions of Complex Macromolecules John B. McDermott,+Cristian Libanati, Concetta LaMarca, and Michael T. Klein* Center for Catalytic Science and Technology, Department of Chemical Engineering, University of Delaware, Newark. Delabare 19716
A Monte Carlo simulation of complex depolymerization reactions was developed. Transition probabilities derived through analysis of simple prototype reaction sequences were used to test reactive sites in a polymer. This stochastic approach allowed facile accounting of the often confounding effects of unequal reactivity, nearest-neighbor substituents, and oligomer molecular weight distributions. T h e simulations matched all available deterministic solutions and further showed that the confounding effects noted above can influence the effectiveness of process strategies derived from model compound information. Chemists and chemical engineers often address problems involving complex reactions such as heavy hydrocarbon upgrading or polymer synthesis. Information about such systems is generally sought from two complementary perspectives. At a practical extreme, experiments with the actual reactants of interest provide the most relevant information, often in global terms, such as product boiling points, gel fraction, or solubility classes. This is because basic kinetics concepts with a precise meaning in simpler systems (e.g., conversion or rate of reaction) are not directly or easily applied. These highly relevant experiments often yield very little fundamental information. This motivates model compound experiments, a fundamental extreme aimed at the resolution of the intrinsic chemistry of the relevant real system. However, during the reaction of the complex system, physical factors extrinsic to the reactive moiety, e.g., substituent effects or transport of the polymeric molecule, can disguise its intrinsic chemistry. The relevance of model compound information is therefore frequently unclear. Chemical modeling is being developed as a formalism through which intrinsic chemistry can be determined and ultimately used in the prediction of the reaction of complex systems. It comprises three components. The first is a data base of intrinsic reaction pathways and kinetics, discerned from model compound and computational chemistry experiments (Neurock et al., 1989). The second is composed of engineering science models that account for the differences (Le., molecular-weight-dependent diffusivities, phase behavior, substituent effects) between model systems and the real system they are intended to mimic. The third is the computational reaction engineering framework that enables the quantitative a priori simulation of complex systems by using the intrinsic data base and engineering science models. For macromolecular substrates with complicated and conversion-dependent distributions of reactive moieties and physical characteristics, such as lignins. coals, or abphaltenes, deterministic approaches can be cumbersome and difficult t o implement. In these instances, a Monte Carlo simulation, which can incorporate stochastic kinetics, is often more convenient because it renders the physical structure of the reactant explicit. This enables not only the use of the model compound intrinsic reactivity information but also allows facile accounting of the physical differences between model and real systems. *Corresponding author. 'Present address: General Electric Company, Corporate Research and Development, Schenectady, NY 13301.
Table I. Pathways and Kinetics for Thermal and Catalytic Reactions of VGE and GGE log E*/kcal thermal reactionsafc A/s9 mol-' (1) VGE acetophenone + guaiacol 5.03 23.9 12) GGE acetovanillone + guaiacol 5.03 20.7 (3) VGE water + vinyl-VGE 1.87 12.5 (4) GGE water + vinyl-GGE 1.87 12.7 ( 5 ) vinyl-VGE guaiacol + char 4.03 15.4 (6) vinyl-GGE guaiacol + vinylguaiacol 4.03 15.4 (char) reactions over catalystajd -(log k ) / s - ' (7) VGE acetophenone + guaiacol 3.35 (8) GGE guaiacol + char 2.86 3.11 (9) VGE VGE keto ether (10) GGE GGE keto ether 3.11 (11) VGE keto ether acetophenone + guaiacol 2.97 (12) GGE keto ether acetovanillone + guaiacol 2.46
------
--
See Figure 1 for chemical structures. bAt 6.67 g of catalyst and 250 "C. 'McDermott (1986) and McDermott et al. (1986). dMcDermott (1986) and Domburg et al. (1974).
Herein we use the depolymerization of a model lignin, poly(veratry1 P-guaiacyl ether), hereafter referred to as polyVGE, as a vehicle to develop and illustrate a model compound based Monte Carlo simulation of complex systems.
PolyVGE as a Model Lignin Lignin is a phenolic copolymer of sinapyl, coniferyl, and coumaryl alcohol monomers. Its mechanism of growth is through coupling of the various resonance structures of the monomeric radicals formed by enzymatic oxidative dehydrogenation (Freudenberg and Neish, 1968). This growth affords several types of inter-phenolic unit linkages, the most predominant being the @-ethersthat arise through coupling of radicals centered at the phenolic oxygen and the $-carbon of the monomers' propanoid side chains. Shown in Figure 1, polyVGE is a P-ether-linked linear polymer that captures many of the essential features of lignin reactivity. Our approach is to use kinetics information discerned from relevant polyVGE model compounds. The thermal and catalytic reaction pathways and kinetics of VGE and its associated reaction intermediates, the structures of which are shown in Figure 1,are summarized in Table I. Pyrolysis serves as a base case. VGE thermal reactions are to the product pairs guaiacol and acetophenone and water and char (vinylguaiacol polymerizes to char very rapidly); a typical selectivity is S = guaiacol/water = 1/40. U'ithin polvVGE, the analogous reaction of VGE to
0888-S885/90/2629-0022$02.50/0 C 1990 American Chemical Society
Ind. Eng. Chem. Res., Vol. 29, No. 1, 1990 23
L
n
PolyVGE Okc
Okc
OH
OH
bO&-& GGE Veratryl-P -guaiacyi ether
OH
Guaiacy1.p -guaiacyl ether Okc
& O w ; ;
Vinyl VGE QY
Vinyl GGE
okc
0
&%;; VGE Keto-Elher
k0; k
Acetophenone
OHb
Acetovanillone
0
GGE Keto-Ether
OH
'Ti Vinylguaiacol
Guaiacol
Figure 1. Chemical structures for polyVGE, its reaction products, and related model compounds.
guaiacol would lead to a polyVGE oligomer whose terminal monomer is hydroxyl substituted. This indicates the relevance of the reaction of GGE, which is structurally similar to VGE save the replacement of a methoxyl by a hydroxyl substituent: its reaction is faster = 1.5) and more selective to guaiacol and acetovanillone (S = 2/3). The study of VGE catalytic dehydrogenation pathways was aimed at the improvement of hydrogen utilization during lignin liquefaction. As described more fully elsewhere (McDermott et al., 1986), the strategy was aimed at prevention of water evolution and char formation. The key step was efficient catalytic dehydrogenation of VGE's aliphatic alcohol to carbonyls, whose subsequent thermal reactions would lead to evolution of lignin oxygen as carbon oxides. Thus, not only would hydrogen remain available for further hydroprocessing but also the lignin residue would be less carbon rich due to rejection of carbon with the lignin oxygen. Experiments with the model compounds showed the catalytic reaction of VGE and GGE to single-ring products to be about 10 times faster than their thermal reactions. Moreover, the undesirable thermal dehydration reaction was suppressed through dehydrogenation to keto ethers, which, in turn, fragment rapidly to single-ring products. However, our use of this information in deterministic modeling approaches was frustrated by the complexity of possible lignin reaction paths and structure. For example, during polyVGE pyrolysis, the reaction of a given VGE link will transform a nearest neighbor into a more reactive GGE link. Thus, the reactivity of the basic interunit link will be conversion dependent. Moreover, the catalytic reaction rates depend on the diffusivities of the reacting species, which in turn depend on the conversion-dependent molecular weight. These factors would render the observable lignin chemistry different from the intrinsic chemistry provided by model compound experiments. This motivated the formulation of the present Monte Carlo approach for the quantitative use of model compound information in complex reaction systems. Monte Carlo Simulation Approach Monte Carlo calculations and the mathematics of Markov chains have been used in a number of probabilistic problems and have been particularly popular in the field
of polymer science. Either or both of these mathematical techniques have been used in the investigation of such topics as polymer conformation (Kinsinger, 1970; Windwer, 1970), copolymer composition and tacticity (Price, 1970), the excluded volume problem (Windwer, 19701, and molecular weight distributions (Lowry, 1970). These techniques provide a vehicle for the use of statistical information about a system. Herein that information pertains to reaction probabilities. The principles of stochastic kinetics were summarized by McQuarrie (1967). Traditional approaches involve derivation of the stochastic master equation, whose solution provides kinetics information. Solution methods range from classical analytical methods to numerical solutions, including Monte Carlo realizations (McQuarrie, 1967; Turner, 1977; Lopez-Serrano, et al., 1980; Gillespie, 1976). Schaad (1963) provided an early example of the use of Monte Carlo calculations in the integration of rate equations. A more recent application of stochastic kinetics using a Monte Carlo computational procedure was provided by Gillespie (1976, 1977a,b). His algorithm defined a system as an entire population of reacting molecules. The state of this system was defined by the number of each type of molecule it comprised. Reactions were transitions of the system from one state to another. The time step required for a reaction was calculated through the logarithmic distribution 7 = ( l / a )In ( l / r ) , where 0 < r < 1 is a random number and a is the sum of all the reaction rate constants. The actual reaction path was selected randomly using probabilities based on the relative rates. The complementary Monte Carlo simulation technique developed herein considers fixed time intervals ( A t )passing in series; rigorous transition probabilities, described below, allow determination of events (reactions) occurring during each At. This approach permits examination of each reactive moiety for a potential reaction, independently of many other moieties and events in the system. The cumbersome list of all possible events in the complex system, inherent to Gillespie's simulation, is thus eliminated in this approach, and very disparate kinetic constants can be handled easily. Logic of t h e Simulation The logic of the present Monte Carlo simulation approach is schematically illustrated in Figure 2. A reaction trajectory of this arbitrarily linear polyVGE is described as a Markov chain, which depicts the transitions of the system, which is, in this example, a collection of nine aromatic units. Its state is determined by the state of the bonds that can join the units. Originally, at t = 0, the nine phenolic rings are joined by eight bonds, but generally, at t # 0, the bonds can be in either an intact or broken state. The transition of the system (the polyVGE molecule) is thus actually through the transitioins of its subsystems (bonds). The simulated reaction of polyVGE was achieved by allowing fixed time intervals, A t , to pass in series. After each, the state of each subsystem was deduced based on its prior state and its transition probability, which depended on both the magnitude of At and associated reaction rate constants. The state of the system was then recorded, and another At was allowed to pass. The subsystems were then tested again, and so on. The thus-implied Markov chain was of dimension sufficient to span the reaction time of interest, after which the initial polymer was reconstructed and a new Markov chain begun. Figure 2 thus illustrates three loops: the innermost loop spans the polymer length, as each bond was tested after each At; the second loop was reaction time, the passage
24 Ind. Eng. Chem. Res., Vol. 29, No. 1, 1990
Construction of Initial Molecule
P
i
!
Test for reaclion of all moieties
At
i
\
t
Yes
Store desired information. e,g, M.W. disi., conversion, product yields.
Average results of the N Markov Chains
,i
Figure 2. Logic of t h e M o n t e Carlo simulation.
of which mapped out a Markov chain reaction trajectory; the outermost loop was sample space, the dimension of which, NMc,usually was large enough to dampen the random fluctuations of a Markov chain. The Monte Carlo simulation results were thus the average of the predictions of a set of NMC independent Markov chains. The polyVGE reactant was thus an ensemble of irreducible subsystems (bonds and functional groups), each having an associated transition probability. PolyVGE itself had no unique transition probability, but rather, its new state was specified once the new states of all of its subsystems were determined. The transition of a subsystem from one state to another occurred when its transition probability was greater than a separately drawn random number. This essential step in the Markov chain simulation required a quantitative accounting of transition probabilities, to which we now turn our attention.
Transition Probabilities The subsystems of polyVGE could undergo, in general, irreversible or reversible, nth-order reactions that could be part of a parallel or series sequence. These prototype reactions were the modules from which the overall reaction
-
of the polymer was deduced. A B. The simulation of the first-order irreversible reaction of A to B averaged a collection of independent Markov chains. A common approximation for the transition probability for a first-order reaction is Pm = k m A t O(At), where O(At) 0 when At 0 (McQuarrie, 1967). In the present approach, where for computational motives it was useful to allow the fixed A t to be large, the exact solution for PAB was developed and exploited. The mathematics associated with Markovian processes (Appendix I) provided the exact functional form of the A to B reaction probability, Pm, as
-
+
-
pm = 1 - e - k . d t
(1)
where k A B is the rate constant for the reaction. The one-step transition probability (Pm)deduced above by matching stochastic and deterministic kinetics will be the basis for the construction of transition probabilities for other simple sequences where such a match is not easily accomplished. A B1and A B2. The probability of the disappearance of A by parallel fiist-order paths was determined by using the form of eq 1with its single rate constant (km)
-
-
Ind. Eng. Chem. Res., Vol. 29, No. 1, 1990 25
II
: ;
I S
,
NT. Numberol tranrdianr per A!
'.'.
0.01
'..
0
.
A
Za
.
-
1.5
kA
-
2
2.5
J
-
.
Thus, P = 1 - exp(-(kml replaced by the sum kml k,)At). During Monte+Car o simulation, a first-drawn random number (RN,) was compared to P, and for RN1 IP, a reaction of A occurred, to a yet unspecified product; for RN1 > P, A remained intact. When the former condition was met, a second random number (RN,) was drawn to determine the identity of the product. When RN2 I SB1= kABl/(kAB + kAB2), i.e., the selectivity or relative probability for the formation of B1, the product was B,; otherwise, the product was B,. For the more general case of M parallel reactions of A to Bl through BM, P = 1 - exp(CzlkmiAt) and SBi = k m i / ~ ~ l k A B iA. reacted when RN1 IP, and RN2 was then placed on a number line of unit length with M segments of length corresponding to the selectivity to each possible reaction product. The reaction product realized corresponded to the location along the line number on which RN2 fell. A B C. Although ostensibly a simple straightforward extension of the previous cases, the accurate formulation of a series reaction sequence by the above methods is deceptively complex. This is because the Monte Carlo simulation of A to B required only one transition in the fixed-sized time step. After its formation from the reaction of A, a B molecule was appropriately unreactive during the remaining part of the At in which it was formed. Following this logic for the A B C sequence would allow the subsequent reaction of B to C only in the next At. The erroneous overprediction of the yield of B could be especially large when kBc >> kAB. A Monte Carlo protocol was therefore developed to allow for more than one transition per At. This required the calculation of the transition probability for A C in one At, as outlined in Appendix 11. Equation 2 is the final result for k A B # kBc. Appendix I1 also summarizes the limiting cases for particular values of the rate constants
k"l
- -
--
-
(kAB
+
kBC
-
m,
kAB
=
kBC)*
Accounting for the possibility of two consecutive reactions in one time step was then accomplished by comparing a random number (RN) with both Pm and PAC. If RN < PAC,A reacted to C; if PAC 5 RN < Pm,A reacted to B; and if RN 2 Pm, A did not react. Figure 3 summarizes the Monte Carlo simulation of A B C for kBC/kAB= 10. This example was chosen to illustrate the influence of step size (At) for both the approximate method, where only one transition per At was allowed, and the exact method, where two transitions were
--
0.5
1
1.5
3
Figure 3. Monte Carlo simulation of the consecutive sequence A B C: effect of the simulation parameters.
+
0
1
Deterministic Solution
Deterministic Solution
~
1
A
0.94
o.88 0.86 0
0.5
-
Symbol
0.1 0.3 0.3 0.5 0.1 A
-
z
0
transitions p e r l
0.98
i
8
'
N T = Numberol
I
2
2.5 kAB
3
3.5
4
4.5
5
t
Figure 4. Monte carlo simulation of reversible first-order reactions: effect of transition order NT.
possible in a given At. When only one transition was allowed, the x2 test showed that the simulation provided an acceptable representation of the kinetics only for small values (
