Energy & Fuels 1994,8, 570-575
570
CPU Issues in the Representation of the Molecular Structure of Petroleum Resid through Characterization, Reaction, and Monte Carlo Modeling Thomas F. Petti,? Daniel M. Trauth, Scott M. Stark, Matthew Neurock,t Muzaffer Yasar, and Michael T. Klein* Center for Catalytic Science and Technology, Department of Chemical Engineering, University of Delaware, Newark, Delaware 19716 Received October 4, 1993. Revised Manuscript Received January 11, 1994’
The optimal usage of CPU resource for the optimization of resid attribute probability density functions was examined. Alternative strategies involving the use of reaction simulation and data, NMR and increased sample size were examined. The present results suggest optimization with NMR data, large molecular sample sizes, and the absence of reaction simulation as the most efficient strategy.
Introduction The vigor in the development of molecule-based petroleum feedstock reaction models is a result of two broad forces. First, being based at the molecular level and therefore embracing fundamental physical-organic chemistry, molecular reaction models are reasonably anticipated to provide better predictions and allow for extrapolation. Second, new questions of unprecedented molecular detail are being asked of reaction models. The Clean Air Act amendments of 1990 ensure that molecular groups’ and individual species’ concentrations will be relevant issues in refinery operations. At a conceptual level, molecule-based reaction models are a set of conservation equations representing the reactor and the rate law. A plug flow or batch model will generally involve the solution of a set of coupled ODES,with rate constants (K) as parameters and initial conditions on the variable species’ concentrations. Generally, analytical chemistry provides the initial conditions, and model compound data or direct optimization provides values for K. Arguably, the situation becomes increasingly challenging as the boiling point of the feedstock increases. In resid reaction modeling, for example, the molecules are not volatile and therefore are not subject to the analytical measurements that provide molecular information for even gas oils. In a modeling sense, the initial conditions are very poorly known. Indeed, it may be further argued that the model parameters (k)are often better known than the initial conditions. For example, the reaction kinetics of alkylbenzenes are likely known as well as the number of alkylbenzenes present in a resid. This paper is aimed at exploring the thus derived notion that resid reaction behavior can assist in the determination of resid structure. We have developed statistical methods to construct a representation of resid molecules from measurements of resid attributes.1-3 Recognizingthat molecule-by-molecule Present address: W. R. Grace & Co., Columbia, MD 21044. t Present address: Eindhoven University of Technology, Eindhoven, The Netherlands. e Abstract published in Advance ACS Abstracts, March 15, 1994. (1) Trauth, D. M.; Stark, S. M.; Petti, T. F.; Neurock, M.; Yasar, M.; Klein, M. T. Prepr. Diu.Pet. Fuel Chem. 1993. t
0SS7-0624/94/250S-0570$04.50/0
measurements in resid feedstocks would be difficult, we sought to define and characterize the attributes of the resid molecules that would allow construction of a statistical representation of a given feedstock. This approach incorporated methods for both the transformation of the basic analytical chemistry into quantitative probability density functions (pdf s) and the subsequent construction of the molecular representation. Figure 1 illustrates how these pdfs were used to construct a resid molecular representation. A series of random numbers were chosen to “sample” the pdf s for each of the structural attributes in order to build individual molecules. For example, a first random number is chosen to decide the molecule type (paraffin, aromatic, or naphthenic), which in this example is a naphthene. The second and third random variables define the number of rings and their configuration which makes up the core of this naphthene. The final two distributions are then sampled to determine the number of aliphatic substituents and carbon chain length of each. This process was repeated many times until a sample population of molecules was created to represent the entire resid. A more detailed account of this process is described el~ewhere.~J Clearly, the ultimate molecular representation is defined by the functional type and the parameter values of each pdf. In this work, we examined the value of the information from a well-defined, standard pyrolysis test in the iterative determination of pdf parameter values. The pdf s were approximated by the chi-square distribution, since its shape retained the qualitative features expected to mimic resid. Thus, two parameter values (mean and minimum) for each attribute provide a quantitative representation of the resid. These parameters were determined via optimization to minimize the differences between the properties of the molecular representation and the measured properties of the resid. In this context, resid reactivity is simply another measurable (2) Neurock, M.; Nigam, A.; Libanati, C.; Klein, M. T. Chem. Eng. Sci. 1990,45, 2083.
(3) Neurock, M.; Nigam, A.; Trauth, D.; Klein, M. T. Molecular Representation of Complex Hydrocarbon Feedstocks through Efficient Characterization and Stochastic Algorithms. Energy -- Fuels, to be submitted for publication.
1994 American Chemical Society
q)
CPU Issues in the Representation of Structures
:K-iK
00
3
Energy & Fuels, Vol. 8, No. 3, 1994 571 2
4
6
8 1 0
b
*of Naphther& Rlngr
/
1 Vumber of Sidechimi
conflgmuan Of Nnphlhenlc Rtngs
5
8
lengih of Sldshalnr
m I
Figure 1. Stochastic construction of resid molecule; pdfs are sampled to determine the value for each structural attribute. Table 1. Structure-Reactivity Polanyi Parameters for Hydrogen Abstraction and &Scission of Paraffins, Alkylaromatics, and Alkylnaphthenes' H-abstraction @-scission compound class E. a E. B 0.42 12.0 0.58 paraffins 14.8 0.39 13.2 0.64 alkylaromatics 17.7 alkylnaphthenes 16.2 0.13 11.7 0.28
resid property, like proton environment (NMR), MW (VPO), SARA fractions (chromatography), and boiling point distribution (simulated distillation). Thus, a reaction model was used to obtain predicted resid reactivity. This was used as additional information in the optimization objective function that provided the optimal pdf parameters. This paper examines the usefulness and CPU implications of taking this approach.
The Reaction Subroutine The reaction of resid molecules was chronicled in terms of the reactions of their reactive sites. These were, in turn, inferred from the reactions of relevant model comp0unds.e A four-component kinetically coupled mechanistic simulation of pyrolysis (4CM), described in detail by Nigam and Klein,7y8 was used as the ultimate basis for modeling the reactions of the set of simulated resid molecules. The four base components in the mechanistic model were pentadecylbenzene (PDB),tridecylcyclohexane (TDC), decyltetralin (DT), and octane (0). The 4CM simulation allowed for dealkylation, dehydrogenation of naphthenic rings, and paraffin and olefin cracking. Concentrations of each of the base 4CM components were used in the radical balances for approximately 30 radical species. The 30 radical components were chosen as a representative basis for the radical fragments which exist under pyrolysis conditions. A set of six Evans-Polanyi relationships were used to compute the kinetics for the governing hydrogen abstraction and /3-scission steps for paraffins, alkylaromatics, and alkylnaphthenics. These optimal structurereactivity parameters are summarized in Table 1and serve as the kinetic input. The mechanistic simulation solvesthe set of radical balance equations to determine the rate of reaction for each of the 4CM components. These results were then used to estimate first order rate constants for each component. (4) Savage, P. E.; Klein, M. T. Znd. Eng. Chem. Res. 1987, 26, 488. (5) Savage, P. E.; Klein, M. T. Ind. Eng. Chem. Res. 1988, 27, 1348. (6) Smith, C. M.; Savage, P. E. AZChE J. 1991,37,1613. (7) Nigam, A,; Klein, M. T. Ind. Eng. Chem. Res. 1993,32, 1297. (8)Nigam, A. Ph.D. Dissertation, University of Delaware, 1992.
Figure 2. Reactive sites on 2-nonyl-4-propyltetralin.
The input to the 4CM simulation required determining the concentration, ci, of each of the four 4CM model molecules. The following set of rules was used to link the set of resid molecules to the 4CM components. In the resid, (1) a substituted aromatic ring was considered a PDB moiety; (2) a substituted naphthenic ring was considered a DT moiety if there were any aromatic rings in the molecule; otherwise, it was considered a TDC moiety; and (3) a paraffin or olefin molecule was considered an 0 moiety. Applying rules 1-3 to the set of resid molecules provided the number of moieties, nij, of each 4CM component i on each resid molecule j . The molecular weight (mwj) of each resid molecule and the experimental resid density ( p ) allowed estimation of the concentration as the total number of 4CM molecules of type i divided by the volume occupied by all of the resid molecules:
-
no. molecules
nij
P
ci =
J=1
i = {PDB, DT, TDC, 0) (1)
no. molecules J-1
mwj
The rate constants determined from the 4CM were subsequently applied to the reactive moieties of the resid molecules on an individual site basis, taking into account modifications for the number of aromatic rings, side-chain length, and selectivities. This is illustrated in Figure 2, where 2-nonyl-4-propyltetralin is labeled with its 13 possible reaction sites. Dehydrogenation may occur at reaction site 1. The results of model compound experiments in the literature9 were used to provided first estimates of the kinetics for dehydrogenation for cyclohexane, tetralin, di- and tetrahydrophenanthrene, and diand tetrahydroanthracene moieties. Reactions of the alkyl side chains presented several additional challenges. First, the length of a side chain was likely to be different from the base 4CM molecule. Relationships developed by Savage and Klein415were used to account for the carbon number dependence of reactions of the alkyl groups. Equations 2 and 3 were used for PDB and DT (or TDC) moieties, respectively, multiplier =
chain length 15
[
1
112
Regarding reaction selectivity, side chains are capable of cracking at several primary locations: at the ring, a to the ring, /3 to the ring, and the minor equiprobable route at all other locations. Factors affecting the cracking location selectivity include moiety type and the number of aromatic rings.@ Returning to the example of Figure 2, the overall rate for the nonyl side-chain cracking was corrected by (9) Poutsma, M. ORNL/TM-10637, Chemistry Division, Oak Ridge National Laboratory, 1987.
Petti et al.
572 Energy & Fuels, Vol. 8,No. 3, 1994
multiplying the DT overall rate constant by the carbon number relation. Individual site rate constants were determined by multiplying the corrected overall rate by the location selectivity. For example, the rate constant for site 2 was found by multiplying the overall rate by the percentage of products formed by cracking at the ring. A similar procedure was followed for the rest of the nonyl side-chain sites and for the propyl side chain. Each atom of every molecule is sampled to determine whether it can be a reactive site. If so, the appropriate rate constants for alkyl cracking, or naphthene ring dehydrogenation, are then computed. This specification of "site" kinetics for every molecule allowed construction of a cumulative probability distribution comprising all possible reaction events (Cki). The reaction simulation subroutine was therefore based on a variable time step Monte Carlo approach,2J0-12where eq 4 provided the time required for an even to occur. Random numbers were
4 CM components
for each of tbe 4 CM comwnents
J I
neventa
At = -ln(l - RN)/
ki
(4)
i=l
chosen to determine which reaction event occurred and the elapsed reaction time. If the net elapsed time was greater than the specified total reaction time, the reaction subroutine was finished and the next set of moleculesbegan reacting. If not, the reacted molecule was updated and the reaction loop was repeated. Figure 3 outlines the main components of the reaction model. The simulation proceeded by reacting the set of resid molecules over the time period of interest. This provided a complete description of the molecular composition versus time. The reaction product molecules were then lumped according to type and boiling point correlations to provide SARA and SimDisfractions. SARAwas computed directly based on each molecule's structure. The final reaction results were represented by the profile of each simulated distillation fraction over time. The boiling point correlations used to determine the SimDis fractions were derived by adding a ring core contribution with an alkyl chain contribution. The ring core component is given here in eq 5 , where NAR refers to the number of aromatic rings and TR to the total number of rings.
+ 486.5[log(TR)l) + [%](71.95 + 4ll.l[log(TR)I) ( 5 )
bpcore = (%)(80.93
The component for aliphatic chains is described via a set of different correlations depending on the nature of the core. These correlations are functions of the side-chain length and in some instances the number of aromatic and naphthenic rings. For example, a molecule with a core of more than two aromatic rings and aliphatic substituents less than five carbon atoms long would use the following correlation for the alkyl substituents, where SCL refers to the side-chain length. bpchain = 9.0 + 7.3(SCL) + (SCLI2 The complete set of correlations used, along with summary
I no Yes
Figure3. The stochastic,kinetically coupled reaction subroutine.
of how they were derived and how well they predict hydrocarbon boiling points, is given e1~ewhere.l~
Using Reaction Simulation for Structure Characterization The reaction simulation was used in the optimization of structural pdf parameters. In order to include reaction information in the optimization, terms were added to the objective function that compared simulated reaction product properties to those measured for the real resid. The properties considered were the simulated distillation weight fractions produced by the pyrolysis of the resid at several reaction times. Batch tubing bomb reactions were performed to pyrolyze an Arabian Light resid feedstock.14 Its initial SARA and SimDis characterizations are provided in part A of Table 3. Products were separated by soxhlet extraction and quantified in terms of gas, maltene, asphaltene, and coke ~
(10) Gillespie, D. T. J. Comput. Phys. 1976,22,403. (11) Gillespie, D. T. J. Phys. Chem. 1977,81, 2340. (12) Neurock, M. Ph.D. Dissertation, University of Delaware, 1992.
(13) Trauth,D. M. Ph.D. Dieeertation, University of Delaware, 1993. (14) Trauth,D. M.; Yaear, M.; Neurock, M.; Nigam, A.; Klein, M. T. Fuel Sci. Technol. Int. 1992, 10, 1161.
Energy &Fuels, Vol. 8, No. 3, 1994 573
CPU Issues in the Representation of Structures 100
Table 2. Attribute Pdf Parameters Optimized with and without Reaction Data
attribute no. of aromatic rings no. of naphthenic rings no. naphthenic rings on aromatic core no. of side chains length of side chains
minimum 1 1
0 0 1
mean with without reaction reaction 4.7 5.5 7.0 7.8 3.2 3.6 4.5 8.2
4.8 7.8
80
#
5 60
B V
'ii 2 a
40
20
solubility lumps. The maltene fraction was then characterized by simulated distillation. The resulting optimization objective function included simulated distillation terms associated with the weight fractions at several reaction times:
0.05 MW,,,
aexp
- Hapred
0.02 H to Celp - SaraWt,,,
0.03
0 0
20
40
60
100
80
Experimental Wt % Figure 4. Comparison of simulation results obtained using the pdf parameters optimized with and without reaction product properties in the objective function. Filled symbols refer to simulation results with reaction. Unfilled symbols refer to simulation results without reaction. Circles, squares, triangles, and diamonds refer to initial SimDis weight fractions, SimDis weight fractions after 30 min of reaction, SimDis weight fractions after 90 min of reaction, and initial SARA fractions. Table 3. A Comparison of Experimental with Simulation Predicted SimDis, SARA,Molecular Weight, and H/C Ratio for an Arabian Light Feedstock.
The results of this approach are considered in the next section, in addition to the cost in terms of CPU time.
analytical measurement
expt
Effect of Reaction Data
gas 90-190 O F 190-380 O F 380-520 O F 520-610 O F 610-800 O F 800-1000 O F 1000 O F + saturates aromatics/resins asphaltenes MW HC
O.oo00 O.oo00 O.oo00 O.oo00 O.oo00 1.9000 9.4000 88.700 27.000 71.000 2.oo00 622 1.60
gas 90-190 O F 190-380 OF 380-520 O F 520-610 O F 610-800 OF 800-1000 O F 1000 O F +
3.8000 0.4oo00 4.1000 5.9000 5.3000 15.500 15.100 50.000
The optimized pdf parameters are summarized in Table 2 for five structural attributes (the number of aromatic
rings per aromatic molecule, the number of naphthenic rings per saturated molecule, the number of substituted naphthenic rings per aromatic molecule, and the number and length of alkyl side chains) for the Arabian Light resid. These pdf parameters were determined by the simulation and optimization algorithm and represent the optimum values determined first using pyrolysis reaction data and then excluding the reaction data. Comparing the two approaches, most of the pdf s optimized using the reaction data were quite similar to those computed without reaction information. However, differences exist. For example, the pdf describing the number of naphthenic rings shifted to a lower mean when reaction data were excluded. The question for optimization of CPU resources was whether the pdf s found using the pyrolysis data resulted in a more accurate representation of the actual resid. Considerable CPU investment was required to simulate the resid reaction, and further analytical effort was necessary to pyrolyze the feedstock resid and analyze the resulting products. Comparison of the properties of the simulated and real feeds probed the accuracy of the pdf s. SARA and SimDis results are summarized in Figure 4, for the sets of pdf s optimized with and without pyrolysis reaction data (Table 2). The results shown in Figure 4 and in Table 3 indicate that the two sets of pdfs (obtained with and without reaction data) represent the feedstock to nearly the same accuracy. The data all lie within an acceptable band about the diagonal; the points for each case appear to be within approximately the same degree of error. This was further confirmed by the calculated standard deviations: 6.1 for
simulation without reaction with reaction A O.oo00 O.oo00 O.oo00 O.oo00 O.oo00 O.oo00 O.oo00 O.oo00 O.oo00 O.oo00 2.7000 3.5000 21.3000 23.600 75.900 72.900 27.500 26.500 70.700 72.700 1.9000 0.8000 691 656 1.57 1.58
B 5.6000 1.3000 6.3000 5.4000 1.9000 7.1000 22.100 50.300
3.7000 1.1000 5.7000 5.3000 1.6000 6.9000 24.700 50.900
C '8.3000 16.800 13.200 0.6oo00 3.oo00 2.5000 11.100 8.9000 11.600 11.500 8.2000 8.6000 7.6000 3.oooo 3.1000 15.000 10.700 12.000 8.7000 19.400 21.700 39.500 27.300 27.800 aThe three c( unns refers to exnerimentc data, simulation predicted without the reaction subroukne, and simulation with the reaction subroutine. The results here refer to (A) the initial feedstock, (B)products after 30 min of reaction, and (C) products after 90 min of reaction. gas 90-190 O F 190-380 O F 380-520 O F 520-610 O F 610-800 O F 800-1000 O F 1000 O F +
the pdf s with the reaction data and 5.6 without the reaction data. The difference between these statistics is not significant. These results led to the conclusion that the additional analytical effort and CPU time required to
Petti et al.
574 Energy & Fuels, Vol. 8, No. 3, 1994 $
i
lo4
distribution, with an average of the true global mean and a sample variance equal to the true variance divided by the sample size. This can be represented as a normal distribution ( N (mean, variance)) describing the predicted property: Propprd
I 0.1
'11"'
10
'
'
""""
100
'c"'"i
'
lo00
'
"""'i
104
0.1 '"J
los
No. of molecules
Figure 5. Analysis of the effect on the optimization objective function (chi-square),ita standard deviation,and the associated CPU time of the number of molecules used in the simulation. perform the optimization using the pyrolysis data was not worth the investment, since significant improvemenb in the simulation accuracy could not be guaranteed.
Sample Size Effects As previously discussed, the feedstock simulation consisted of stochastically sampling the various attribute pdf s many times to construct a representative population of molecules. The CPU time required to build the molecules in the feed was proportional to the number of molecules simulated. For this reason, it was efficient to simulate the fewest molecules possible that accurately described the resid. The importance of the sample size is evident when it is recognized that the properties calculated for the simulated feed were global averagestaken over each molecule present. For example, the average molecular weight of the feed was determined by adding the individual molecular weights of each molecule in the simulation and dividing by the total number of molecules. Increasing the number of molecules that were built, increased the number of times that the pdf s were sampled, which resulted in a more accurate estimate of the average molecular weight. Clearly, the generation of only a single molecule permitted the possibility that the molecular weight would be far from the "true" distribution mean. At the other extreme, simulation of an infinite number of molecules would guarantee convergence to the true mean by the central limit theorem. However, the requirement of a reasonable CPU investment conflicted with our accuracy goals, demanding a compromise for the sample size. Figure 5 illustrates the effect of sample size on the CPU time and the objective functionaveragevalue and standard deviation. This was probed by performing 50 replicate simulations using the same set of pdf parameters for each sample size. The standard deviation (variability) of the objectivefunction decreased with an increase in the sample size. The average value of the objective function also changed as the number of molecules increased. This was a significant result since it suggested that a minimum sample size was required to perform an accurate simulation for reasons other than stochastic variability. Scrutiny of the nature of the objective function (eq 6) explained these sample size effects. The objective function components can be thought of as normal random variables. Even if the true simulated properties were not normally distributed, the central limit theorem states that each sampled average property will tend toward a normal
-
N(Prop,,
a2/n)
(7)
where Propprd is the random variable describing the predicted property (e.g., MW), Proptrueand u2 are the mean and standard deviation of the actual simulated property (corresponding to infinite sample), and n is the sample size used to determine Propprd. Therefore, the larger the sample size, the smaller the sample variance; in other words the averaged property is a better estimate of the true simulated property, Prop,,,. Given that the predicted properties were nearly normal, then the distribution of the objective function would follow the chi-square distribution, x2(mean),(hence the objective function's name). Since the variance of each of the predicted properties was inversely proportional to the sample size, we expected the average of the objective function (Objective) to be inversely proportional to the sample size: Objective
a
l/n
(8)
This is the behavior reported in Figure 5. The initial slope of the objective function average on the log-log plot was found to be approximately -1; the average decreased to an asymptotic valuerougly with the reciprocal sample size. For the chi-square distribution, the variance is proportional to the mean. The expected behavior of the standard deviation was therefore inversely proportional to the square root of the sample size: s (Objective) a l / h
(9)
This is also seen in Figure 5. The standard deviation decreased linearly on the log-log plot with a slope of approximately -0.5 (inverse square root dependence). This analysis suggested that an optimal sample size would be between 10oO and 10 OOO molecules. This would ensure that a reasonable property average was obtained, corresponding to the leveling off of the objective function average. Additionally, the larger the sample size, the lower the expected variability. The standard deviation analysis suggested that using a sample size of 10 000 ensured that the standard deviation of the objective function would be less than one. This was a significant improvement over the standard deviation for the 1000 molecule simulations, which was around 5.
Conclusions The use of experimental and simulated pyrolysis data in the objective function may provide insight into the structural character of the resid molecules. However, in order to use this information, additional simulation time was required to react the molecules stochastically. Additional experimental effort was also necessary to characterize the reaction products of the real feed. This effort may not be warranted because the accuracy of the resulting pdf s was not significantly improved. The CPU time required to react the feedstock was better invested in using a larger sample size of simulated
Energy & Fuels, Vol. 8, No. 3, 1994 575
CPU Issues in the Representation of Structures molecules. The sample size had a dramatic effect on the outcome of the optimization objective function both in average and variability (standard deviation). A sample size of 10 000 molecules offered an accurate estimate of the feedstock properties at a reasonable computation expense. Although the reaction model did not prove useful for pdf parameter optimization, the fact that the predicted reaction behavior did not adversely affect the optimization offers a validation of the reaction model. This model can therefore be used in studying the reaction behavior of resids that have been previously characterized in terms of structural attribute pdf s.
Nomenclature and Symbols 4CM Ci
DT hi, mwi MW(VP0)
four-component kinetically coupled mechanistic simulation of pyrolysis concentration of each of the four component moieties in the feed or reacted fraction decyltetralin kinetic rate constants molecular weight of each molecule molecular weight as measured by vapor phase osmometry
ni
N, n NAR
0 Objective
ODE Pdf PDB
&Objective) SARA SimDis TDC TR MW(VP0)
At U2 u2/ n
P X2
number of each of the four component moieties in the feed or reacted fraction number of molecules, sample size number of aromatic rings per core octane the objective function ordinary differential equations probability density function pentadecylbenzene random variable describing predicted property mean of the actual predicted property random number the standard deviation of the value of the objective function saturate, aromatic, resin and asphaltene fractions of a petroleum feedstock simulated distillation fractions tridecylcyclohexane total number of rings molecular weight as measured by vapor phase osmometry variable reaction time step standard deviation variance experimental resid density objective function (chi-squared probability distribution)
