Advanced Chemistry Surrogate Model Development within C3M for

Article pubs.acs.org/IECR

Advanced Chemistry Surrogate Model Development within C3M for CFD Modeling, Part 1: Methodology Development for Coal Pyrolysis Dirk Van Essendelft,*,† Tingwen Li,†,‡ Philip Nicoletti,†,‡ and Terry Jordan† †

National Energy Technology Laboratory, Morgantown, West Virginia 26505, United States URS Corporation, Morgantown, West Virginia 26505, United States

‡

ABSTRACT: A generalized method was developed to implement complex chemistry from third-party computational coal chemistry tools within multiphase, reacting computational fluid dynamic (CFD) codes. This method involves generating response surfaces that represent the computational coal chemistry tool output and using them to build a comprehensive surrogate model which can be easily incorporated into a CFD code. In this first part of a series of work, the method was applied to coal pyrolysis of Powder River Basin coal using PC Coal Lab (PCCL) and Carbonaceous Chemistry for Computational Modeling (C3M) to generate the series of response surfaces which form the basis of the surrogate model that can be implemented in CFD. This method has the benefit that local environmental variables (such as heating rate, final pyrolysis temperature, local pressure, and particle diameter) can be incorporated as part of the chemistry model in CFD on a cell by cell, time step by time step basis which creates the potential to significantly increase the accuracy of and reduce uncertainty in modeling coal chemistry within CFD. As demonstrated with this surrogate modeling method, the modeler is alleviated of the responsibility to predetermine a heating rate and reactor temperature for pyrolysis. This particular part of the work focuses on the theory and development of the surrogate model method. The methods used to generate the response surfaces, build them into a coherent surrogate model, and estimate the thermochemical properties of pseudocomponents are discussed. Future work and publication will focus on verification within CFD, validation against experimental data, multidimensional surrogate model development, and application of the surrogate model methodology to other reaction chemistry.

■

INTRODUCTION The thermal pyrolysis of coal is among the most complex reaction systems in modern science.1−8 The chemical/physical structure of coal is very complex and there is a large degree of variation within coals. Not only does coal significantly vary between ranks, it can have wide variations within the same coal seam. These variations and complexities make the task of accurately modeling coal pyrolysis difficult. Furthermore, coal pyrolysis has been shown to be quite sensitive to environmental factors such as heating rate, pressure, and ultimate temperature.9−11 That said, coal pyrolysis remains an extremely important part of most coal-based energy conversion processes such as combustion and gasification. In recent times, several computer codes have been developed that attempt to predict bulk kinetics and gas composition during pyrolysis (e.g., CPD, FG-DVC, and PCCL).1,7,8,12−35 These tools rely on large databases of coal experiments or structural information supplied by the user to generate kinetic expressions and gas yields. Most of the experimental information available was conducted on either wire grid experiments with a prescribed temperature profile or with a drop tube experiment. Within reason, these computational tools have proven to be useful as predictive tools, though some level of disagreement exists between them.5,36−40 Numerous efforts have been made to integrate kinetics from these detailed pyrolysis models into CFD.36−38,41−48 The key challenge to implementing complex pyrolysis sub models such as CPD, FG-DVC, and PCCL is the dependence of rate and yield on particle temperature history and the lack of the ability to know or reasonably estimate the temperature history before a CFD simulation is started. This has lead researchers to take an © 2014 American Chemical Society

iterative approach (either manual or automated) in which kinetics are obtained from one of the sub models on a best guess basis to start followed by CFD simulations and subsequent refinements in the temperature history input to the sub model followed by further CFD runs and further refinements. This process is continued until the temperature history stabilizes to some defined tolerance. In recent times, at least two groups have applied algorithms to automate the iterative solutions. Hashimoto et al. have developed a “tabulated−devolatilization−process model (TPD model)”.37,38 The TPD model involves generating a table in which temperature history and reaction parameters are correlated before running the CFD. The kinetic parameters are either taken directly from the sub model or are calculated from fits to their output. The iteration is conducted between the CFD and the tabulated results. Vascellari et al. took a similar approach by substituting a preprocessed table with direct application calls to the sub model software and minimizing scheme using a genetic algorithm to find the global minima during the iterations.36 Although these approaches arguably produce accurate results, the requirements for multiple CFD runs in an iterative loop make them too costly to be used in industrial scale work where a single CFD simulation can take multiple months on a super computer. It is desirable to utilize and develop a predictive chemistry model that does not require iterations of CFD runs. Recently, Maffei et Received: Revised: Accepted: Published: 7780

October 8, 2013 February 19, 2014 March 26, 2014 March 26, 2014 dx.doi.org/10.1021/ie402678f | Ind. Eng. Chem. Res. 2014, 53, 7780−7796

Industrial & Engineering Chemistry Research

Article

Figure 1. Schematic diagram of surrogate model generation algorithm in C3M.

Chemical Reaction Model. In order to allow the chemical stoichiometry to vary in response to application input variables, it was necessary to define the chemistry on a per gas basis rather than on a volatile consumption basis. To do this, coal was modeled as a set of solid species consisting of ash, moisture, carbon, hydrogen, oxygen, nitrogen, and sulfur. The amount of each species can be determined from the coal’s ultimate and proximate analysis. Defining coal in this way allows for easy element tracking and the definition of simple elementary reactions within the CFD model. For example, the formation of methane during pyrolysis is defined as follows:

al. proposed a simplified single step model in which the reaction and stoichiometric parameters are allowed to vary with coal composition.48 Although this model can be directly applied to CFD without the use of a detailed pyrolysis code (and thus no need for iteration), it does not account for environmental factors such as heating rate, temperature, and pressure. A more ideal solution would be to update the CFD model with chemistry on a cell-by-cell basis at each time step during the simulation with information obtained directly from a detailed pyrolysis models. This would eliminate the need to preemptively determine the heating rate and reactor temperature or conduct additional iterations. One method would be to simply make an application call at each cell every time step. However, this is excessively time-consuming and may not work at all due to stability concerns in many of the pyrolysis model codes. An alternative is to generate a closed form surrogate model that can represent the computational pyrolysis model within the CFD model. This can be done by creating a design space that covers the range of values expected within an industrial model and fitting expressions to the resulting parameters. It is now possible to do this within the National Energy Technology Laboratory’s Carbonaceous Chemistry for Computational Modeling software (C3M). The software is available at mfix.netl.doe.gov on the World Wide Web.

C(coal) + 4H(coal) → CH4(gas)

(1)

Every gas evolved during pyrolysis can be defined in a similar manner and a generalized chemical equation can be written αC(coal) + β H(coal) + γ O(coal) + δ N(coal) + εS(coal) → Cα HβOγ NδSε

(2)

The only pyrolysis reaction that poses a challenge in this methodology is tar formation. Tar is not a pure compound but is made up of a wide variety of compounds. Therefore, it is not possible to come up with an elementary chemical formula. It is possible, however, to define an elemental composition for tar (ultimate analysis), Ui, where i is the elemental species C, H, O, N, or S, and calculate the reaction stoichiometry from that. The sum of all Ui’s is unity. Therefore, the stoichiometry in eq 2 can be defined for tar

■

SURROGATE MODEL GENERATION METHODOLOGY The developers of C3M have taken the “black box” approach to surrogate model generation as seen in Figure 1. First, the application inputs must be specified such as reactor temperatures, heating rates, pressures, particle sizes, and so forth. These can be in the form of a uniform or normal distribution by specifying the minimum and maximum range or the mean and standard deviation, respectively. In addition, the sampling method must be specified. This information is then fed through an uncertainty quantification (UQ) engine such as PSUADE to generate the sample design. Each point in the design space represents a specific set of input parameters and variables for the simulation model (pyrolysis model in this case). C3M can then run the sample design by making application calls to the pyrolysis model for each sample design point and collecting the results in a formatted output file. Once complete, a closed form surface expressions for each response variable can be generated. Once the collection of surfaces has been generated, they can be combined to form a comprehensive surrogate model.

α=

β=

γ=

δ=

ε=

M Wtar*UC A WC

(3)

M Wtar*UH A WH

(4)

M Wtar*UO A WO

(5)

M Wtar*UN A WN

(6)

M Wtar*US A WS

(7)

Here, Mw and Aw are molecular and atomic weight, respectively. 7781

dx.doi.org/10.1021/ie402678f | Ind. Eng. Chem. Res. 2014, 53, 7780−7796


Article

The rate for each gas species reaction can be defined from the pyrolysis model ri = ρε x A e−EaTot / RTs(XFinal − XChar)Yi s s CHONS Tot

(8)

Here, ri is the rate of gas production rate for a given species, i, ρs is the solids density, εs is the solids volume fraction, and xCHONS is the mass fraction of the C, H, O, N, and S species in the coal/ Char on a local bases. ATot and EaTot are the pre-exponential and apparent activation energy for the total rate of pyrolysis predicted by the pyrolysis computational model. R and Ts are the gas constant and solids temperature, respectively. X is the yield of char and the subscripts “Final” and “Char” indicate the final char yield predicted by the surrogate model and the instantaneous char yield at any point within the reacting system. Yi is the yield of a given species predicted by the surrogate model. The instantaneous char yield can be calculated from the ratio of the mass fraction of char to ash at any point in the system relative to the ratio at the inlet of the coal injection. This is expressed in eq 9. It should be noted that eq 9 is only valid if the ash in the coal is unreactive in the system. If this is not the case, or cannot be approximated as such, one would have to derive a different measure of char yield

( ) =1− ( ) xChar xash

XChar

xChar xash

Figure 2. General information flow diagram for application of the pyrolysis surrogate model in CFD.

detailed below. With the heating rate at each cell and a measure of reactor temperature calculated, the surrogate model variables (Atot, EaTot, XFinal, and Yi) can be calculated for each cell or particle. The only steps remaining are to evaluate the reaction rates for each product gas (eq 8) and to calculate the local char yield via eq 9. With all of the variables in eq 8 now available on a cell/ particle basis, the reaction rates can be calculated for each cell/ particle. Once the rates are calculated, the CFD code can move on to the chemistry solver and the next time step can be evaluated. Evaluating Heating Rates in CFD. It is straightforward to calculate the heating rate in the Lagrangian approach where temperature history of each coal particle/parcel is tracked. As no temperature history of individual coal particle is available in the Eulerian approach, only a mean heating rate of coal particles in each control volume can be estimated. In the Eulerian specification of the flow field, the heat rate is calculated as

Dry

Dry,Initial

(9)

In this model, the variables Atot, EaTot, XFinal, and Yi are all closed form surface expressions that were fit to the output of the pyrolysis model within C3M. These variables could be functions of heating rate, reactor temperature, reactor pressure, and particle diameter depending on what is chosen as application inputs to the UQ engine within C3M. With the use of the surface expressions, every variable in eq 8 can be calculated in each cell at each time step within a CFD model. Assembling the Information in a Coherent Surrogate Model. Within many reacting systems, the reactor pressure and inlet particle size are known with a relatively high degree of confidence and often do not need to be part of the surrogate model. The final reactor temperature and heating rate (temperature profile) are harder to predict accurately and cannot be reasonably predicted ahead of time. Therefore, many simulations will only need to have the heating rate and reactor temperature as application inputs for the pyrolysis reactions. With this in mind, the generalized procedure for applying the surrogate model in a CFD system can be shown schematically in Figure 2. In this example, the surface expressions that make up part of the surrogate model are functions of the local heating rate and the final reactor temperature. At the start of each time step, a defined domain section within the simulation is used to average the solids temperature and calculate a representative reactor temperature. It is up to the user to define this section, but the reactor temperature calculated should represent the final temperature of the solids during pyrolysis. Also at the start of each time step, the solids heating rate should be evaluated in some way within each cell or particle. In most Lagrangian CFD codes, the heating rate can be directly calculated from the particle temperature at the current and previous time step as well as the known time step. It is also possible to estimate the heating rate for Eulerian codes from the energy equation as

HR =

dTs ∂T = s + us·∇Ts dt ∂t

(10)

where Ts is the solid temperature and us is the solid velocity vector. From eq 10, it is clear that the mean heating rate of particles in each control volume is related to the temporal variation of temperature and the spatial convection. Necessity To Time Filter Heating Rate and Reaction Temperature. Because the heating rate is a strong function of the change in solid temperature from the current and past time steps and directly affects the speed of the reaction and thus the energy generated in the reaction, there is a strong possibility of forming high frequency, discrete feedback noise in the transport variables when applying this method. This is especially true when the reactions are endothermic. To understand this, imagine a single cell in which coal is being heated from ambient temperatures to gasification temperatures with a condition where the sum of all the reaction energies is endothermic. At time zero, the previous solid temperature is initialized to the current solid temperature. The heating rate is calculated to be zero, which relates to a very small reaction rate, and thus, a very small cooling source is input into the energy equation. At the same time, assume that some convection occurred that was greater than the reaction energy. At the next time step, solving the energy equation results in the solid temperature rising. Because of the structure of the method, the reactions cannot 7782



Article

remaining solution is to cut the time step. It may also prove to be advantageous to create a θ function that is responsive to the local heat of reaction such that the amount of filtering is minimized in each cell. Finally, it may also be important to apply filtering to the calculation of reactor temperature depending on how the reactor temperature is estimated and how stable the calculated value is with respect to time. It is conceivable to model a reaction system in which the solid temperature in the evaluation domain varies significantly with time. This scenario is often encountered in cases where the hydrodynamics are unsteady. In this case, it is undesirable to force the entire reaction system to respond to the instantaneous domain temperature but rather to respond to a low frequency filtered result (essentially a time averaged result). In cases like this, eq 11 can be applied with a large θ and an input signal of the instantaneous domain temperature.

respond until the following time step. Therefore, the heating rate due only to convection was likely too high. With this higher heating rate, the reaction rate calculated on the next time step is also artificially high, which could create a large source of cooling in the energy equation. The convection heat source could now be much smaller than the reaction energy. If so, solving the energy equation will result in a temperature drop. This results in a negative heating rate, which results in a negligible reaction rate and, thus, a negligible cooling source in the energy equation. The process can now repeat indefinitely, creating noise. This is obviously not physical but is an aberration of the numerical technique. Of course, the magnitude of the oscillation and impact on the simulation depends on the scale of the reaction source term in the energy equation. If the scale of the reaction source term is small (i.e., nearly zero heat of reaction), the feedback is negligible and has almost no impact on the simulation. The obvious and rigorous fix to this problem is to limit the maximum fluid time step to something close to the reaction time scale (something on the order of 0.1−1 μs). This will ensure that the reaction energy source to the energy equation is sufficiently small. However, reducing the fluid time step is not desirable, as the computational cost will rise in proportion. A more tactful approach may be to apply a first order, low-pass, infinite impulse response (IIR) filter within the heating rate calculation. For discrete systems, the first order IIR filter can be written as49 F(n) = θF(n − 1) + (1 − θ )f (n)

■

EXAMPLE SURROGATE MODEL For this example, Powder River Basin (PRB) coal was considered. The ultimate analysis was 70.89% C, 4.50% H, 22.81% O, 1.30% N, and 0.51% S on a dry ash-free basis. The proximate analysis was 45.64% fixed carbon, 44.63% volatile matter, 6.38% ash, and 3.34% moisture. PCCL was used as the pyrolysis submodel in C3M and predicts gas yields for CH4, CO, CO2, H2, H2O, H2S, C2H6, C3H6, C2H4, HCN, and tar. This gives a total of 11 elementary reactions to add to the system according to eq 2 and a total of 14 closed form surface expressions that make up the surrogate model. The surrogate model was generated within C3M with the reactor temperature range set between 1000 and 1280 K and the heating rate specified between 25 and 5000 K/s. The distribution for each was chosen as uniform and the sampling method was Monte Carlo. The sampling size was set to ∼2000 data points. Figure 3 shows the design space filled with the sample points that were run in PCCL. Care should be taken by any user in choosing the design space. The operating conditions in the

(11)

Here, f is an input signal (in this case the calculated variable heating rate) and F is its output filtered value. Both input and output variables have finite values at discrete, unit steps, n. θ is a constant related to the frequency cutoff of the filter and will produce a stable filter when the value is between zero and one. To grasp the effect of θ on a variable value, it is instructive to consider a decay signal in which the initial value is unity and all following inputs are zero. It becomes clear then that the F(n) will take the following form: F(n) = θ n

(12)

The time evolved in the discrete system is simply the number of time steps multiplied by the size of the time step, Δt. It follows that θ can be defined in terms of some fractional decay value, Fd, the time step, and the decay time to that fractional decay value.

θ = FdΔt / t

(13)

It is up to the user to decide what length of time is acceptable to reach a given fractional decay value and is a problem specific issue. As a starting point, it is probably acceptable to reach an Fd of 0.1 in the time scale that it would take for fluid to travel through a single cell. For example, if the simulation runs at relatively stable time step of 1 ms, the cell size is about 1 cm, and the maximum expected velocity is 1 m/s, Δt/t would be 0.1 and θ would be ∼0.794. Of course, this value may have to be adjusted once the simulation starts. It is ideal to do the minimum amount of filtering necessary to prevent temperature oscillation in field of the solution. Becsause the oscillation is a function of the sum heat of reaction for all reactions, the necessary filtering may differ depending on the reaction system. If the effects of filtering are proven to be large, the only

Figure 3. Design space created in C3M. 7783



Article

Figure 4. Example response surface from C3M for the total yield of volatile gases.

Table 1. Surface Fit Parameters for the First Seven Response Variables ATot A B C D E F G H I J R2

2.418 −2.523 1.278 1.661 7.774 −5.046

× × × × × ×

EaTot 10−02 1000 1000 10−04 10−01 1001

0.9792

−2.621 × −7.186 × −4.429 × 3.413 × 7.968 × 4.275 × −2.126 × 1.136 × −1.178 × 3.740 × 0.9818

YTot 1001 10−04 10−07 10−11 10−11 10−06 10−09 10−01 10−04 10−08

−2.907 × −1.731 × 8.201 × −9.068 × −6.297 × 2.356 × −7.719 × 7.652 × −5.872 × 1.531 × 0.9954

YCO 1002 10−02 10−07 10−12 10−10 10−05 10−09 10−01 10−04 10−07

−2.808 × 1.081 × 1.834 × −3.526 × 1.470 × −2.894 × 7.423 × 5.401 × −3.195 × −7.071 × 0.9965

YCO2 10−01 10−05 10−09 10−13 10−12 10−08 10−12 10−04 10−08 10−11

−6.561 × −2.359 × 1.937 × −2.965 × 9.503 × −9.919 × 2.128 × 5.223 × −4.584 × 1.314 × 0.9668

Y H2 10−03 10−06 10−09 10−13 10−13 10−09 10−12 10−04 10−07 10−10

8.362 × −8.481 × −1.503 × −1.890 × 2.973 × 1.643 × −8.365 × −2.233 × 1.961 × −5.579 × 0.9979

YH2O 10−01 10−06 10−10 10−14 10−13 10−08 10−12 10−03 10−06 10−10

3.493 × −3.531 × 3.430 × −3.281 × −3.242 × −1.275 × 6.009 × 3.636 × −3.172 × 8.881 × 0.9617

10−03 10−08 10−11 10−14 10−14 10−10 10−13 10−07 10−09 10−12

2.437 × 4.055 × 3.265 × −2.700 × 2.046 × −6.816 × −4.382 × −6.769 × 6.371 × −1.229 × 0.9984

10−02 10−06 10−09 10−13 10−13 10−08 10−12 10−04 10−07 10−11

Table 2. Surface Fit Parameters for the Second Seven Response Variables YTar A B C D E F G H I J R2

6.001 × 5.094 × −8.940 × 1.314 × −3.562 × 4.531 × −8.407 × −1.541 × −3.062 × 1.549 × 0.9891

YCH4 10−01 10−06 10−09 10−12 10−12 10−08 10−12 10−04 10−07 10−10

−1.852 × −2.975 × 8.558 × −1.795 × 8.346 × −1.238 × −1.442 × 7.684 × −6.870 × 2.027 × 0.9537

YC2H4 10−01 10−06 10−10 10−13 10−13 10−09 10−12 10−04 10−07 10−10

−2.359 × −3.350 × 4.702 × −4.971 × 2.672 × −2.358 × 9.189 × 1.525 × −1.354 × 3.927 × 0.9748

10−02 10−07 10−10 10−14 10−14 10−09 10−13 10−04 10−07 10−11

YC2H6 5.329 × 1.679 × 9.487 × −9.141 × 1.383 × −3.565 × 1.464 × 7.930 × −5.381 × 9.900 × 0.9540

simulation should fall inside the design space in order to ensure a valid solution. This must be done because one cannot be

10−03 10−06 10−11 10−15 10−14 10−09 10−12 10−06 10−09 10−13

YC3H6 −9.952 × −3.815 × 3.550 × −3.849 × 3.370 × −1.599 × 5.743 × 9.902 × −8.834 × 2.569 × 0.9719

10−03 10−07 10−10 10−14 10−14 10−09 10−13 10−05 10−08 10−11

YH2S 5.197 × 7.036 × 8.063 × −1.395 × 5.499 × −6.816 × 1.345 × 9.501 × 5.490 × −2.861 × 0.9413

YHCN 10−02 10−07 10−11 10−14 10−13 10−10 10−13 10−05 10−08 10−11

assured that the surrogate model is valid outside of the design space. Because of this, care should be taken to limit any 7784



Article

Figure 5. Response of H2S mass fraction to reactor temperature and heating rate. The blue markers are PCCL data points. The red and green surface is the surrogate model generated in C3M.

and also prevented the rate from becoming negative when the values were close to zero. However, if the range of heating rates and reactor temperatures were relatively narrow, the cubic expression would also fit the PCCL data for ATot.

calculated reactor temperatures and heating rates to the maximum and minimums of the design space before they are used in any surrogate model calculations. The limiting method applied, or something similar, is necessary to maintain stability in the CFD code to ensure that unexpected values do not cause numerical instabilities. PCCL requires five other input variables to run the devolatilization besides heating rate and temperature. These variables are initial solids temperature, oxygen percentage, total run time, particle diameter, and pressure. The variables were set to 25C, 0%, 4s, 100 μm, and 1.55 MPa, respectively. It should be noted that any of the seven variables can be set as application inputs to the UQ engine within C3M, which will create a highdimensional design space. However, C3M cannot currently produce fits with any more than two independent parameters. High-dimension fitting will be a focus of future work and publication. That said, a design space that covers reactor temperature and heating rate covers the most significant variability and uncertainty in the reaction system and still demonstrates the viability of using surrogate models in CFD. Figure 4 shows an example response surface produced within C3M. The example is the response of PCCL’s prediction of the total yield of volatile gases. The blue stars are data points from PCCL and the red mesh is the response surface. The parameters for the closed form surface equations were determined by least-squares regression and are given in Tables 1 and 2. With the exception of ATot, all of the response surfaces were cubic with the form given in eq 14

F″(HR , Tr) = |A(Hr ) + B|C *|D(Tr) + E|F

Here, F′ and F″ are the response surface variables in the first row of Table 1 and Table 2. A through J correspond to the values in each row of Tables1 and 2. Hr and Tr are the heating rate in Kelvin per second and the reactor temperature in Kelvin, respectively. YTot is the total volatile yield and can be used to calculate the final char yield as in eq 16 XFinal =

Yi ,Tot 100

(16)

The correlation coefficients were quite high indicating a good fits to the PCCL output. However, several fits were lower than 0.98, which is somewhat less than desirable. Upon investigation, it was found that the lower correlation coefficients were not due to a poor surrogate model form but were due to “stair stepping” and scatter in the PCCL data. An example of this is shown in Figure 5. The stair stepping is likely not a real physical trend and appears to be linked to some rounding errors in the numerical techniques within PCCL. However, the smooth surface appears to adequately capture the major physical trends even if it does not precisely capture the direct PCCL output. PCCL does predict the tar composition, and it does vary slightly in response to heating rate and temperature. However, the variations are only a few tenths of a percent. For simplicity, the median tar composition was used and the stoichiometry was calculated from that. Equation 17 gives the tar reaction stoichiometry used in the reaction system

F ′(HR , Tr) = A + B(Hr ) + C(Hr )2 + D(Hr )3 + E(Hr )2 (Tr) + F(Hr )(Tr) + G(Hr )(Tr)2 + H(Tr) + I(Tr)2 + J(Tr)3

(15)

(14)

ATot was fit with a power expression of the form in eq 15. The power expression was used for ATot because it fit the model with much wider circumstance as compared to the cubic expression

15.6386C + 11.8368H + 2.66645O + 0.2466N + 0.0308S → tar 7785

(17)



Article

Mason and Gandhi equation for estimating the heating value of the fuel and standard equations for the heat of combustion. The Mason and Gandhi equation is less well known than the Dulong, Boie, Grummel and Davis, or Mott and Spooner equations for the calculation of the heating value of fuels.51−56 However, Mason and Gadhi developed a new formula by fitting experimental heating value data from 775 coals across all coal ranks.51,55 They demonstrated a substantial improvement in predictive capability when compared to the previous models.51 Translating the Mason and Gandhi equation to SI units, the dry ash laden higher heating value for a fuel can be estimated as

The molecular weight of tar was set to 246.6. Admittedly, there are many more significant figures in eq 17 than what normal scientific practice would allow. However, they are necessary to maintain the mass balance tolerances within most CFD codes, and without them, the code will produce errors and fail. Special Treatment of Tar. During the testing and development of this method, it was found that it was necessary to split eq 17 into five separate equations, one for each element. This was necessary because very small mass imbalances of minor species in the definition of either the coal or the stoichiometry of eq 17 could cause large product gas composition errors. In particular, this was problematic with sulfur in the case of PRB because if the coal was defined with slightly less sulfur than it was supposed to have or if the stoichiometry in eq 17 was slightly too high, it would cause the tar yield to stop early which would leave significant amounts of carbon, hydrogen, and oxygen in the coal/char for other competing gas products to consume. Splitting the tar reaction into five parts allows minor species to reach their full yield without affecting the major species. In this example, the reaction can be split up split up as follows: 5 × 15.6386C → TarC

(18)

5 × 11.8368H → TarH

(19)

5 × 2.66645O → TarO

(20)

5 × 0.2466N → TarN

(21)

5 × 0.0308S → TarS

(22)

HHVCHONS = 34094.51C̃ + 132298.2H̃ − 11985.9(Õ + Ñ ) + 6838.44S ̃ − 1530.51Ã

(23)

C̃ , H̃ , Õ , Ñ , S̃, and Ã are the mass fractions of the elements and ash in the coal on a dry ash laden basis. The reported standard deviation for this equation was 300.05 kJ/kg (∼1% of HHV) across all ranks of coal tested. The standard deviation for the Grummel and Davis equation and Mott and Spooner equation was about 1.4 times larger. The Boie and Dulong equations were in excess of 1.7 times as large.51 The heat of formation for the fuel can then be calculated according to ΔHf_CHONS = −HHVCHONS + +

The reaction rate for each chemical equation (eq 18 to 22) was that of the tar in eq 8 divided by five. It was necessary to define the system this way so that both the same mass of tar was formed as well as the same number of moles of gas (volume of gas) was generated per unit reaction. This means that the effective molecular weight of the tar in the system was 1233 g/ mol (5 × 246.6 g/mol), which is perfectly acceptable in the computational model. The other technical detail necessary for this method to work well is to ensure that the diffusion coefficients and the discretization method/order for the five tar species are exactly the same so that they advect together at all times once they are created in the system. In this way, these five species together can represent the single tar species defined in eq 17. Thermodynamic Property Estimation for CHONS Compounds. The elements in each CHONS compound need special treatment to derive their thermochemical properties such that when combined, they calculate the correct values for the single compound they represent. All other compound properties can be obtained from a standard thermochemical database, such as the Burcat database.50 With the heat of formation and specific heat for every compound involved in the pyrolysis reaction, a simple energy balance can be done for all reactions represented in eq 2. The pyrolysis reaction energy is then the sum of all the reaction energies weighted by their respective mass fraction. The energy balance calculations are handled within almost every CFD software and the only relevant information to this discussion is the heat of formation and specific heat for nonpure compounds such as the coal/char and the tar. Heat of Formation. The heat of formation for each of the CHONS species was estimated using a combination of the

1000C ΔĤ f_CO2(g) A WC

1000S 1000H ΔĤ f_SO2(g) + ΔĤ f_H2O(l) A WS 2A WH (24)

Equation 24 is the energy balance for a combustion reaction where all carbon goes to CO2, all sulfur goes to SO2, all nitrogen goes to N2, and all hydrogen goes to liquid water. Here, ΔHf_CHONS is the heat of formation of the fuel in kilojoules per kilogram, HHVCHONS is the higher heating value of the fuel in kilojoules per kilogram, ΔĤ f is the standard heat of formation for the pure substance indicated by the trailing subscript in kilojoules per mole, and AW is the atomic weight of the atom in the trailing subscript in grams per mole. Substituting eq 23 in eq 24 yields an equation in which the heat of formation is a function of a constant multiplied by the mass fraction of coal CHONS on a dry ash laden basis ⎛ 1000ΔĤ ⎞ f_CO2 ΔHf_CHONS = ⎜⎜ − 34094.51⎟⎟C̃ A WC ⎝ ⎠ ⎛ 1000ΔĤ ⎞ f_H2O + ⎜⎜ − 132298.2⎟⎟H̃ 2A WH ⎝ ⎠ + 11985.9(Õ + Ñ ) ⎛ 1000ΔĤ ⎞ f_SO2 + ⎜⎜ − 6838.44⎟⎟S ̃ A WS ⎝ ⎠ − 1530.51Ã

(25)

The conversion of any of the ultimate and proximate analysis values to dry ash laden values is as follows: 7786



Article

Table 3. Burcat Database Cp/R Values for CHONS Solids temperature range (K)

A

B

C

D

E

0−1000 1000−6000

6.18202141 × 10−12 −4.51413728 × 10−15

- 1.58440183 × 10−08 7.46618648 × 10−11

1.22151875 × 10−05 −4.54119512 × 10−07

−5.22309144 × 10−04 1.22115543 × 10−03

4.96951836 × 10−01 1.71771804 × 1000

Ux̃ =

Ux 1+A

Ã =

A 1+A

Burcat format, eq 30 multiplied by AW/R and fit with two fourthorder polynomials that were split at a temperature of 1000 K.50,57 Table 3 lists the parameter values that can be used in MFIX. The R2 value for the low and high temperature regions were 0.9999 and 0.9945, respectively. Again, the number of significant figures is far more than can be supported by experimental measure, but the Burcat has room for eight trailing zeros and there are no issues with keeping them within a numeric simulation. The equation to calculate Cp/R is from the values in Table 3 is

(26)

Here, the tilde represents the ash laden fraction, and Ux is the element in the ultimate analysis (C, H, O, N, or S). By substitution and the fact that the ultimate analysis sums to unity, eq 25 can be rewritten as ΔHf_CHONS =

⎞ ⎛ 1 ⎜ 1000ΔĤ f_CO2 − 34094.51 − 1530.51A ⎟⎟C ⎜ 1 + A⎝ A WC ⎠ +

⎞ ⎛ 1 ⎜ 1000ΔĤ f_H2O − 132298.2 − 1530.51A⎟⎟H ⎜ 1 + A⎝ 2A WH ⎠

Cp R

1 (11985.9 − 1530.51A)(O + N) + 1+A

(27)

Here, C, H, O, N, and S are the values from the ultimate analysis and A is the value of ash from the proximate analysis. Substitution of the heat of formation for the pure gases, the atomic weights, and the proximate ash for the PRB coal yields Hf_PRB = −62968.7C − 191622H + 11175.37(O + N) − 15244.9S

Cp R

(28)

Hf_tar = −66887.8C − 203748H + 11985.88(O + N) − 16119.7S

(29)

Because the heat of formation calculation in most CFD codes follows an ideal mixing law (ie: simple linear combination), the coefficients in eqs 28 and 29 are the heat of formation for each element in the CHONS compound. The remaining step is to work the coefficients into a form that each CFD code will recognize. In the case of MFIX, the code uses the Burcat format, so the coefficients have to be converted from kilojoules per kilogram to joules per mole and then divided by the gas constant.50,57 Specific Heat. For the solid CHONS compounds, the two characteristic temperature Einstein model proposed by Merrick was used.58 This model has shown very good agreement with a variety of coal and coal-like compounds. The two characteristic temperatures represent the characteristic oscillation frequency of the atoms in the graphitic-like plane and that normal to it.58 The two-parameter Einstein model is succinctly summarized in two equations ⎛ 1800 ⎞⎤ R ⎡ ⎜⎛ 380 ⎟⎞ ⎟ Cp = + 2f ⎜ ⎢⎣f ⎝ ⎠ ⎝ T ⎠⎥⎦ AW T f (x ) =

Cp = 4.22 × 10−3T

(33)

(34)

Here, the units for specific heat are joules per gram per Kelvin. Li and Eisermann both cite Hyman and Kay as the source of this information.60,63,64 However, after careful treatment of units, it was found that there was a mistake in the understanding of the originally reported units. Hyman and Kay report their units as Btu per pound. Using those units, one arrives at a different equation than that reported by Li and Eiserman as follows: Cp = 2.35 × 10−3T

(35)

If one were to assume that Hyman and Kay reported their units as calories per gram rather than Btu per pound, one obtains eq 34. Even if corrected, the data was taken for the condensed phase rather than the vapor phase and cannot be used in this case. In order to make a more accurate estimate of the specific heat of tar, standard relations used in petroleum processing were applied. Kesler and Lee report a correlation between molecular weight, mean-average boiling point, and specific gravity65,66

(30)

exp(x) ((exp(x) − 1)/x)2

= 8.57253933 + 8.45802261 × 10−3T

The literature is lacking in property measurements for coal- or biomass-derived tars in the gas phase. Several measurements of tars and pitches were made in the condensed range in the 1940s and 1950s.59−62 However, it is clear that the original research was conducted on condensed phase samples at relatively low temperatures.60 Furthermore, there has been a propagation of reporting errors in the more recent literature due to a unit conversion error. Both Li and Eisermann report the specific heat of tar as follows after making an assumption that the specific gravity of the tar was 1.17

Setting the ash to zero gives the equation for the tar Δ

(32)

Because the dependence on R and AW is removed in the transformation of eq 30 to the Burcat format, every element in the solid CHONS compound has the same polynomial coefficients. The ash specific heat was obtained from Merrick’s work as well.58 He reports a linear equation with respect to solid temperature in Celsius. Transforming the equation to temperature units of Kelvin, assuming a molecular weight of 120 g/mol, and transforming to the Burcat format, one obtains

⎛ ⎞ 1 ⎜ 1000ΔĤ f_SO2 ⎟S 6838.44 1530.51 A + − − ⎟ 1 + A ⎜⎝ A WS ⎠

Δ

= AT 4 + BT 3 + CT 2 + DT + E

(31)

Here, the Cp is in joules per gram per Kelvin and the T is the temperature of the solid in Kelvin. To put the model in the 7787



Article

of less than 0.8. For reduced boiling temperatures higher than 0.8, Lee and Kesler suggest using the following equation:66

MW = −12272.6 + 9486.4SG + (8.3741 − 5.99175SG)Tb ⎛ 222.466 ⎞ 107 + (1 − 0.77084SG − 0.02058SG2)⎜0.7465 − ⎟ Tb ⎠ Tb ⎝

ω = −7.904 + 0.1352K w − 0.007465K w2 + 8.359Tbr

⎛ 17.335 ⎞ 1012 + (1 − 0.80882SG − 0.02226SG2)⎜0.3228 − ⎟ Tb ⎠ Tb3 ⎝

+

(36)

Tc = 189.8 + 450.6SG + (0.4244 + 0.1174SG)Tb 5

10 Tb

(37)

0.0566 ln(Pc) = 5.689 − SG ⎛ 4.1216 0.21343 ⎞ Tb ⎟ − ⎜0.43639 + + ⎝ SG SG2 ⎠ 103 2 ⎛ 1.182 0.15302 ⎞ Tb ⎟ + ⎜0.47579 + + ⎝ SG SG2 ⎠ 106 3 ⎛ 9.9099 ⎞ Tb ⎟ − ⎜2.4505 + 2 10 ⎝ SG ⎠ 10

(38)

A 0 = 0.11828K w − 1.41779

(43)

A1 = −(6.99724 − 8.69326K w + 0.27715K w2 )10−4

(44)

A 2 = −2.2582 × 10−6

(45)

B0 = 1.09223 − 2.48245ω

(46)

B1 = −(3.434 − 7.14ω)10−3

(47)

B2 = −(7.2661 − 9.2561ω)10−7

(48)

⎡ (12.8 − K w )(10 − K w ) ⎤2 C=⎢ ⎥ ⎣ ⎦ 10ω

(49)

T ≤ 1000K

+ 4.8227 × 10−3T − 2.2371 × 10−6T 2 Cp = 2.3400 T > 1000K

(39)

(50)

Here, the specific heat is in joules per gram per Kelvin and temperature is in Kelvin. Figure 6 shows eq 50 graphically. Typical specific heats for hydrocarbons range from ∼1 J/g/K at about 300 K to around 2.5−3 J/g/K at 1000 K.70,71 Therefore, this approximation method appears to produce reasonable estimations for gas phase thermodynamic properties of tar. All the necessary chemical and thermodynamic equations have been defined for the chemistry model and it can now be implemented within a CFD code.

6.09648 6 + 1.28862ln(Tbr) − 0.169347T br Tbr 15.6875 6 − 13.4721ln(Tbr) + 0.43577T br Tbr

ln(Pbr) − 5.92714 + 15.2518 −

(42)

Cp = −2.4558 × 10−1

For the tar mixture, the Watson factor was calculated to be ∼9.30. According to Watson, characterization factors between 12.5 and 13 represent purely paraffinic mixtures and values below 10 represent aromatic mixtures.69 This is fitting as the molar hydrogen to carbon ratio for the tar is ∼0.757, which indicates that it has a relatively high degree of aromatic character. If the H:C ratio was closer to 2 or greater, it would indicate a paraffinic character. The acentric factor can now be estimated according to Lee and Kesler66,67 ω=

Cp = [A 0 + A1T + A 2 T 2 − C(B0 + B1T + B2 T 2)]

Here, Cp is in J/g/K. The equation is generally valid up to about 1200C according to Chang. However, eq 42 is an inverted parabola with a maximum between 1000K and 1200K. Most substances do not decrease their specific heat with temperature because the degrees of freedom available to the substance only increase with temperature. To approximate higher temperature behavior, the value of the specific heat above 1000K was assumed to be that calculated by eqs 42 to 49 at 1000K. Substituting the values calculated in eqs 39 and 40, one obtains the function that can be used in the CFD simulation.

Using the assumed specific gravity of 1.17 and the calculated mean-average boiling point of 715 K, one obtains a critical temperature of ∼974 K and a critical pressure of 23.58 bar. The Watson factor, Kw, can now be calculated for the hydrocarbon mixture. The Watson factor is used to classify the type of hydrocarbons in the mixture based on the contents of paraffinic, aromatic, and naphthenic compounds65,68,69 (1.8Tb)1/3 Kw = SG

(41)

In this case, the reduced boiling temperature was calculated to be ∼0.734, so eq 40 was used to calculate an acentric factor of ∼0.634. It is now possible to calculate the specific heat. Chang details the method used in Aspen Hysis to calculate the specific heat from the Watson factor and the acentric factor65

Here, MW is the molecular weight of the hydrocarbon mixture in grams per mole, SG is the specific gravity, and Tb is the meanaverage boiling point in units of Kelvin. In this case, the molecular weight was reported by PCCL as 246.6 g/mol and a specific gravity of 1.17 was assumed based on previously reported literature.60,62−64 Numerically solving for Tb, one obtains a value of ∼715 K, which is reasonable. Lee and Kesler also report correlations for critical temperature and pressure based on specific gravity and mean-average boiling point66,67

+ (0.1441 − 1.0069SG)

(1.408 − 0.01063K w ) Tbr

■

APPLICATION IN CFD The surrogate chemistry model was implemented in MFIX (Eulerian−Eulerian) to compare the exit gas composition with that predicted by the chemistry surrogate model. A simple pseudo 2D transport gasifier was chosen as a test case because it

(40)

Here, Pbr is the reduced boiling point pressure (atmospheric pressure divided by the critical pressure) and Tbr is the reduced boiling point temperature (Tb divided by the critical temperature). This equation is good for reduced boiling temperatures 7788



Article

Fresh coal was brought in through a side inlet at 298 K with a mass flow of 1.0 g/s. At the same inlet, nitrogen was brought in at a flow rate of 0.1 g/s and a temperature of 298 K. The gas volume fraction was set to 0.995. The inlet defined to be between 32 and 34 cm from the bottom of the reactor. The coal was assigned to the first solid phase with a particle diameter of 100 μm and a density of 2.85 g/cm3. Recycle ash was brought in near the bottom of the reactor and was assigned to a second solid phase with a particle diameter of 100 μm and a density of 2.85 g/cm3. These solids were brought in at a composition of 99.55% ash and 0.45% carbon by mass. The gas and solids came in at 1280 K. These solids were brought in at a flow rate of 100.0 g/s with nitrogen at 0.2 g/s at the same temperature. The gas volume fraction was set to 0.5. The inlet was located between 14 and 18 cm from the bottom of the reactor. The outlet of the system was set to be a pressure outlet near the top of the reactor. The outlet pressure was set to 1.55 MPa and was set to be between 290 and 294 cm from the bottom of the reactor. All other boundaries were set to be no-slip walls. Initial Conditions. The entire volume of the system was initialized to be 1280 K and 1.55 MPa. The composition of the system was 98% by volume of nitrogen and the remainder was ash from the second solid phase. There were no constituents from the first coal phase in the system initially. Filtering Heating Rate and Reactor Temperature. For this study, the cells at the outlet were used to calculate the reactor temperature used as input to the chemistry surrogate. The values were spatially averaged at each time step to obtain a single reactor temperature. Equation 11 was used with a variable cutoff frequency defined in eq 51 because MFIX utilizes a variable time step

Figure 6. Tar specific heat predicted by the present model.

runs quickly and was sufficient to demonstrate how the surrogate model concept could be applied successfully in CFD. Geometry and Mesh. The reactor consisted of a simple rectangular geometry that was 3 m tall and 10 cm wide (see Figure 7). The geometry was gridded at a ∼1 cm resolution (10 cells in the x direction and 320 in the y direction). The grid was set to be a single cell deep. Boundary Conditions. Hot nitrogen was brought into the system at 1280 K at the bottom of the reactor. The entire width of the bottom was made to be an inlet. The inlet was set to be a mass flow inlet with a value of 2.76 g/s.

θ = 1 − 0.2 × Δt

(51)

Here, Δt is the time step in seconds. It should also be noted that the maximum time step was limited to 1 ms which ensured that the value of θ was always between 0 and 1, which meant that the filter was stable. Further, this gave the filter a half-life of about ∼3.5 s. It was beneficial to apply a slow filter because the desire is to have the reaction system respond to some representative average temperature reactor temperature rather than the highly variable instantaneous solids exit temperature. An IIR filter was used for filtering the heating rate, but the θ was set to be a constant and not variable as in eq 51, which was used for the filtering of the reactor temperature. The sensitivity of this constant was investigated and discussed in further detail below. Remaining CFD Parameters. All other transport, momentum, and energy parameters and equations were set to the default values in MFIX that can be found in the MFIX manual.72 The second public release of MFIX in 2013 was used. Copies of the mfix.dat, usr1.f, and usr_rates.f files used in these simulations can be obtained through correspondence with the lead author.

■

CONSTANT HEATING RATE AND REACTOR TEMPERATURE SIMULATION The first test performed was one in which the heating rate and reactor temperatures were artificially forced to constant at values of 1000 K/s and 1100 K respectively in the chemistry submodel. A successful test would show a constant normalized product gas composition in which the mass fraction of the gas species was

Figure 7. Reactor configuration diagram. 7789



Article

Table 4. Product Gas Yields from the Surrogate Model (SM) and the 30 to 45 Second Time Averaged CFD Outlet (CFD) SM CFD

YCO

YCO2

Y H2

YH2O

YTar

YCH4

YC2H4

YC2H6

Y C3 H 6

YH2S

YHCN

0.171 58 0.171 58

0.180 06 0.180 05

0.009 34 0.009 34

0.161 74 0.161 73

0.299 52 0.299 51

0.094 13 0.094 13

0.031 29 0.031 29

0.008 48 0.008 48

0.025 17 0.025 17

0.008 69 0.008 69

0.010 01 0.010 01

identically defined by eq 14 and the constants in Tables 1 and 2. The results of this calculation are shown in Table 4 as well as the outlet normalized product distribution from the CFD averaged over the last 15 s. Figure 8 shows the outlet, normalized product gas distribution as a function of time.

Figure 9. Tar composition at the outlet of the reactor (constant heating rate and reactor temperature).

Figure 8. Time series of the normalized product gas distribution (constant heating rate and reactor temperature).

The product gas distribution is indeed the same as the surrogate model. The averaged gas composition matches the surrogate model to five decimal places in most cases and to four decimal places in every case. Further, the product gas composition is constant when normalized to include only the product gases. The initial spike is a numerical aberration that occurs when the total product gas composition is near the lower simulation limit and should be ignored. The tar mass fraction is also consistent with what was defined in eqs 18 to 22. The tar mass fraction used to calculate the stoichiometry was 0.761C, 0.048H, 0.173O, 0.014N, and 0.004S. The time averaged CFD results matched those compositions to five decimal places. Figure 9 shows the elemental tar composition at the outlet of the reactor as a function of time. Again, it is constant as expected. Furthermore, the fact that both the total mass fraction of tar within the overall product gas distribution as well as the elemental composition of the tar match the surrogate model predictions indicate that the special treatment of tar did indeed produce the same mass as if eq 17 had been used rather than eqs 18 to 22. Figure 10 shows the char yield (equivalent to the fractional gas yield) at the exit of the gasifier in the CFD simulation and that predicted by the surrogate model. No solids make it to the exit of the gasifier until ∼3.75 s. The equations were set up such that if there were no solids in a given control volume, the char yield was set to the predicted value so that the calculated rate would be zero. The messy spike at about five seconds is likely the result of the chemistry solver overstepping the reaction due to very low solids concentration. Appreciable amounts of solids make it to the exit at about 6.14 s. From then on, the char yield matches the predicted values.

Figure 10. Char yield at the exit of the reactor (constant heating rate and reactor temperature).

The slow transition from no solids to some appreciable concentration of solids is also reflected in the first phase solids exit temperature shown in Figure 11. The temperature transitions from the initial condition of 1280 K to the pseudosteady-state temperature of ∼1160 K in the same time region as the spike in Figure 10. Furthermore, the significant variation in temperature about the mean is indicative of quite unsteady hydrodynamics that occur in the transport gasifiers and supports the need to filter the solids exit temperature to obtain estimates of the reactor temperature. In this case, the reactor temperature was artificially set to a constant value, which was about 60 K lower than the average solids temperature at the exit of the reactor at pseudo-steady-state. This test gives support to the notion that the surrogate model has been implemented correctly in the CFD model because the CFD model is able reproduce the predictions of the surrogate model to a substantial level of accuracy. Furthermore, this simulation will produce results that are very similar to the results 7790



Article

Figure 11. Char temperature at the exit of the reactor and the reactor temperature (constant heating rate and reactor temperature).

Figure 12. Time series of the normalized product gas distribution (variable heating rate and reactor temperature).

produced without the surrogate model. That is to say the results should not be significantly different from a simulation in which a single, fixed chemistry pyrolysis reaction had been produced from a single PCCL run in C3M with the same fixed heating rate and reactor temperature. The only possible way the results can differ is if the response values from the singular PCCL run do not lie directly on all of the response surfaces in the surrogate model discussed here. The failure of the surrogate model to exactly reproduce any particular PCCL run can occur as a result of at least two reasons: First, the response surfaces are suboptimal fits to the response variables in PCCL data. Second, the response variables from PCCL themselves contain error (random or systematic). The first can be addressed by changing the response surface expression to suit the data form. The second cannot be addressed and represents an inherent error in the pyrolysis model. That said, if the response surface adequately represents the response variable from the pyrolysis model, the surface should represent the median value in an error distribution at any point in the domain. Furthermore, it is possible to associate a probability of obtaining that median value and, thus, an uncertainty related to the surrogate model itself relative to the pyrolysis model output. Therefore, there is value in using the surrogate model methodology even if the heating rate and reactor temperature are held constant because one would have knowledge about the uncertainty in the predicted values.

time ∼2.6 s, the first product gases make it to the exit of the gasifier. At this time, the tar mass fraction was ∼28%. At ∼5.8 s, the tar mass fraction dropped to ∼24%. Following this, the tar composition slowly rose to ∼27% over the next 15 s. The gas that makes it out of the reactor first was created at the start of the coal injection. Under these conditions, the heating rate and reactor temperature would be their highest. These conditions would be true because the coal that is injected first experiences the greatest convective heat transfer. The comparatively high heat transfer is due to the fact that as the injection of coal continues, the pathway traveled by the coal becomes comparatively cooler because the coal ahead in the path has cooled the surroundings and, thus, reduced the rate of convection experienced by subsequent solids. The first coal that is injected does not experience the influence of any solids coming before it and the heating rate remains high. At the same time, the cells at the exit experience the relatively high initial temperature of 1280 K because nothing cooler has made it to the exit. The cooler pathway is established within a couple of seconds of injection, which is faster than the time it takes solids to get to the exit of the reactor. Under these conditions, the heating rate is relatively low and the calculated reactor temperature remains high. At 5−6 s, solids begin to exit the reactor and the calculated reactor temperature drops over the next 15 s as seen in Figure 13. Understanding the conditions under which the tar was generated allows for understanding of why the tar composition changed. It was observed that the generation conditions changed in succession: (1) high heating rate, high reactor temperature, (2) low heating rate, high reactor temperature, and (3) low heating rate, low reactor temperature. The corresponding changes in tar yield can be seen in Figure 14, which demonstrates consistency between the observed reaction conditions and the tar product yields. Further evidence that the changing conditions are producing reasonable results is found by examining the profiles of product gases that exhibit strong correlations to either heating rate or reactor temperature but not both. Carbon dioxide, water, and methane all exhibit strong negative correlations with heating rate (−0.859, −0.875, and −0.810, respectively) and relatively weak negative correlation with reactor temperature (−0.274, −0.270, and −0.345, respectively). On the basis of observations of the transient reaction conditions, the expectation was the gases

■

VARIABLE HEATING RATE AND REACTOR TEMPERATURE SIMULATION In this simulation, the heating rate in each cell was estimated using eq 10 with the IIR filter in eq 11 with the filter parameter set to 0.98. The reactor temperature was calculated by averaging the phase one solids temperature at the cells at the exit of the gasifier and putting them through the IIR filter in eq 11 with the filter parameter calculated in eq 51. In this case, the normalized product gas composition is not expected to be constant and should shift in response to the influence of heating rate and reactor temperature. These shifts are evident in Figure 12, which shows the normalized product gas distribution as a function of time. Figure 12 shows that the gas composition is not constant and that the most widely varying species fraction is tar. Taking tar as the example to understand these changes, it is observed that at a 7791



Article

Figure 15. Char yield at the exit of the reactor (variable heating rate and reactor temperature). “CFD” is the instantaneous char yield calculated at the exit of the gasifier. “SM” is the maximum char yield calculated by the surrogate model.

Figure 13. Char temperature at the exit of the reactor and the reactor temperature (variable heating rate and reactor temperature).

should quickly rise in concentration between 3 and 5 s and then hold their value in Figure 12, which is what was observed. Hydrogen production exhibits almost no correlation with heating rate (−0.064) and a near perfect positive correlation with reactor temperature (0.994). If the understanding of the reaction conditions is correct, hydrogen should exhibit no change in gas yield under 5 s and should exhibit a gradual decrease over the following 15−20 s, which is what was observed. The maximum char yield also has a weak correlation to heating rate (−0.127) and a strong positive correlation to reactor temperature (0.969). The result is a slow drop in the maximum char yield as the calculated reactor temperature drops as seen in Figure 15.

Figure 15 shows the char yield at the exit of the gasifier (CFD) and the maximum char yield calculated by the surrogate model at the exit. It is reasonable to think that the instantaneous char yield should not exceed the maximum char yield in Figure 15. Under most circumstances, it should not. However, if the solids concentration is very low somewhere in the field of the solution, it is possible for the chemistry solver to overstep the yield and still meet the convergence criteria. Examination of the solution between 23.4 and 23.8 s reveals that an upset condition caused solids to enter and mix into the headspace above the exit of the reactor which caused the solids concentration in the headspace to rise above the dilute solids limit (but only slightly), which created a region near the exit in which the solids were

Figure 14. Response surface of tar yield. 7792



Article

overconverted due to the chemistry solver overstepping its solution. This region was then flushed out of the reactor before any mixing could occur to dilute the overconverted solids. Clearly, the issue is not in the definition of the chemistry but in the methods used to solve the chemical equations. Such overpredictions can be understood as numerical falsities that have a minor impact on the solution as a whole. It appears that these falsities appear as high frequency spikes in the time series data. Furthermore, because these events are isolated to areas where there is extremely low solids concentration, the impact on mass balance is very minimal as the amount of mass being converted is very low. It should also be noted that the lower frequency and high amplitude variations of the measurable quantities at the exit indicates that the system is poorly mixed, which is consistent with the hydrodynamic nature of the transport reaction systems. The tar composition was identical to the previous test, which is expected as the defined ultimate analysis was the same and there was no influence of heating rate or reactor temperature in the mass balance equations. Perhaps the most important feature of this method is that the char composition is allowed to vary and consist of more elements than simply carbon as was assumed in previous coal models.2,3,73,74 This is quite important because if char is only fixed carbon, then all the other elements in coal must be artificially pushed to a pseudospecies. Because all product gases in the pyrolysis reactions were defined as pure gases except for pseudospecies (like tar), the pseudospecies ended up containing excessively high amounts of minor elements by atomic balance. This is especially true for lower rank coals like PRB because they contain high amounts of minor elements. In many cases, the elemental balance would produce tars with hydrogen-to-carbon ratios far exceeding four, which is impossible as there can be no more than four hydrogens attached to one carbon. In most reasonable accounts, the hydrogen to carbon ratio should be less than two as this represents fully saturated hydrocarbons, and it is known that coal derived heavy hydrocarbons contain significant aromatic structure which would further lower the hydrogen to carbon ratio. This false atomic balance makes estimation of thermodynamic properties from chemical composition impossible and one must assume that tar behaves like a (set of) representative substance(s). In the surrogate methodology demonstrated in this work, char can exist as a mixture of atomic species, which eliminates the problems discussed above. Figure 16 shows the mass fractions of char exiting the reactor as a function of time. The unsteady composition of the tar is indicative of both the unsteady hydrodynamics as well as the changing reaction conditions experienced by coal/char in the system. The time-averaged field values of gas yield, heating rate and pyrolysis rate over the last 20 simulation seconds are shown in Figure 17. The results are sensible. The highest average heating and reaction rates as well as the lowest gas yield are near the injection site. The gas yield increases and the reaction rate decreases as the solids move up the riser, which is expected. Of particular note is the relatively low averaged heating rate in the area around the injection site as compared with the edges of the solid phase plume seen on initial injection of the solids (see Figure 19, right side), which also supports the trends seen in the initial exit gas composition. The heating rate on the plume edge is consistently higher (in the low thousands of Kelvin per second) than that following the plume and also lower than the

Figure 16. Char composition at the exit of reactor (variable heating rate and reactor temperature).

Figure 17. Time averaged field variables.

time averaged results (in the low to mid hundreds of Kelvin per second), which is consistent with the analysis above. This test demonstrates that the chemical surrogate model has been implemented correctly as it produces results that are sensible and within the bounds of the surrogate model.

■

EFFECT OF FILTERING HEATING RATE In this test, the filtering constant, θ, was varied in eq 11 to assess the impact on the simulation. The impact was measured by analyzing the time averaged output gas concentration and by visually inspecting the heating rate field values. It is conceptually simple to discuss the appropriate filtering constant in terms of the inverse filter speed (the exponent in eq 12), which is the number of discrete steps needed to get the output filter value to be equal to Fd in a unit decay case. Doubling n makes the filter twice as slow. It is thus possible to set the speed of the filter in a conceptual sense. For example, if 7793



Article

inverse filter speed, the heating rate fields will appear to be generally smooth and continuous. An example of each case is shown in Figure 19. As a rule of thumb, the inverse filter speed should be increased until a smooth heating rate field is observed and then doubled to ensure that the simulation is operating above the critical inverse filter speed.

one were to set the speed of the filter to reach 90% of its final output value in 10 discrete steps (time steps in the case of CFD), one should set the filter constant to be ∼0.794. The inverse filter speed was varied between ∼10 to ∼180 (filter constants of 0.8 to 0.9875 with Fd set to 0.1). The time averaged output gas concentration as a function of the filter speed (n) is shown in Figure 18.

■

CONCLUSIONS AND FUTURE WORK

The method outlined above represents one way to comprehensively represent the fidelity of a complicated pyrolysis submodel within CFD without the need to conduct iterative CFD simulations or make in-simulation application calls. The method relies on the ability to represent the pyrolysis submodel through a series of response surfaces that together form a comprehensive surrogate model. The surrogate model is designed to give the CFD code the ability to adapt and respond to changing reaction conditions on the fly without the need to predefine a thermal pathway in which the pyrolysis is supposed to occur. In doing so, the method removes much of the guessand-check methodology previously used. A fair criticism of this method is that filtering key variables may be necessary and that the filtering constants cannot be determined before running a simulation. However, these constants can be adjusted during the progress of the simulation without the need to rebuild the chemistry model. With proper adjustment during the course of the simulation, the filtering impact should be minimal. It should also be noted that this method may not produce realistic values during transient operational shifts, such as reactor startup, but should produce reliable results for evaluating steady state or pseudosteady state transient simulations. The trouble during transient operational shifts can stem from the evaluation of final reactor temperature. If the evaluation domain is significantly downstream of important reactant inlets, the entire reaction system can lag behind the flow system until a new steady state is reached. If this is important in the simulation, perhaps a different method of evaluating reactor temperature should be developed. The set of tests discussed in this paper demonstrate that the surrogate model methodology has been implemented correctly and functions as expected. Under constant conditions, the computational model will produce the same output as the surrogate model predicts. Under varying conditions, the gas compositions, yield, and char elemental balance produce results which are reasonable given what is known about the surrogate models and the conditions in the gasifier. Further, it was demonstrated that the necessary filtering of the heating rate has little impact on the simulation as long as the inverse filtering constant is above a critical value. In short, the method has been verified to produce reasonable results. At the current time, this method is undergoing continued testing at NETL. Initial results are promising. The results of these tests will be published as parts of future works in this series. In addition, this method is being integrated into C3M and will be available to the public in the next release of C3M. C3M will be able to directly write the reaction equations described in this work within user defined functions for both MFIX and Fluent to minimize or eliminate user error. Furthermore, the general methodology is being applied to gasification submodels as higher-dimensionality response models are being developed. The hope of this body of work is to develop a completely predictive gasification model that minimizes the amount of educated guessing that the user must do and thereby minimizes error and improves model accuracy.

Figure 18. Time-averaged exit gas composition (30−45s) as a function of the inverse filter speed.

As shown, there is a critical inverse filter speed at about n = 25. Above this inverse filter speed, the simulation results appear to be relatively stable and constant. Below this inverse filter speed, the results appear to be much more unstable and less repeatable. The state of the inverse filter speed in the simulation relative to the critical inverse filter speed is easily observable in the continuity of the heating rate field as seen in Figure 19. If the inverse filter speed is too low, one will observe highly discontinuous, checkerboard-like characteristics in the heating rate fields. This kind of condition is obviously false and nonphysical. If the inverse filter speed is above the critical

Figure 19. Heating rate fields from n = 10 (left) and n = 45 (right) simulations at near the coal injection site just after the start of the coal injection. The colored cells are those which have solid fuel in them. 7794



■

Article

(14) Fletcher, T. H.; Kerstein, A. R.; Pugmire, R. J.; Solum, M. S.; Grant, D. M. Chemical Percolation Model for Devolatilization 0.3. Direct Use of C-13 Nmr Data to Predict Effects of Coal Type. Energy Fuels 1992, 6 (4), 414−431. (15) Solomon, P. R.; Hamblen, D. G.; Carangelo, R. M.; Serio, M. A.; Deshpande, G. V. General-Model of Coal Devolatilization. Energy Fuels 1988, 2 (4), 405−422. (16) Solomon, P. R.; Hamblen, D. G.; Carangelo, R. M.; Serio, M. A.; Deshpande, G. V. Models of Tar Formation during Coal Devolatilization. Combust. Flame 1988, 71 (2), 137−146. (17) Solomon, P. R.; Serio, M.A.; Hamblen, D.G.; Yu, Z.Z.; Charpenay, S. Advances in the Fg-Dvc Model of Coal Devolatilization. Abstr. Pap.Am. Chem. Soc. 1990 199. (18) Zhao, Y. X., Serio, M.A.; Solomon, P.R. Modeling the devolatilization of large coal particles. Coal Science, Vols I and II, 1995. 24: p 833−836. (19) Wojtowicz, M. The FG-DVC Coal Pyrolysis Model User’s Guide; Advanced Fuel Research: East Hartford, CT, 2005. (20) Niksa, S.; Liu, G. S.; Hurt, R. H. Coal conversion submodels for design applications at elevated pressures. Part I. devolatilization and char oxidation. Prog. Energy Combust. Sci. 2003, 29 (5), 425−477. (21) Liu, G. S.; Niksa, S. Coal conversion submodels for design applications at elevated pressures. Part II. Char gasification. Prog. Energy Combust. Sci. 2004, 30 (6), 679−717. (22) Chen, J. C.; Niksa, S. Coal Devolatilization during Rapid Transient Heating 0.1. Primary Devolatilization. Energy Fuels 1992, 6 (3), 254−264. (23) Chen, J. C.; Castagnoli, C.; Niksa, S. Coal Devolatilization during Rapid Transient Heating 0.2. Secondary Pyrolysis. Energy Fuels 1992, 6 (3), 264−271. (24) Niksa, S.; Kerstein, A. R. The Distributed-Energy Chain Model for Rapid Coal Devolatilization Kinetics 0.1. Formulation. Combust. Flame 1986, 66 (2), 95−109. (25) Niksa, S. The Distributed-Energy Chain Model for Rapid Coal Devolatilization Kinetics 0.2. Transient Weight-Loss Correlations. Combust. Flame 1986, 66 (2), 111−119. (26) Niksa, S.; Kerstein, A. R. Flashchain Theory for Rapid Coal Devolatilization Kinetics 0.1. Formulation. Energy Fuels 1991, 5 (5), 647−665. (27) Niksa, S. Flashchain Theory for Rapid Coal Devolatilization Kinetics 0.2. Impact of Operating-Conditions. Energy Fuels 1991, 5 (5), 665−673. (28) Niksa, S. Flashchain Theory for Rapid Coal Devolatilization Kinetics 0.3. Modeling the Behavior of Various Coals. Energy Fuels 1991, 5 (5), 673−683. (29) Niksa, S. Flashchain Theory for Rapid Coal Devolatilization Kinetics 0.4. Predicting Ultimate Yields from Ultimate Analyses Alone. Energy Fuels 1994, 8 (3), 659−670. (30) Niksa, S. Flashchain Theory for Rapid Coal Devolatilization Kinetics 0.5. Interpreting Rates of Devolatilization for Various Coal Types and Operating-Conditions. Energy Fuels 1994, 8 (3), 671−679. (31) Niksa, S. Flashchain Theory for Rapid Coal Devolatilization Kinetics 0.6. Predicting the Evolution of Fuel Nitrogen from Various Coals. Energy Fuels 1995, 9 (3), 467−478. (32) PCCL User Guide and Tutorial, version 4.1; Niksa Energy Associates: Belmont, CA, 2008. (33) Niksa, S.; Kerstein, A. R.; Fletcher, T. H. Predicting Devolatilization at Typical Coal Combustion Conditions with the Distributed-Energy Chain Model. Combust. Flame 1987, 69 (2), 221− 228. (34) Niksa, S. Predicting the distributions of volatile products from any coal with FLASHCHAIN (TM). Prospects for Coal Science in the 21st Century 1999, 657−660. (35) Niksa, S. Rapid Coal Devolatilization as an Equilibrium Flash Distillation. AIChE Journal 1988, 34 (5), 790−802. (36) Vascellari, M.; Arora, R.; Pollack, M.; Hasse, C. Simulation of entrained flow gasification with advanced coal conversion submodels. Part 1: Pyrolysis. Fuel 2013, 113, 654−669.

AUTHOR INFORMATION

Corresponding Author

*D. Van Essendelft. E-mail: [email protected] Notes

The authors declare no competing financial interest. Disclaimer. This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference therein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed therein do not necessarily state or reflect those of the United States Government or any agency thereof.

■

REFERENCES

(1) Fletcher, T. H.; Hardesty, H.D. Compilation of Sandia Coal Devolatilization Data: Milestone Report; Sandia National Laboratory: Albuquerque, NM, 2002. (2) Guenther, C., Shahnam, M.; Syamlal, M.; Longanbach, J.; Cicero, D.; , Smith, P.V. CFD Modeling of a Transport Gasifier. In Proc. - Annu. Int. Pittsburgh Coal Conf., 19th; International Pittsburgh Coal Conference: Pittsburgh, PA, 2002. (3) Guenther, C., Syamlal, M.; Longanbach, J.; Smith, P.V.; CFD Modeling of a Transport Gasifier Part II. In Proc. - Annu. Int. Pittsburgh Coal Conf., 20th; International Pittsburgh Coal Conference: Pittsburgh, PA, 2003. (4) Niksa, S.; Liu, G.-s.; Hurt, R. H. Coal conversion submodels for design applications at elevated pressures. Part I. devolatilization and char oxidation. Prog. Energy Combust. Sci. 2003, 29 (5), 425−477. (5) Pollack, M. Integration of a detailed Pyrolysis Submodel for use in CFD Simulations of Coal Gasification and Combustion. In Department of Mechanical, Process and Energy Engineering; Technische Universität Bergakademie Freiberg: Freiberg, Germany, 2012. (6) Solomon, P. R.; Hamblen, D. G.; Carangelo, R. M.; Serio, M. A.; Deshpande, G. V. General model of coal devolatilization. Energy Fuels 1988, 2 (4), 405−422. (7) Solomon, P. R.; Serio, M. A.; Despande, G. V.; Kroo, E. Crosslinking reactions during coal conversion. Energy Fuels 1990, 4 (1), 42− 54. (8) Solomon, P. R.; Serio, M. A.; Suuberg, E. M. Coal pyrolysis: Experiments, kinetic rates and mechanisms. Prog. Energy Combust. Sci. 1992, 18 (2), 133−220. (9) Stubington, J. F.; Sasongko, D. On the heating rate and volatile yield for coal particles injected into fluidised bed combustors. Fuel 1998, 77 (9−10), 1021−1025. (10) Tremel, A.; Spliethoff, H. Gasification kinetics during entrained flow gasification - Part I; Devolatilisation and char deactivation. Fuel 2013, 103, 663−671. (11) Gibbins-Matham, J.; Kandiyoti, R. The Effect of Heating Rate and Hold Time on Primary Coal Pyrolysis Product Distribution. Fuel 1987, 318−323. (12) Grant, D. M.; Pugmire, R. J.; Fletcher, T. H.; Kerstein, A. R. Chemical-Model of Coal Devolatilization Using Percolation Lattice Statistics. Energy Fuels 1989, 3 (2), 175−186. (13) Fletcher, T. H.; Kerstein, A. R.; Pugmire, R. J.; Grant, D. M. Chemical Percolation Model for Devolatilization 0.2. Temperature and Heating Rate Effects on Product Yields. Energy Fuels 1990, 4 (1), 54− 60. 7795



Article

(59) Chemistry of Coal Utilization: Second Supplimentary Volume; Elliott, M. A., Ed.;John Wiley & Sons, Inc.: New York, 1981. (60) Hyman, D.; Kay, W. B. Heat Capacity and Content of Tars and Pitches. Ind. Eng. Chem. 1949, 41 (8), 1764−1768. (61) Chemistry of Coal Utilization; Lowry, H. H., Ed.; John Wily & Sons, Inc.: New York, 1945. (62) Chemistry of Coal Utilization: Supplemetary Volume; Lowry, H. H., Ed.; John Wily & Sons, Inc.: New York, 1963. (63) Eisermann, W.; Johnson, P.; Conger, W. L. Estimating Thermodynamic Properties of Coal, Char, Tar and Ash. Fuel Process. Technol. 1980, 3 (1), 39−53. (64) Li, C.; Suzuki, K. Tar property, analysis, reforming mechanism and model for biomass gasificationAn overview. Renewable and Sustainable Energy Reviews 2009, 13 (3), 594−604. (65) Chang, A. F.; Pashikanti, K.; Liu, Y.A. Refinery Engineering: Integrated Process Modeling and Optimization; Wiley: New York, 2013. (66) Kesler, M. G.; Lee, B. I. Improve Prediction of Enthalpy of Fractions. Hydrocarbon Processing 1976, 55 (3), 153−158. (67) Lee, B. I.; Kesler, M. G. A generalized thermodynamic correlation based on three-parameter corresponding states. AIChE J. 1975, 21 (3), 510−527. (68) da Silva, J., Leon, J.; Arredono, H.; , and Molina, A.. On the Exergy Determination for Petroleum Fractions and Separation Process Efficiency. In Brazilian Congress of Thermal Sciences and Engineering, 14th ; ABCM: Rio de Janeiro, Brazil, 2012. (69) Watson, K. M.; Nelson, E. E.; Murphy, G. B. Characterization of Crude oil Fractions. Ind. Eng. Chem. 1935, 27 (12), 1460−1464. (70) Burgess, D. R. Thermochemical Data. In NIST Chemistry WebBook, NIST Standard Reference Database Number 69; Linstrom, P.J., Mallard, W.G., Eds.; National Institute of Standards and Technology: Gaithersburg MD, 2012. (71) CRC Handbook of Chemistry and Physics, 82nd ed.; Lide, D. R., Ed.; CRC Press LLC: Boca Raton, FL, 2001. (72) Multiphase Flow with Interphase eXchanges. https://mfix.netl. doe.gov/ (accessed Dec 2013). (73) Guenther, C., Syamlal, M.; Longanbach, J.; , and Smith, P.V. Two-Fluid Model of an Industrial Scale Transport Gasifier. In AIChE Annual Meeting 2003: San Francisco, CA, 2003. (74) Syamlal, M.; Bissett, L.A. METC Gasifier Advanced Simulation (MGAS) Model. National Technical Information Service: Springfield, VA, 1992.

(37) Hashimoto, N.; Kurose, R.; Shirai, H. Numerical simulation of pulverized coal jet flame employing the TDP model. Fuel 2012, 97, 277−287. (38) Hashimoto, N.; Kurose, R.; Hwang, S. M.; Tsuji, H.; Shirai, H. A numerical simulation of pulverized coal combustion employing a tabulated-devolatilization-process model (TDP model). Combust. Flame 2012, 159 (1), 353−366. (39) Jupudi, R. S.; Zamansky, V.; Fletcher, T. H. Prediction of Light Gas Composition in Coal Devolatilization. Energy Fuels 2009, 23, 3063−3067. (40) Chaundhari, K., Development of Advanced Coal Devolatilization and Secondary Pyrolysis Kinetics Models for CFD (and process simulation) Codes. In Department of Chemical Engineering; West Virginia University: Morgantown, WV, 2010. (41) Lu, X. J.; Wang, T. Water-gas shift modeling in coal gasification in an entrained-flow gasifier - Part 2: Gasification application. Fuel 2013, 108, 620−628. (42) Alvarez, L.; Gharebaghi, M.; Jones, J. M.; Pourkashanian, M.; Williams, A.; Riaza, J.; Pevida, C.; Pis, J. J.; Rubiera, F. CFD modeling of oxy-coal combustion: Prediction of burnout, volatile and NO precursors release. Appl. Energy 2013, 104, 653−665. (43) Yan, B. H.; Cheng, Y.; Jin, Y.; Guo, C. Y. Analysis of particle heating and devolatilization during rapid coal pyrolysis in a thermal plasma reactor. Fuel Process. Technol. 2012, 100, 1−10. (44) Jovanovic, R.; Milewska, A.; Swiatkowski, B.; Goanta, A.; Spliethoff, H. Sensitivity analysis of different devolatilisation models on predicting ignition point position during pulverized coal combustion in O-2/N-2 and O-2/CO2 atm. Fuel 2012, 101, 23−37. (45) Yan, J. H.; Zhu, H. M.; Jiang, X. G.; Chi, Y.; Cen, K. F. Analysis of volatile species kinetics during typical medical waste materials pyrolysis using a distributed activation energy model. J. Hazard. Mater. 2009, 162 (2−3), 646−651. (46) Li, X. G.; Wang, J. A.; Wang, F. L.; Ren, C. L. Dust Explosions Modeling. Proc. Int. Conf. Meas. Control Granular Mater., 8th 2009, 358−364. (47) Backreedy, R. I.; Habib, R.; Jones, J. M.; Pourkashanian, M.; Williams, A. An extended coal combustion model. Fuel 1999, 78 (14), 1745−1754. (48) Maffei, T.; Frassoldati, A.; Cuoci, A.; Ranzi, E.; Faravelli, T. Predictive one step kinetic model of coal pyrolysis for CFD applications. Proc. Combust. Inst. 2013, 34 (2), 2401−2410. (49) Oppenheim, A. V.; Schafer, R.W. Discrete-Time Signal Processing; Englewood Cliffs: Prentice-Hall, Inc., 1989. (50) Burcat, A.; Ruscic, B. Third millenium ideal gas and condensed phase thermochemical database for combustion (with update from active thermochemical tables); DOE Scientific and Technical Information: Washington, DC, 2005. (51) Mason, D. M.; Gandhi, K. Formulas for Calculating the Heating Value of Coal and Coal Char: Development, Tests and Uses. In ACS Fall Meeting; American Chemical Society: Washington, DC, 1980; p 235−244. (52) Boie, W. Fuel Technology Calculations. Energietechnic 1953, 3, 309−316. (53) Grummel, E. S.; Davis, I. A. A new Method of Calculating the Calorific Value of a Fuel from its Ultimate Analysis. Fuel 1933, 12, 199−203. (54) Mott, R. A.; Spooner, C. E. The Calorific Value of Carbon in Coal. Fuel 1940, 19, 226−231 , 242−251.. (55) Coal conversion Systems Technical Data Book; Springfield, Va: Institute of Gas Technology, 1978. (56) Selvig, W. A.; Gibson, F.H. Calorific Value of Coal. In Chemistry of Coal Utilization;Lowry, H. H., Ed.; John Wiley: New York, 1945; Vol. 1, p 139. (57) Syamlal, M., Rogers, W.; O’Brien, T. MFIX Documentation Theory Guide, U.S.D.o. In Energy; DOE: Morgantown, 2013. (58) Merric, D. Mathmatical Models of the Thermal Decomposition of Coal: 2. Specific Heats and Heats of Reaction. Fuel 1983, 62 (5), 540−546. 7796


Advanced Chemistry Surrogate Model Development within C3M for

Recommend Documents