High fidelity modeling of flare combustion and emission detection

The combustion of hydrocarbons by industrial flares emits greenhouse gases into the atmosphere along with trace amounts of ozone-forming highly-reactive volatile organic compounds and human carcinogens. Continuous and autonomous monitoring of flare emission is challenging due to low species concentrations and varying environmental conditions. We outline the development of a novel combination of modeling and laboratory measurements to ascertain the optimal combination of detectors to monitor the ozone-forming compounds and human carcinogens. The modeling component couples a high-fidelity computational fluid dynamic model of the combustion process with a first-principles radiation transport model for remote detection. This model provides a consistent and complete description of both the combustion and detection process. Predictions made by the model will be tested in an instrumented and controlled flare facility, established at the Air Force Institute of Technology. We outline the coupling of the computational fluid dynamic code (ARCHES), developed by the University of Utah, with the radiation transport model (SAMM®) and provide results from the composite forward model. We then outline the details of the retrieval algorithm used to interpret laboratory measurements.

High-Fidelity Modeling of Flare Combustion and Emission Detection R. Panfili*a, H. Dothea, P. Vujkovic-Cvijina, X. Tana, R. Kennetta, R. Taylora, L. Bernsteina, P.J. Smithb, J.N. Thornockb, K. Grossc, J. Seeboldd aSpectral Sciences, Inc., 4 Fourth Avenue, Burlington, MA 01803; bUniversity of Utah, Salt Lake City, UT; cAir Force Institute of Technology, Wright-Patterson AFB, OH 45433; dChevron Research and Technology Company (Ret.), Richmond, CA 94802, ABSTRACT The combustion of hydrocarbons by industrial flares emits greenhouse gases into the atmosphere along with trace amounts of ozone-forming highly-reactive volatile organic compounds and human carcinogens. Continuous and autonomous monitoring of flare emission is challenging due to low species concentrations and varying environmental conditions. We outline the development of a novel combination of modeling and laboratory measurements to ascertain the optimal combination of detectors to monitor the ozone-forming compounds and human carcinogens. The modeling component couples a high-fidelity computational fluid dynamic model of the combustion process with a first-principles radiation transport model for remote detection. This model provides a consistent and complete description of both the combustion and detection process. Predictions made by the model will be tested in an instrumented and controlled flare facility, established at the Air Force Institute of Technology. We outline the coupling of the computational fluid dynamic code (ARCHES), developed by the University of Utah, with the radiation transport model (SAMM®) and provide results from the composite forward model. We then outline the details of the retrieval algorithm used to interpret laboratory measurements. Keywords: Industrial Pollution, Remote Monitoring 1. INTRODUCTION Modeling and simulation of complex physical phenomena can be a valuable component of any design and development process. The construction of a model and the use of its simulations can substantially reduce the costs associated with the development of new technology. Design concepts can be evaluated in a virtual environment at a fraction of the cost of the assembly and field testing of an equivalent system. Automotive design, drug discovery and petroleum exploration are just a few examples of applications within industries making use of modeling and simulation in their design process. We seek here to model the remote detection of emissions from the combustion of hydrocarbons by an industrial flare. The model will allow us to determine the fitness for purpose of existing monitoring technology as well as provide us with a means to explore emerging sensor technology. It will also enable us to explore the detection process for a large range of environmental conditions and flare designs. As a result, the model will provide an economical complement to controlled field measurements. Those field measurements, including the most recent Marathon Oil Field Tests [1] and Texas Commission on Environmental Quality tests in Tulsa [2], provide important, validated data sets. Validated data sets obtained from field measurements will always be an important validation component for numerical models but they are unlikely to obviate the utility of mathematical modeling. The field measurements will almost always be constrained by the choice of hardware to test and, critically, the environmental conditions present during the tests. Consequently, the ability to extrapolate from those data sets to more general scenarios can always be called into question. Numerical modeling does not suffer from such limitations and enable us to explore flare tip designs and environmental conditions which are difficult to explore experimentally. While the concept of modeling the detection process is applicable for a host of technologies, we focus specifically on the remote detection of emission products in this effort. The remote detection process involves first capturing the photons emitted by the combustion products and then resolving the molecular distribution that most likely leads to those emissions. This model utilizes several well established codes to provide the relevant physics. The ARCHES code [3] is used to describe the combustion process. The SAG code [4], augmented with available weather data, is used to describe *[e-mail]; phone [phone]; fax [phone]; spectral.com SAMM Computer Software is a registered trademark owned by the United States Government as represented by the Secretary of the Air Force the background environment. The SAMM code [5] uses the molecular concentrations from ARCHES and the atmospheric profiles from SAG and performs the radiation transport to the detector. The composite output o f these three codes is then used as predictors for an inversion or retrieval algorithm. This inversion algorithm [6], conceptually similar to the method to be used for data analysis on the NASA Orbiting Carbon Observatory mission [7], is discussed in greater detail. A validated modeling and prediction effort consists o f four interconnected elements. The first is the construction o f the model itself, some elements o f which are discussed here. Next is the process o f using the model to analyze measurements whose results can be independently validated. Remote sensing measurements acquired during the recent combustion flare tests in Tulsa from two separate perspectives serve this purpose. In the third step, the model is employed to make predictions for physical scenarios outside of the range in which measurements have previously been taken. The fourth and final step is to perform experiments consistent with conditions outlined in the previous step. This step is repeated until a sufficient parameter space is covered to give confidence in the predictive powers of the model. We are in the process o f establishing a controlled and fully instrumented combustion flare laboratory for this purpose. This paper is separated into this introduction and four additional sections. In section 2 we provide a brief review o f the science o f the combustion and radiation transport processes. In section 3 we give details o f the inversion algorithm used to predict combustion products from radiance measurements. In section 4 we illustrate the process with a selection of example computations. Finally, in Section 5 we give conclusions o f this work to date and outline future directions for this project. 2. BACKGROUND The key data set we operate upon is the amount o f emitted light which is directed along a line-of-sight to a detector. When this information is obtained from a set o f experiments, it can be operated upon directly. When this information is obtained from a simulation, the results can be recorded at arbitrarily high spatial, temporal, and spectral resolution. The results can later be made to simulate those o f a specific detector by convolution with the instrument response function and instrument noise characteristics followed by a degrading o f the output to conform to a detector's spatial, spectral, and temporal characteristics. It is from this data collect (or simulation output) that an inversion algorithm must then provide a good estimate of the temperature distribution and profile o f molecular concentrations which led to the observed radiance values. Explicitly solving the radiation transport equation also requires us to define the atmosphere through which photons originating from the combustion process propagate. It is safe to assume that some very basic information regarding the measurement scenario is always available during the detection process. This includes details o f the measurement location (latitude and longitude), date, and time. These inputs are used to define a baseline model atmosphere using SAG [4] and the solar irradiance using a Kurucz model [8]. If no additional information is provided, SAG can provide a reasonable estimate of atmospheric temperatures and densities for a given location and time o f day. If specific weatherrelated data, such as local temperature, pressure and humidity at the measurement location, is available it can be incorporated into the atmospheric profile in a consistent manner. The ambient atmosphere generated by SAG provides the radiation transport model with a vertically stratified atmosphere, schematically represented in Figure 1. An observer can be arbitrarily placed within that vertically stratified atmosphere, and incident radiation can be observed for an arbitrary viewing angle. Output from the combustion model is embedded within this global environment and represented as a cube placed between two atmospheric layers. Atmospheric properties (temperature and molecular densities) for the local region are supplied by the computational fluid dynamics (CFD) simulation. In the event that a radiator is naturally occurring in the atmosphere but has not been recorded by the CFD simulation background values are used. Figure 1. Illustration of the vertically stratified atmosphere utilized by the radiation transport algorithm. The cube inserted between two layers represents the local flare combustion region whose properties are supplied by the CFD simulation. An example o f the type o f information used by the inversion algorithm is shown in Figure 2. In this example, an image o f the radiance due to the ambient atmospheric background is spatially resolved in the left-hand image. The registered band pass radiance o f each pixel can further be spectrally resolved, as is shown in the right-hand image. It is the spectrally resolved radiance o f a single pixel which is used by the inversion algorithm. 1.4x10 1.2x10 1.0x10 8.0x 10 2.0x 10 70 0 Frequency [cm-1] Figure 2. Broadband integrated radiance of the atmospheric background (left) and spectral radiance of an individual pixel within the scene (right). The requirements o f our measurements, namely the need to measure the products o f the process, dictate that the sensor must observe the plume far enough away from a flame that combustion chemistry has completed. Properly modeling the propagation o f light from the atmosphere and the flare requires a radiation transport capability to propagate molecular emission from the combustion products to the sensor. The SAMM model is capable o f solving the radiation transport equation for arbitrary viewing geometries. The direct solution o f the physical and chemical processes involved in combustion dynamics and radiation transport (the forward problem) is unique in the sense that a specific set of initial conditions yield a specific result. A more detailed description of how SAMM is being utilized to model radiance emitted from flare combustion products is available [9]. Those results are then used as inputs in the iterative processes of retrieving temperature and density profiles from the observed radiance. 3. INVERSION MODEL The basic equations necessary to construct the inversion model are the path radiance (I) and first derivatives o f the path radiance. We first express the path radiance in a general manner and then proceed to discretize that expression to conform to the equations directly being solved in the inversion model. 3.1 Radiance Path radiance arising from absorption and emission in the ambient atmosphere can be defined by the radiation transport equation. The radiance from a local radiance source S a , observed by a sensor at a distance R, including LOS path radiance, is given by, R I a = S «T« + J T« (r) £ J m,t (r)Kaj. (r)p ( r V r (1) i A summation over all o f the transition lines from all species that contribute to the absorption at a is implied in Eq. (1), where Kai (r) is the frequency dependent extinction coefficient due to an absorption line i , and p" is the corresponding population o f the lower state o f the transition, at the distance r. The transmittance at r is defined by Ta(r )) = - *e-- aa ) (2) where the optical depth is given by (r) = J 2 * a .i (r' " ' (r') dr ' (3) i 0 We can immediately identify the first term on the right-hand side o f Equation (1) as the transmission at the end o f the LOS. This can be expressed as: T a - Ta(R) (4) Next, we consider the differentiation o f key quantities in the expression for radiance. Differentiating the transmission gives an expression the change in transmission as a function o f the change in optical depth. dTa a = - e -tadTa a = T.,dr a a (5) ' Likewise, we can further expand the differentiation o f optical depth by using the definition provided in Equation 3. Differentiating that expression yields the following relation. dTa = Z dTa,i. = Z Ka,i ( r ) p (r^ i (6) i We can now take these expressions and substitute them into the second term on the right-hand side o f Equation 1. This gives us the following path radiance expression. R J ( J a (r))T«(r)X Ka (r )p "( r ) r = J (J a (r 'd*a = - J (J a (r)dTa (7) 0 where an averaged source function has been defined as, Z J a (r )KoAr p ( r ')dr J a (r ) = ( ) "( d Z K< ‘*(r )p i(r ) r V (8) i Discretizing the LOS into path segments k, and denoting the segment averaged source function as J ) , the second term on the right-hand side o f Equation 1 becomes, Z J ')(e-'(k-"- e-'")) k So Eq. (1) can now be rewritten as, (9) I® = S J a + ^ ( j ka ) (l - e-A*k )e-Tk-1 (10) k Where A *k = *k - *k-1 (11) Now we recall that the source function, for the transition L o f species i, is defined in general by, cl p 'l(r V g L J l (r ) = * (12) p'L(r )/g'L - p ( r )lgL where p are the population densities, g are the degeneracy factors, the primes and double primes denote the upper and lower state o f the transition L, and the species index has been suppressed for clarity purposes. The expression in Equation (12) is general, and is the one appropriate to describe both local thermodynamic equilibrium (LTE) and nonLTE conditions. Under LTE conditions, the source function reduces to the black body function for all species and the expression in Equation (12) simply reduces to the black body function at temperature T(r), J l (r ) = B®(T ( ( ) = ^ n 1 1 e ' v!-1 ( 13) We can now substitute the black body function into the source term for the discretized path radiance expression. I® = + ^ j l ) (1 - e -AZk ) k k = S me-*L + £ B m(Tk )(l - e -A*k) -*k-1 (14) 3.2 Retrieval In the retrieval problem, one seeks to retrieve species densities and temperatures along the radiance path, by minimizing the expression, F (P) = Z ( 7 * - 1 ( L , P))2 J (15) where I (®, p) = I® is given by Eq. (5) and the set o f parameters P are the segment densities and temperatures p t (m) and Tm. The minimization can be performed by algorithms such as the Gauss-Newton method [10], or the LevenbergMarquardt [11,12] method. Both methods enjoy several advantages over alternative approaches. For example, both methods require only first derivatives in their analysis, avoiding the need to determine higher-order derivatives. Additionally, both methods do not rely on a local linearity of the solution, a basic assumption made by methods based on series expansion. The Levenberg-Marquardt further possesses advantages over Gauss-Newton for the current problem at hand. Specifically, it is more likely to converge to an accurate solution for a poor initial guess than the Gauss-Newton method. This comes at a cost o f slightly longer convergence time. Given the context o f this problem, specifically that combustion emission products may span a large parameter range, trading a little extra computation time for greater robustness represents a sensible design decision. This, and other retrieval algorithms, generally requires derivatives o f the equations being optimized. In both algorithms, the derivatives, dI® / dpt (m) and dI® / dTm are required. The evaluation o f the derivatives o f the LTE radiance I® , involves application of the chain rule, and ultimately require the partial derivatives with respect to density and temperature o f the optical depths and o f the temperature derivatives o f the black body function , drk (L)/ dpt (m ), dxk (® )/ dTm and dB® / dTmwhere k and m are layer indices. Setting the lower state densities, p "(m)= pi(m)e~E The expression for optical depth becomes (16) Tk( k '_ Z Z (17) \ i p (m)e_E'(i)kTmDm m_l i,i Where E" (i) is the energy o f the lower state o f the line transition i , Dmis the length of the path segment m . Recall that path transmission is directly related to optical depth Ta = e~Ta . Therefore, derivatives o f the transmission with respect to an arbitrary variable (A) can be reduced to the derivative o f the optical depth with respect to that same variable. dTk - k = -d e _T Ta _ _ dlkk e dA dA dA _ _T dikk dA (18) The variable N is the total number o f path segments and the derivatives of optical depth for each layer along the line-ofsight with respect to lower state density can be expressed as. dTk( k _ dPi(m) D v (m) ( ) _E'(iVkTm m '*■ai m i,i m <k 0 m >k (19) We can also analytically compute the derivative o f the optical depth with respect to the temperature at each layer d*k ( k ) _ dTm D, i,i _E-(i)kTm E ''(i) J.m) 2 ai ( ) - dTm kT m <k (20) m >k 0 The final term which needs to be differentiated is the source function. Under LTE conditions, that term reduces to the black body, or Planck, function. The Planck function has no explicit molecular density dependence, leaving only differentiation with respect to temperature. As with the optical depth and transmission terms, the temperature derivative o f the black body function can be expressed analytically. ha haikTk -1 dTm 0 1 m _k kBTk Tk (21) m^k 4. DEMONSTRATION We demonstrate the retrieval algorithm by applying it to determine the densities and temperatures from path radiance calculations given by the SAMM code. The atmospheric profile is provided by SAG. The scenario chosen is a threesegment path from observer to source. The radiator in this demonstration is CO2, and the band pass selected is the 2.06^m region. This spectral region is one o f the two bands used by the Orbiting Carbon Observatory [6] to determine carbon dioxide path quantities in the atmosphere. The "observed" radiances are generated by a forward calculation using a set o f densities and temperatures for three 1-km segments. Only the major isotope o f carbon dioxide is included. The values used to generate the observed radiances are densities o f 7.88x1015, 7.15 x1015, and 6.46 x1015 molecules/cm3, and temperatures o f 285 K, 272 K, and 268 K for the first, second and third segment, respectively. Two preliminary tests of the retrieval algorithm were performed. The first case involves keeping the temperatures constant (at the initial values used to generate the observed radiances) and retrieving the densities. The second case involves keeping the densities fixed and retrieving the temperatures. Both tests require starting from some reasonable initial guess o f the quantities to be retrieved. The results are shown in Figures 3 and 4 for the two cases. The convergence o f the inversion algorithm for density retrieval as a function o f the number o f iterations performed is illustrated in Figure 3. Two separate initial guesses were used and the Levenberg-Marquardt method converged to the correct solution in both cases. Notice as well that the initial guess did not correctly predict the basic trend o f densities along the line-of-sight. This did not impede the algorithm from iteratively reaching the correct solution. Segment Segment Figure 3. Illustration of convergence for density retrievals, starting from two different guesses. The temperatures are fixed at initial values described in the text. The convergence of the inversion algorithm for temperature retrieval as a function o f the number o f iterations performed is illustrated in Figure 4. Once again, the Levenberg-Marquardt method converged to the correct solution despite being given an initial poor guess and incorrect basic trend o f the temperature profile. Ite r Ite r Ite r Ite r Ite r Ite r Ite r 1 2 3 4 5 6 7 • Ite r Ite r Ite r Ite r Ite r 8 9 10 11 12 Ite r 13 Ite r 14 Segment Figure 4. Illustration of convergence for temperature retrievals. The densities are fixed at initial values described in the text. The figures illustrate the convergence o f the results toward the initial values listed above. It can be seen that the density retrieval converges much faster. We note again that the derivative o f radiance with respect to the source term does not explicitly depend on density. Derivatives o f radiance with respect to temperature are thus expected to be the more non­ linear o f the two and this may be related to why convergence of the temperature profile requires more iterations than the convergence o f the density profile. An illustration o f the "observed" and initial guess spectra from the temperature retrieval test is provided in Figure 5. While the initial temperature guess was no greater than 50 K from the correct result, the composite radiance discrepancy was as large as a factor of two. 3.0x10-12 2 .5 x 1 0 -12 1 2 .0 x 1 0 -12 ' E „o 0 1.5x10-12 1 0 ° § 1.0x10 i2 5 .0 x 1 0 -13 0.0 4800 4850 F re a u e n c v (cm-1) Figure 5. Comparison of observed spectrum and initial guess of the temperature retrieval test. The spectral residuals (observed spectrum minus modeled spectrum) as functions of iterations are shown in Figure 6. Notice that the scale o f the residuals has fallen by two orders o f magnitude between the first and seventh iteration and continues to fall further as we increase the number o f iterations performed. Freqency (cm-1) Freqency (cm"1) Figure 6. Spectral residuals as functions of iterations. Note that the y axis scale in (b) is 100 times smaller than that in (a). 5. CONCLUSIONS AND FUTURE DIRECTION The use o f high-fidelity models can be useful for exploring the feasibility o f measuring combustion emission products. We have illustrated a method to consistently model the combustion and measurement process within a single simulation. The approach utilizes the output o f a model o f flare combustion process and combines it with an atmospheric background and radiation transport capability. The direct computation o f the radiance o f the combined combustion product and background environment has been previously demonstrated [9]. Those computations are used as iterative inputs to an inversion algorithm whose purpose is to determine the temperature and molecular concentration profile along a line-of-sight from observed radiance. The effectiveness o f this inversion, or retrieval, algorithm has been demonstrated here. The inversion algorithm completes the major elements necessary to predict hydrocarbon concentrations at a specific point in space temporal interval. This capability will next be combined with a temporal and volumetric integration of the flow out o f the emission plume region to give the averaged combustion efficiency and trace species production for the entire combustion process. The results o f this modeling capability w ill be tested against data acquired as part o f the Texas Commission on Environmental Quality flare study. Additional tests o f the modeling capability w ill be performed at a combustion flare laboratory being established for this specific purpose at the Air Force Institute o f Technology. ACKNOWLEDGEMENTS Funding for this research comes from the Department o f Energy through Contract No. DE-SC003373. support has been provided to this project by Dr. Jeffrey A. Mercier o f Sandia National Laboratory. Substantial REFERENCES [1] Cade, R., Evans, S., et al., "Performance Test o f a Steam-Assisted Elevated Flare with Passive FTIR ," Marathon Petroleum Company, LLC Final Report, 2010. [2] Allen, D.T., and Torres, V.M., "TCEQ 2010 Flare Study Final Report ," PGA No. 582-8-862-45-FY09-04 (2011). [3] Shroll, R. M., S. Adler-Golden, J. W. D uff and J. H. Brown, "Users' Manual for SAG-2, SHARC/SAMM Atmosphere Generator," AFRL-TR-03-1530, 2003. [4] Spinti, J.P., Thornock, J.N., Eddings, E.G., Smith, P.J., and Sarofim, A.F., "Heat Transfer to Objects in Pool Fires," in Transport Phenomena in Fires, WIT Press, Southampton, U.K. (2008). [5] H. Dothe, H., Duff, J.W., Panfili, R., Gruninger, J.H., Kennett, R.G., Acharya, P.K., Berk, A., and Bernstein, L.S., "Users' Manual for SAMM® Version 3," SSI-TR-599, 2010. [6] Connor, B.J., Boesch, H., Toon, G., Sen, B., Miller, C., and Crisp, D., "Orbiting Carbon Observatory: Inverse method and prospective error analysis, Journal o f Geophysical Research 113, D05305 (2008); doi:10.1029/2006JD008336 [7] Crisp, D., Atlas, R.M., Breon, F.-M., Brown, L.R., Burrows, J.P., Ciais, P., Connor, B.J., Doney, S.C., Fung, I.Y., Jacob, D.J., Miller, C.E., O'Brien, D., Pawson, S., Randerson, J.T., Rayner, P., Salawitch, R.J., Sander, S.P., Sen, B., Stephens, G.L., Tans, P.P., Toon, G.C., Wennberg, P.O., Wofsy, S.C., Yung, Y.L., Kuang, Z., Chudasama, B., Sprague, G., Weiss, B., Pollock, R., Kenyon, D. and Schroll, S., "The Orbiting Carbon Observatory (OCO) mission," Advances in Space Research 34, 700-709 (2004); doi:10.1016/j.asr.2003.08.062 [8] Kurucz, R. (1993), available from http://kurucz.harvard.edu. [9] Panfili, R., Vujkovic-Cvijin, P., Tan, X., Kennett, R, Taylor, R., Dothe, H., Bernstein, L., Smith, P.J., Thornock, J.N., Gross, K., and Seebold, J., "Remote Detection o f Volatile Organic Compound Emissions from Combustion Flares ," 33rd Annual Review o f Atmospheric Transmission Models Meeting (2011). [10] Hartley, H.O., " The Modified Gauss-Newton Method for the Fitting o f Non-Linear Regression Functions by Least Square," Technometrics 3, 269-280 (1961). [11] Marquardt, D.W., "An algorithm fo r least-squares estimation o f nonlinear parameters ," Journal o f the society for Industrial and Applied Mathematics 11, 431-441 (1963). [12] Levenberg, K., "A method for the solution o f certain non-linear problems in least squares," Quart. Appl. Math. 2, 164-168 (1944).