ARTICLE pubs.acs.org/est
Numerical Solution of the Polanyi-DR Isotherm in Linear Driving Force Models David M. Lorenzetti* and Michael D. Sohn Indoor Environment Department, Lawrence Berkeley National Laboratory, Berkeley, California, United States
bS Supporting Information ABSTRACT: The PolanyiDubininRadushkevich isotherm has proven useful for modeling the adsorption of volatile organic compounds on microporous materials such as activated carbon. When embedded in a larger dynamic simulation—e.g., of wholebuilding pollutant transport—it is important to solve the sorption relations as quickly as possible. This work compares numerical methods for solving the Polanyi-DR model, in cases where transport to the surface is assumed linear in the bulk-to-surface concentration differences. We focus on developing numerically stable algorithms that converge across a wide range of inputs, including zero concentrations, where the isotherm is undefined. We identify several methods, including a modified Newton-Raphson search, that solve the system 34 times faster than simple bisection. Finally, we present a rule of thumb for identifying when boundary-layer diffusion limits the transport rate enough to justify reducing the model complexity.
’ INTRODUCTION Sorption plays an important role in the transport and fate of volatile and semivolatile organic compounds in buildings. For example, these compounds may sorb to room surfaces and furnishings,1 to activated-carbon filters,2 and even to particles trapped on fibrous media filters in the ventilation system.3 In other applications, sorption filters may be deployed to protect against toxic industrial chemicals and chemical warfare agents.4 Sorption can change airborne pollutant profiles, both reducing peak concentrations, and prolonging exposure at low concentrations. Therefore, sorption models may be needed to assess the health consequences of chemicals with intermittent sources, or with nonlinear dose-response characteristics. For example, suppose the health response to a chemical agent depends on R the toxic load cndt, where c is the airborne concentration.5 Then cutting c{t} in half reduces the toxic load by a factor 0.5n. For an agent such as HCN, with its estimated n = 2.7, this amounts to an 85% reduction. Physical adsorption—the reversible binding of gas-phase sorbate molecules to a sorbent surface—typically is fast compared to the processes that transport the sorbate to and from active sites on the surface.6,7 Therefore sorption models usually assume that the gas and liquid phases remain in equilibrium at the sorbent surface. One class of model, the sorption isotherm, relates the sorbed concentration in a closed isothermal system at equilibrium, qe, to that in air, ce.2,810 The nonlinear PolanyiDubininRadushkevich isotherm has proved useful for modeling the sorption of volatile and semivolatile organic compounds on microporous materials, including r 2011 American Chemical Society
activated carbon.6,9,1113 This isotherm has qe ¼ qmax eðDlnfcmax =ce gÞ
2
ð1Þ
where cmax and qmax are the maximum gas- and liquid-phase concentrations, respectively. Dimensionless parameter D = R0T/ E, where R0 is the universal gas constant (J/mol/K), and T is the absolute temperature (K). The adsorption energy, E (J/mol), characterizes the sorbent-sorbate system. A substantial theory exists for estimating E, either by comparing to a reference vapor,9,12 or as a function of system properties such as the micropore width 11 and the molecular weight of the sorbate.14 Concentrations may have any consistent units, and cmax may be a vapor pressure or a saturation concentration. Sorbed concentrations relate to the liquid density F of the sorbate as qmax = WoF, where Wo is the specific volume of micropores in the sorbent.12 In addition to the isotherm, dynamic sorption models must account for the time rate of contaminant transport to and from the sites where adsorption takes place. External or film diffusion drives sorbate through the air-side boundary layer. Internal or intraparticle diffusion accounts for transport within the sorbent granules, including gas-phase diffusion in the pores, and surface diffusion from site to site on the pore surfaces.6,7,15 The simplest transport model treats these processes as linear in the concentration difference between the sorbent surface and the surrounding bulk regions. A linear model of external diffusion Received: July 11, 2011 Accepted: September 29, 2011 Revised: September 29, 2011 Published: September 29, 2011 10091
dx.doi.org/10.1021/es202359j | Environ. Sci. Technol. 2011, 45, 10091–10095
Environmental Science & Technology
ARTICLE
gives a mass transport rate w ¼ kc ðcb cs Þ
ð2Þ
to the surface, where cb and cs are the bulk and surface concentrations in air. The rate constant kc depends on the area of the sorbent-air interface, and on details of the boundary layer.6,16 Similarly, the linear driving force model gives the rate of internal diffusion as w ¼ kq ðqs qb Þ
ð3Þ
away from the surface, where qs and qb are the surface and bulk sorbed concentrations.6,7,17 For more information on sorption modeling, including alternate formulations to the linear models used here, see.7,15,1820 Applying the equilibrium assumption at the surface gives qs = f{cs}, where f represents the isotherm model of interest. Assuming no mass accumulates at the surface, the rates of external diffusion, eq 2, and internal diffusion, eq 3, must be equal, so w ¼ kc ðcb cs Þ ¼ kq ðf fcs g qb Þ
ð4Þ
In a typical dynamic simulation, eq 4 must be solved for cs. Consider, for example, a lumped-parameter model that divides the physical domain into control volumes, such as room air and sorbent material. A set of ordinary differential equations describes the time rate of change of contaminant mass in each control volume, in terms of the transport rates between them.21 Thus the mass in each volume, and the corresponding bulk concentrations are assumed known at any given time. Then to advance the simulation through time, the solver must find the mass transport rates, w. In eq 4, this requires solving for cs. This work investigates the solution of eq 4 when f is given by the Polanyi-DR model of eq 1. The resulting nonlinear equation has no closed-form solution for cs, so a dynamic simulation must use numerical methods to find the mass transport rate (alternately, it may solve a more complicated differential-algebraic system 22). Particular challenges include the fact that the isotherm is undefined at zero concentration, the high accuracy requirement in the context of iterative differential equation solvers,23 and finiteprecision effects associated with defining and ensuring convergence. Furthermore, because the solution is required so often in the course of solving the dynamic system, it is critical to solve eq 4 as quickly as possible. Therefore, we require fast as well as robust solution methods.
’ MODELING This section describes our modifications to the Polanyi-DR model, and the numerical treatment used to ensure convergence to a solution. Recasting eq 4 in residual form, we define r ¼ kc ðc cb Þ þ kq ðf fcg qb Þ
ð5Þ
Given bulk concentrations cb and qb, we seek c = cs such that r = 0 (i.e., such that the model predicts zero accumulation at the sorbent surface). To apply the residual to the Polanyi-DR isotherm, substitute eq 1 for f. Unfortunately, the isotherm is undefined at cs = 0. Since a general simulation tool cannot adopt the case-specific strategy of avoiding zero initial conditions,6 we modify the isotherm to give qe = 0 at ce = 0. Since sorption theory predicts linear (Henry’s law) behavior at trace concentrations,10 we introduce a linear
regime up to some transition point ct: 8 qt > < 0 e ce e ct : ce ct qe ¼ f fce g ¼ 2 > : ct < ce e cmax : qmax eðDlnfcmax =ce gÞ ð6Þ Given ct, setting qt according to eq 1 ensures continuity at the transition point. In the absence of experimental data to guide the selection of ct, consider setting it to provide a continuous derivative at the transition (a noncontinuous derivative may cause problems solving the larger dynamic system in which the sorption model appears). Since 8 qt > : < 0 e ce e ct dqe df ct ð7Þ ¼ ¼ qe 2 > dce dce : ct < ce e cmax : 2D lnfcmax =ce g ce the required point is ct ¼ cmax e1=2D
2
ð8Þ
The resulting “smooth” transition point has no physical interpretation, and in practice may even prevent the use of the Polanyi-DR isotherm. For example, at 20 °C, a hypothetical system with E = 10 kJ/mol would have a smooth derivative at ct = 2.2 104cmax. For benzene, this corresponds to 59 ppm in air. In human health applications, where 10 ppm benzene represents significant exposure,24 using this cutoff would mean modeling sorption entirely using Henry’s law. Solution in Linear Regime. To solve eq 5 with the modified Polanyi-DR isotherm, first check whether the solution lies in the linear, Henry’s law, regime. To do so, consider the trial solution c+ = ct. If r+ g 0, then cs e ct, and eq 6 gives cs ¼
kc cb þ kq qb qt kc þ kq ct
ð9Þ
If, on the other hand, r+ < 0, then cs > ct, and the solution lies in the nonlinear regime. Solution in Nonlinear Regime. Now suppose the solution lies in the nonlinear regime, where the standard Polanyi-DR isotherm holds. This problem poses a challenge to many rootfinding methods, due to the sharp curvature it can develop. Figure 1 shows normalized residual curves, which range from linear in cs (when internal diffusion limits the overall mass transport rate), to highly nonlinear (when diffusion through the air boundary layer limits the transport rate). We tested a number of numerical methods, modifying standard approaches such as regula falsi, Newton-Raphson, and polynomial interpolation in order to guarantee convergence and improve the solution speed. We also tried inverse polynomial fits, rational function interpolation, and Ridder’s method, plus variations on these, as guided by testing. Our main criterion for considering a method was that it should guarantee a solution in exact arithmetic. See the Supporting Information. All the tested algorithms bracket the root. That is, they maintain bounds cl and cr such that cl < cs < cr. A trial solution c+ must fall in the bracket, and it replaces either cl or cr, depending 10092
dx.doi.org/10.1021/es202359j |Environ. Sci. Technol. 2011, 45, 10091–10095
Environmental Science & Technology
ARTICLE
Figure 1. Families of residual curves for thePolanyi-DR isotherm in the linear driving force model. The residual scaling depends on the rate constants and the bulk concentrations. Here, qmax = cmax.
on the sign of r+. The bracket guarantees convergence, and it also helps to define convergence. The search for cs stops when either jrþ j e jwþ jτrel w
or
cr cl e Δcmin
ð10Þ
The first criterion requires a small error relative to the estimated flow, w+. The second, backup, criterion requires a small bracket, allowing termination in cases where finite-precision effects mean the first criterion cannot succeed—for example, due to a large kq, or when qb ≈ f{cb} so that w ≈ 0. A later section describes setting the relative tolerance τrel w and minimum bracket width Δcmin. To set the initial bracket bounds, we let cl = ct and cr = cmax. This constrains the search to the nonlinear regime. Furthermore, since cs must lie between the bulk concentrations, testing c+ = cb and c+ = f 1{qb} can narrow that initial bracket considerably. The tested methods may adjust c+ in order to be sure of narrowing the bracket, even in the face of finite-precision effects. We observed two problems in particular. First, in finite-precision arithmetic, adding a small change to an existing bound may fail to change the bound. Thus every method except bisection requires c+ g cl + Δcmin. If not, c+ is incremented by Δcmin. A similar test is applied at the right-hand bound. This adjustment also avoids taking unnecessarily small steps. Second, changing one bound by Δcmin may fail the “small bracket” test of eq 10. For example, suppose cr was calculated as cl + Δcmin. Finite-precision effects can still allow cr cl > Δcmin. In this case, the solver would fail to recognize convergence, even though it has effectively bracketed the solution. The end result can be an infinite loop, in which the same trial point is calculated repeatedly. To avoid this, we include a final check that cl < c+ < cr. If not, we take the bracket midpoint. To control roundoff error, we use cþ ¼ cl þ ðcr cl Þ=2
ð11Þ
rather than the more conventional (cl + cr)/2. This final bisection also guarantees a valid trial point for the inverse quadratic method, which, even in exact arithmetic, may jump the bracket.
’ NUMERICAL TESTING To ensure each solution method can solve any reasonable problem, we generated 100 000 sorbent-sorbate systems, by random sampling from wide ranges of possible model parameters. While not all the resulting systems have a real-world analog, our goal was to challenge the solution methods, in order to test their utility for a general-purpose simulation tool. To compare the speed of the methods, we solved each sorbent-sorbate system with a number of bulk concentration inputs, defining 100 equilibrium problems for each system. Thus each solution method was compared on the basis of its robustness and speed at solving eq 4 for one million distinct sorption transport problems. Sorbent-Sorbate System. To generate each system, we sampled the sorbent micropore volume, Wo, from a uniform distribution on 0.30.7cm3/g,6,13,25 and the adsorption energy from 330 kJ/mol.12,13 The sorbate saturation concentration was sampled from a collection of temperature-vapor pressure relationships.24 To tie the sampled values to real properties, we used curves for ammonia (H3N), carbon tetrachloride (CCl4), benzene (C6H6), hexane (C6H14), hydrogen cyanide (HCN), methanol (CH4O), propane (C3H8), and toluene (C7H8). However, the particular sorbates do not matter so much as the range provided for testing (9 103 e cmax e 8 104g/m3). In addition, each sorbate has a representative liquid density (g/cm3), used to find qmax = WoF.12,24 Because the diffusion rate constants depend on the size of the system, we assumed an external surface area of 1 m2. This allowed us to treat w as a mass flux density (g/s/m2), as well as a mass flow rate (g/s), thus making kc analogous to a deposition velocity. Accordingly, each system picks kc between 5 106 and 0.1 m/s.16 To provide a representative range of values for kq, we defined kq ¼ kc ðcmax =qmax Þ10k
ð12Þ
At k = 0, internal and external diffusion have roughly equal control over the mass transport rate. For k , 0, external diffusion is fast relative to internal diffusion, bringing the surface concentration close to that in the bulk air. In Figure 1, k = 2 produces a straight line, because kc is so large compared to kq that the linear term kc(c cb) dominates eq 5. Conversely, k . 0 means internal diffusion is fast relative to external diffusion. This boundary layer diffusion controlled case, in which the external boundary layer limits the mass transport rate, is more challenging numerically, since the dominant term of eq 5 now includes the isotherm relation. For example, k = 2 yields the upper curves in Figure 1. Note that if kq is sufficiently greater than kc, a more tractable model results from assuming the sorbent surface is always at the bulk sorbed concentration, and calculating w = kc(cb f1{qb}) directly.6,15 For testing, we sampled uniformly on 3 e k e 3. Note that, in a real system, these parameter choices are not independent. For example, the adsorption energy and rate constants depend, in part, on the system temperature and sorbate. Since our purpose was to challenge the solution methods, we did not restrict the range of possible systems by introducing such dependencies. Problem Definition. For each sorbent-sorbate system, we solved the linear driving force model for 100 problems. A problem comprises a choice of bulk concentrations. In order to test the full 10093
dx.doi.org/10.1021/es202359j |Environ. Sci. Technol. 2011, 45, 10091–10095
Environmental Science & Technology
Figure 2. Average runtime of each solution method, relative to bisection, for machines A (b) and B(O). The sorting order is based on the average for both machines. The Supporting Information describes the methods denoted by the acronyms.
range of model inputs, we set cb to zero 1% of the time, and to cmax 1% of the time. Another 2% of the problems placed 0 e cb < ct. The remaining tests had ct e cb < cmax, in order to provide data on the speed of the solution methods in the nonlinear regime. Similarly, 1% of the problems had qb = 0, 1% at qmax, and 2% below the linearization cutoff. Another 4% of the time, qb was chosen to put the corresponding equilibrium gas-phase concentration within (1% of the selected cb. Finally, the remaining problems randomly chose qt e qb < qmax. Again, the goal was to define mostly challenging problems, while still testing the full range of possible inputs. Tolerances. The tests reported here used τrel w = 5000ε in eq 10. Machine epsilon, ε, is the smallest positive number that, added to 1, yields a number different from 1. Thus a machine with ε = 2 12 . 1016 has τrel w = 1 10 abs The minimum bracket size Δcmin = crτrel c + τc . Basing the tolerance on the right-hand bound ensures the machine can resolve a change Δcmin to either bound. The relative tolerance abs rel τrel c = 10ε. The absolute tolerance, τc , guards against crτc ≈ 0, and is set three decades smaller than the machine precision. For example, if double-precision numbers effectively have 15 digits, 18 . then τabs c = 10
’ RESULTS AND DISCUSSION Figure 2 shows the time each nonlinear solution algorithm takes to run the tests, normalized by the time required for a bisection search (bis). For each algorithm, we ran the complete suite of tests five times; each point in the figure averages the three runs that remain after discarding the fastest and slowest. We report results for two different machines. Machine A has an x86 processor; machine B is a PowerPC. A different vendor supplied the Fortran compiler on each machine. Most of the algorithms achieve between 25% and 35% of the bisection runtime. The exact ranking depends somewhat on the test conditions. For example, increasing k makes NewtonRaphson perform relatively worse, by increasing the nonlinearity of the model relations, while increasing the number of problems in the linear regime tends to reduce timing differences between
ARTICLE
the methods. However, in general the results of Figure 2 should hold across a wide range of applications. Many of the algorithms approach the solution from one side, repeatedly replacing the left or the right bound. Attempts to break this behavior met with varying success. For example, applying Aitken’s method to regula falsi lets it periodically replace the left bound, reducing the runtime from 63% (rf) to 27% (rfa) of that of bisection. On the other hand, Aitken slows down Newton-Raphson (nr), because replacing the right bound does not affect the calculation of the next iterate. Adjusting c+, in order to avoid placing the trial point within Δcmin of one of the bracket bounds, helps prevent one-sided convergence. Although introduced to overcome finite-precision effects, it also speeds up many of the algorithms, by avoiding the repeated addition of very small increments to one bound. Other implementation details, such as which values to cache, and using the stable bisection of eq 11, also matter. In particular, controlling numerical instability also tends to speed up the solution. For example, finding the larger-magnitude root of the quadratic Ax2 + Bx + C = 0, then dividing into C/A to find the smaller-magnitude root, provides a better estimate of that smaller root than does the standard quadratic formula. This results in both faster runtimes (by reducing the number of iterations needed to converge), and better mass flow estimates (by better estimating the root in the final iteration). The Supporting Information contains more implementation details. The two machines and compilers produce broadly similar results, suggesting that these conclusions should apply to other implementations of the same algorithms. In general, machine B has slightly better times relative to bisection; we speculate that its greater number of processor registers disproportionately benefits the more complicated algorithms. All of the algorithms maintain a bracket on the solution. The tests here initialize the bracket using the bulk concentrations. In a dynamic transport simulation, cs is unlikely to change much from call to call of the solver. Therefore one might further narrow the initial bracket, and speed up the solution, by checking the last concentration calculated for the same surface. We recommend Newton-Raphson as the best generalpurpose alternative. Although some methods run faster, NewtonRaphson is particularly straightforward, with a simple formula and relatively few critical implementation details. However, this recommendation could change for a dynamic simulation that initialized the bracket using the last surface concentration, as suggested above. Since Newton-Raphson always approaches the solution from lower concentrations, it may not be the best choice when the surface concentration is falling with time. Finally, if external (internal) diffusion controls the mass transport rate, then reformulating the model to fix the surface concentration at the bulk sorbed (air) concentration provides the fastest solution of all. As a rule of thumb, this approximation should hold for |k| > 2 in eq 12.
’ ASSOCIATED CONTENT
bS
Supporting Information. Further implementation details. Because most of the solution methods have several possible formulations, shows the particular one tested for Figure 2, and describes the rationale behind each choice. This material is available free of charge via the Internet at http:// pubs.acs.org.
10094
dx.doi.org/10.1021/es202359j |Environ. Sci. Technol. 2011, 45, 10091–10095
Environmental Science & Technology
’ AUTHOR INFORMATION Corresponding Author
*Phone: +1 (510) 486-4562; fax: +1 (510) 486-6658; e-mail:
[email protected].
’ ACKNOWLEDGMENT This research was funded in part by the U.S. Defense Threat Reduction Agency, and was performed under U.S. Department of Energy contract no. DE-AC03-76SF00098. ’ REFERENCES (1) Singer, B. C.; Hodgson, A. T.; Hotchi, T.; Ming, K. Y.; Sextro, R. G.; Wood, E. E.; Brown, N. J. Sorption of organic gases in residential rooms. Atmos. Environ. 2007, 41, 3251–3265. (2) Axley, J. W. Modeling sorption transport in rooms and sorption filtration systems for building air quality analysis. Indoor Air 1993, 3, 298–309. (3) Hyttinen, M.; Pasanen, P.; Kalliokoski, P. Adsorption and desorption of selected VOCs in dust collected on air filters. Atmos. Environ. 2001, 35, 5709–5716. (4) Linders, M. J. G.; Baak, P. J.; van Bokhoven, J. J. G. M. Exploratory investigation of the risk of desorption from activated carbon filters in respiratory protective devices. Ind. Eng. Chem. Res. 2007, 46, 4034–4039. (5) ten Berge, W.; Zwart, A.; Appelman, L. Concentrationtime mortality response relationship of irritant and systemically acting vapours and gases. J. Hazard. Mater. 1986, 13, 301–309. (6) Axley, J. W. Tools for the analysis of gas-phase air-cleaning systems in buildings. ASHRAE Trans. 1994, 110, 1130–1146. (7) Weber, W. J.; Smith, E. H. Simulation and design models for adsorption processes. Environ. Sci. Technol. 1987, 21, 1040–1050. (8) Kleineidam, S.; Sch€uth, C.; Grathwohl, P. Solubility-normalized combined adsorption-partitioning sorption isotherms for organic pollutants. Environ. Sci. Technol. 2002, 36, 4689–4697. (9) Lordgooei, M.; Kim, M. Modeling volatile organic compound sorption in activated carbon. I: Dynamics and single-component equilibrium. J. Environ. Eng. 2004, 130, 212–222. (10) Pikaar, I.; Koelmans, A. A.; van Noort, P. C. Sorption of organic compounds to activated carbons. Evaluation of isotherm models. Chemosphere 2006, 65, 2343–2351. (11) Stoeckli, H. F. Microporous carbons and their characterization: The present state of the art. Carbon 1990, 28, 1–6. (12) Wood, G. Affinity coefficients of the Polanyi/Dubinin adsorption isotherm equations: A review with compilations and correlations. Carbon 2001, 39, 343–356. (13) Wood, G. O. Review and comparisons of D/R models of equilibrium adsorption of binary mixtures of organic vapors on activated carbons. Carbon 2002, 40, 231–239. (14) Duchowicz, P. R.; Casta~neta, H.; Castro, E. A.; Fernandez, F. M.; Vicente, J. L. QSPR prediction of the DubininRadushkevich’s k parameter for the adsorption of organic vapors on BPL carbon. Atmos. Environ. 2006, 40, 2929–2934. (15) Pei, J.; Zhang, J. Modeling of sorbent-based gas filters: Development, verification and experimental validation. Build. Simul. 2010, 3, 75–86. (16) Nazaroff, W. W.; Cass, G. R. Mass-transport aspects of pollutant removal at indoor surfaces. Environ. Int. 1989, 15, 567–584. (17) Delage, F.; Pre, P.; Le Cloirec, P. Mass transfer and warming during adsorption of high concentrations of VOCs on an activated carbon bed: Experimental and theoretical analysis. Environ. Sci. Technol. 2000, 34, 4816–4821. (18) Jørgensen, R. B.; Dokka, T. H.; Bjørseth, O. Introduction of a sink-diffusion model to describe the interaction between volatile organic compounds (VOCs) and material surfaces. Indoor Air 2000, 10, 27–38.
ARTICLE
(19) Lodewyckx, P.; Wood, G. O.; Ryu, S. K. The Wheeler-Jonas equation: A versatile tool for the prediction of carbon bed breakthrough times. Carbon 2004, 42, 1351–1355. (20) Murakami, S.; Kato, S.; Ito, K.; Zhu, Q. Modeling and CFD prediction for diffusion and adsorption within room with various adsorption isotherms. Indoor Air 2003, 13, 20–27. (21) Axley, J. W. Multi-zone dispersal analysis by element assembly. Build. Environ. 1989, 24, 113–130. (22) Brenan, K.; Campbell, S.; Petzold, L. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations; Society for Industrial and Applied Mathematics: Philadelphia, 1996. (23) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes: The Art of Scientific Computing, 3rd ed.; Cambridge University Press: New York, 2007. (24) CRC Handbook of Chemistry and Physics, 72nd ed.; Lide, D. R., Ed.; Chemical Rubber Company: Boca Raton, FL, 19911992. (25) Li, L.; Quinlivan, P. A.; Knappe, D. R. U. Predicting adsorption isotherms for aqueous organic micropollutants from activated carbon and pollutant properties. Environ. Sci. Technol. 2005, 39, 3393–3400.
10095
dx.doi.org/10.1021/es202359j |Environ. Sci. Technol. 2011, 45, 10091–10095
