Low-Order Modeling of Internal Heat Transfer in Biomass Particle

May 11, 2016 - Comparisons of transient temperature trends predicted by our one-dimensional method with three-dimensional simulations of woody biomass...
3 downloads 23 Views 4MB Size
Download PDF
Article pubs.acs.org/EF

Low-Order Modeling of Internal Heat Transfer in Biomass Particle Pyrolysis Gavin M. Wiggins,*,† Peter N. Ciesielski,‡ and C. Stuart Daw† †

Oak Ridge National Laboratory, 2360 Cherahala Boulevard, Knoxville, Tennessee 37932, United States National Renewable Energy Laboratory, 15013 Denver West Parkway, Golden, Colorado 80401, United States

‡

S Supporting Information *

ABSTRACT: We present a computationally efficient, one-dimensional simulation methodology for biomass particle heating under conditions typical of fast pyrolysis. Our methodology is based on identifying the rate limiting geometric and structural factors for conductive heat transport in biomass particle models with realistic morphology to develop low-order approximations that behave appropriately. Comparisons of transient temperature trends predicted by our one-dimensional method with threedimensional simulations of woody biomass particles reveal good agreement, if the appropriate equivalent spherical diameter and bulk thermal properties are used. We conclude that, for particle sizes and heating regimes typical of fast pyrolysis, it is possible to simulate biomass particle heating with reasonable accuracy and minimal computational overhead, even when variable size, aspherical shape, anisotropic conductivity, and complex, species-specific internal pore geometry are incorporated.

1. INTRODUCTION Pyrolysis is a promising thermochemical pathway under investigation for biomass conversion into liquid fuel forerunners, which are frequently referred to in the literature1 as bio-oils or tars. Biomass pyrolysis is typically carried out in bubbling or circulating fluidized-bed reactors, but other types of reactors are also used.2,3 An extensive body of research for many reactor designs has identified the importance of particle heating rate on the yields of liquids, light gases, and char.4−7 In general, the maximum yields of bio-oil are achieved with “fast” pyrolysis conditions, which involve very high particle-scale heating rates of up to 1000−10 000 °C/s, and peak temperatures near 500 °C.1,8−10 Consequently for bio-oil production, there is considerable interest in developing accurate models that can predict fast pyrolysis heating rates for industrially relevant biomass feedstocks as functions of the material properties and process conditions. Previous experimental and analytical studies have demonstrated that the most significant intrinsic (intraparticle) factors controlling heating rate are particle size and shape, composition (including moisture), density, and microstructure.7,11−28 It is not surprising that particle size should be important since it directly affects the distance through which external heat must be conducted from the surface to the particle center. However, a feature of biomass that is often not fully appreciated (and often not readily evident from reported biomass particle size data) is that biomass typically fragments into particles that are highly asymmetric due to its high aspect-ratio and directional porosity; therefore, thickness can vary significantly with orientation. This is illustrated in Figure 1, which includes micrographs of typical biomass feedstocks after mechanical grinding. Studies have revealed that overall particle shapes can vary substantially between feedstock species as well as with the method of comminution.29 The dominant particle geometry displayed by many types of biomass is roughly cylindrical, although a range of shapes and aspect ratios are usually present © XXXX American Chemical Society

Figure 1. Micrographs of common pyrolysis feedstocks highlight the nonspherical geometry of the particles as well as variations in shape and size. Particles were ground with a Wiley knife mill to pass a 1 mm screen mesh.

and can vary depending on the level of screening. Such variations in shape impart large differences between surfacearea to volume ratios of the particles, which in turn affects the heating rates experienced by the particles due to differences in surface area available for heat transfer. Some recent modeling investigations of pyrolysis particle heating have assumed particle shapes that conform to ideal geometries such as cubes, cylinders, spheres, or tetrahedrons.12,17,19,22,23,30,31 This facilitates the development of onedimensional (1-D) approximations for intraparticle heat transfer but requires designation of an appropriate length scale to represent the dominant transport direction. In some situations, these heat transfer models have been successfully combined with simple kinetics to simulate laboratory experiReceived: March 7, 2016 Revised: May 10, 2016

A

DOI: 10.1021/acs.energyfuels.6b00554 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels ments.8 Lu et al.16 proposed a procedure for using 1-D simulation that, at least theoretically, can be used with any overall particle shape. Other investigators have utilized higherdimensional particle models to account for the effects of multidimensional transport.6,8,24,27,32,33 In the case of 2-D particle models that include reaction kinetics, it has been reported that the predicted particle conversion rates and char yields tend to be higher and lower, respectively, compared to 1D models and experiments. These discrepancies become larger as particle size increases.6,32 Particle heating rates are also strongly affected by intraparticle structure and composition, since these characteristics dictate the thermal properties of the particle and can vary substantially among different types of biomass. A major source of this variation is interspecies differences in porosity due to the cellular structure inherent to plant tissue, which was utilized by the living organism for transport of water and nutrients. For example, straw, perennial grasses, and hardwoods display vastly different tissue structure, which manifests as large differences in the materials’ properties such as density (320−720 kg/m3 drybasis for softwoods and hardwoods34,35 ) and thermal conductivity. Ultimately the bio-oil produced from these feedstocks under the same pyrolysis conditions varies significantly in terms of yield and composition.21 Since the thermal properties of wood are largely affected by its porosity, interspecies variations in wood microstructure impart differences to its thermal properties. In addition, the microstructure of wood is highly directional, causing the thermal properties to vary directionally in relation to the grain.36 Experiments have shown that wood thermal conductivity along the grain is an average of 1.8 times greater than conductivity across the grain.34 Such directional effects cannot be directly accounted for in 1-D particle models. Thus, when microstructural differences are significant, it may be necessary to utilize more computationally intensive 2-D or 3-D models for resolving all the details of particle heating. To date, there are only a few examples of models that include particle anisotropy.8,18,32,35 A much larger number of proposed particle heating models do not consider this issue.5,12,14,15,17,19,20,22,23,25,26,30,37−47 An additional complication in analyzing intraparticle heating is the variation of thermal properties over time. Such variations impart nonlinearities into the simulation governing equations which arise from several phenomena: (1) the temperature dependence of the heat capacity and thermal conductivity of the material, (2) changes in the composition of the material as it undergoes pyrolytic transformation, and (3) changes in the physical structure (including overall particle shape, size, and microstructure) as pyrolysis proceeds. Because of the complexity of these effects and lack of reliable experimental measurements, many particle heating models developed so far assume constant internal structure and properties.12,23,25,26 However, some recent studies have attempted to include the impact of these effects.18,27,35 Our objective in this study is to focus on the conductive heat transfer component of biomass particle heating during fast pyrolysis. Our methodology is based on comparing the predictions of a low-order 1-D heat conduction model to the behavior of biomass particle models with more complex and anatomically accurate structure in 3-D simulations. Consequently, the low-order model is coupled to a pyrolysis kinetics scheme to investigate the effect of reaction heat on particle heat-up. We expect this approach will lead to useful rules-of-

thumb for selecting an appropriate equivalent characteristic dimension and average thermal properties to use for a 1-D model. We consider this a first step in constructing a computationally efficient, but reasonably accurate model of the dominant particle-scale processes that limit the initial yields and composition of the nascent bio-oil exiting biomass particles in a pyrolysis reactor. Our ultimate goal is to combine the completed particle-scale model with other models that account for reactor-scale effects and thus simulate process-scale performance for candidate pyrolysis concepts with realistic biomass feedstocks. Such models are needed to reduce the ambiguity in optimizing bio-oil production processes that can economically handle a variety of biomass feedstocks and blends. We are particularly interested in identifying lower-order approximations that provide reasonably accurate results without expensive computational resources.

2. TECHNICAL APPROACH Our strategy for establishing a low-order methodology for estimating intraparticle biomass heating rate was based on the hypothesis that, by choosing the correct characteristic dimension and average properties, it should be possible to approximate three-dimensional (3-D) heat conduction in realistic biomass particles with a one-dimensional (1-D) model. We pursued this strategy by investigating various assumptions for configuring 1-D approximations of the intraparticle heat transfer to see how well they replicate the 3-D behavior. Finally, we selected which of the 1-D approximations yielded the closest match to the 3-D simulation results. We also combined the resulting 1-D heat transfer model with simplified pyrolysis kinetics to see how well it compared with published experimental measurements of particle temperature transients at fast pyrolysis conditions. 2.1. Assumptions About Biomass Thermal Properties. Given the primary goal of understanding and modeling heat conduction in realistic biomass particles, it is important to establish appropriate parameter ranges for particle geometry and thermal properties to use in computational simulations. Unfortunately, the assumptions made for these parameters in previously reported studies4,5,9,11,23,26,27,37,38,42−45,48−51 have varied considerably and have not led to results that can be easily compared (see Table S1 in the Supporting Information). To reduce ambiguity in this study, we chose woody biomass with typical characteristics taken from recent national lab studies18 and the Wood Handbook.34 Even though we assumed properties typical for wood in this study, we expect that the same general approach is applicable to other types of biomass. Our recent studies of particle-scale heat transfer within detailed 3-D models of biomass particles with explicit microstructure have demonstrated that effects of biomass porosity can be reasonably approximated in heat transfer simulations by employing bulk, isotropic material properties, as long as external particle shape is accurate.18 According to the Wood Handbook, wood properties such as bulk density (ρ) and bulk thermal conductivity (k) vary between species, while bulk heat capacity (Cp) is practically independent of species.34 Expected ranges of specific gravity (SG) and thermal conductivity of various hardwood and softwood species are summarized in Table 1. We note here that the above thermal conductivity values are bulk values, therefore we expect the values to reflect averaging among the different internal directions relative to the grain. However, there is one significant constraint on the thermal conductivity of wood as illustrated in Figure 2, where it can be seen that there is a linear relationship between bulk thermal conductivity and bulk specific gravity. This relationship implies that the bulk thermal properties of woody biomass are largely controlled by the cell walls, where all of the particle mass is concentrated. Interestingly, the intrinsic density of woody biomass cell walls is relatively constant (approximately 1.5 g/ cm3) across wood types, even though bulk properties vary among species.1 The linear relationship in Figure 2 can be approximated as B

DOI: 10.1021/acs.energyfuels.6b00554 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels Table 1. Range of Specific Gravity (SG) and Thermal Conductivity (k) for Various Hardwood and Softwood Speciesa property

hardwoods

softwoods

SG (−) kdry (W/m·K) k12%mc (W/m·K)

0.35−0.78 (0.57) 0.09−0.17 (0.13) 0.10−0.21 (0.16)

0.31−0.62 (0.44) 0.08−0.14 (0.11) 0.09−0.17 (0.13)

a

Calculated from data in Wood Handbook. 34 Dry thermal conductivity at 0% moisture is shown as kdry, and wet thermal conductivity at 12% moisture content was written as k12%mc, while mean property values are given in parentheses.

Figure 3. Three-dimensional solid geometries of particle models with shapes representative of a milled pine feedstock. (inset) Zoomed view showing geometry of small particle models. between a pair of parallel lines that are tangent to the largest projected outline of the particle. A logarithmic scale was used to extrapolate the generic particle dimensions to generate a range of particle sizes as Feret diameters (DF) spanning from 200 μm to 20 mm. These models were employed in 3-D simulations of heat conduction (as presented in Figure 3). Other effects due to feed moisture variations and particle expansion/shrinkage are of practical importance when comparing different process designs, but we did not consider these here to reduce the computational complexity. Therefore, the simulation results discussed in this article assume dry wood (no moisture) and no variation of particle size during heating. 2.3. Computational Simulation of 3-D Heat Conduction. As described above, realistic 3-D shapes were considered based on experimental imaging measurements of milled pine.18 Transient intraparticle heat conduction was simulated for these particles assuming the Fickian heat conduction with the thermal properties discussed previously. Three-dimensional heat conduction can be represented as

Figure 2. Comparison of thermal conductivity (dry and 12% mc) to specific gravity for various hardwood and softwood species. Each symbol represents data for an individual wood species as provided in Table 4−7 of the Wood Handbook.34 kcell wall(dry) = 0.1966 SG + 0.01776 kcellwall(12%mc) = 0.2394 SG + 0.02042

(1)

Extrapolating the above relationships to completely nonporous particles implies that the thermal conductivity of cell walls should be approximately 0.32 W/m·K for dry wood and 0.39 W/m·K for wood with 12% moisture. Equation 2 summarizes the assumed relationship between the bulk heat capacity of dry wood Cpo (kJ/kg·K) and temperature T (K).34

C po = 0.1031 + 0.003867 T

ρC P

∂T = ∇·(k∇T ) ∂t

(3)

where ρ is the bulk density, Cp is the bulk heat capacity at constant pressure, k is the bulk thermal conductivity, and T is temperature. The parameters ρ and k were held constant, while Cp was assumed to vary with temperature according to eq 2. No shrinkage or expansion of particles was included. A convection boundary condition was applied to the surfaces of the particles such that

(2)

According to the Wood Handbook, increases in temperature result in corresponding increases in bulk thermal conductivity at a rate of 2−3% per 10 °C increase.34 This implies that, to a first approximation, the thermal conductivity of the wood during heating can be treated as effectively constant. Consequently, we utilize a constant, species dependent, bulk thermal conductivity for models discussed in this paper. While we only considered dry wood (no moisture content) for purposes of the present discussion, relationships that account for the effect of moisture content and/or particle size changes on bulk heat capacity and thermal conductivity are available for a variety of wood species.34,52 2.2. Assumptions About Particle Shape and Size. External particle shape is important because realistic wood particles are typically far from spherical, and the nature of the external surface determines the available area for heat exchange with the surroundings. Detailed image analysis of milled pine feedstock18 have revealed that most particles have an elongated shape as depicted in Figures 1 and 3. For purposes of simulation, we utilized typical aspect ratios obtained from the image analyses of milled pine as functions of the particle Feret diameter (DF). In this study, the Feret diameter is the longest distance

− n · (− k∇T ) = h(T∞ − T )

(4)

where n is the unit normal vector at the surface, h is the convective heat transfer coefficient, and T∞ is the external temperature. As explained above, average bulk values were assumed for k based on our previous study.18 Finite element simulations were performed using COMSOL 4.3b with the CFD and CAD Import modules. Solutions of these transient simulations were obtained using a fully coupled, direct solver with a relative error tolerance of 1 × 10−3. 2.4. Computational Simulation of 1-D Heat Conduction. Transient 1-D intraparticle heat conduction was simulated using a general 1-D form of the heat equation that has been applied to slabs, cylinders, and spheres.53 1 ∂ ⎛⎜ b ∂T ⎞⎟ ∂T kr = ρC p b ⎝ ⎠ ∂ r ∂ r ∂t r

(5)

where r is the 1-D spatial coordinate (m), b is the shape factor (0 slab, 1 cylinder, 2 sphere), T is temperature in kelvin (K), k is thermal C

DOI: 10.1021/acs.energyfuels.6b00554 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels conductivity (W/m·K), ρ is density (kg/m3), Cp is heat capacity (kJ/ kg·K), and t represents time in seconds (s). We assumed the same values for thermal conductivity and density as for the 3-D simulation, the variable heat capacity from eq 2, and no shrinkage or expansion of the particle. Boundary conditions include symmetry at the particle center (r = 0, eq 6) and convective heat transfer at the particle surface (r = R, eq 7). ∂T =0 ∂r k

at

r=0

∂T = h(T∞ − TR ) ∂r

at

(6) r=R

(7)

An implicit finite-difference approximation scheme was applied to eqs 5−7 to numerically integrate them with variable thermal properties 54

[1 + 2(1 + b)Fo]T0n + 1 − 2(1 + b)FoT1n + 1 = T0n

(8) Figure 4. Analytical and numerical solutions of transient heat conduction in a solid sphere. Lines represent surface (Ts) and center (Tc) temperature profiles from numerical model (num) while symbols represent profiles from analytical solution (ana).

⎛ ⎛ b⎞ b⎞ − Fo⎜1 − ⎟Tin−+11 + (1 + 2Fo)Tin + 1 − Fo⎜1 + ⎟Tin++11 = Tin ⎝ ⎠ ⎝ 2i 2i ⎠ (9)

⎡ ⎛ b ⎞⎤ n + 1 n+1 ⎜ ⎟⎥T − 2FoTM − 1 + ⎢1 + 2Fo 1 + Bi + Bi M ⎝ ⎣ 2M ⎠⎦ ⎛ b ⎞⎟ n T∞ = TM + 2FoBi⎜1 + ⎝ 2M ⎠

Table 2. Characteristic Length Scales for Irregularly Shaped Particles (10)

where M is the numerical spatial index corresponding to the particle surface, where r = R. The center particle temperatures are represented by eq 8, bulk average internal particle temperatures by eq 9, and surface temperatures by eq 10. The Fourier number is Fo = αΔt/Δr2, thermal diffusivity is α = k/ρCp (m2/s), the Biot number is Bi = hr/k, and T∞ represents the external temperature in kelvin (K). The subscript i denotes spatial position within the particle measured incrementally from the center (i = 0 at r = 0) to the surface (i = M at r = R). Time increments are represented by the n superscript where current temperature Tn advances in time to temperature Tn+1 at a time step of 0.01 s. Note that Fo and Bi can vary with temperature, providing solutions for the more general case where the biomass thermal properties are changing during heat up. The system of discretized equations developed from eqs 8−10 can be solved in matrix form A·x = b where A is the coefficient matrix, x is the unknown temperature vector at the next time step, and b is the known temperature vector at the present time. The coefficient matrix is a tridiagonal banded matrix that was solved in Python using the SciPy linear algebra function “scipy.linalg.solve_banded”. The analytical solution for 1-D transient heat conduction in a solid sphere with constant properties (eq 11) was used to verify the 1-D numerical model developed from eqs 8−10. Figure 4 illustrates the agreement for the surface and center temperature profiles for a spherical particle with d = 1 mm, ρ = 540 kg/m3, h = 350 W/m2·K, k = 0.12 W/m·K, and Cp = 1800 J/kg·K. ∞

θ=

∑ Cnexp(−ζn2Fo)Dn , n=1

Dn =

1 sin(ζnr ), ζnr

Cn =

diameter

description

volume equivalent surface area equivalent surface area to volume

sphere diameter with same volume as particle sphere diameter with same external surface area as the particle sphere diameter with same external surface area to volume ratio as the particle, DSV = 6DCH ratio of particle volume to surface area overall particle height minimum particle dimension largest dimension orthogonal to DH longest distance between pair of parallel lines tangent to the largest projected particle outline

DS = (S/π)1/2 DSV = (D3V/D2S) DCH = V/S DH DW DL DF

characteristic height width length Feret diameter

A range of particles with Feret diameters (DF) of 200 μm to 20 mm were simulated at fast pyrolysis conditions for loblolly pine (softwood) and white oak (hardwood). The 3-D solid particle geometry was constructed from realistic wood particles then imported into Comsol for transient heat conduction simulation. The intraparticle temperature profiles obtained from the three-dimensional Comsol simulation were compared to the low-order model. Tables 3 and 4 summarize the parameters used for the 1-D and 3-D simulations for two representative hardwood and softwood particles at dry conditions (0% moisture content, x = 0). The value of the convective heat transfer coefficient, h, used for both the 1-D and 3-D simulations, was selected from the range of values expected for fast pyrolysis in bubbling fluidized beds,26 as were the initial and final particle temperatures, To and T∞, respectively.10,36

4[sin(ζn) − ζncos(ζn)] , 2ζn − sin(2ζn)

1 − ζncot(ζn) = Bi

symbol DV = (6V/ π)1/3

3. RESULTS AND DISCUSSION 3.1. Comparisons of 1-D and 3-D Heat Conduction. To confirm the consistency of the 1-D and 3-D heat conduction computations, we compared the predictions of both methods for the same spherical particle. Figure S1 in the Supporting Information illustrates the resulting temperature profiles from each approach for a 1 mm diameter sphere evaluated at ρ = 540 kg/m3, k = 0.12 W/m·K, Cp = 1800 J/kg·K, and h = 350 W/m2· K. The initial temperature of the solid sphere was 293 K with a surface exposed to an ambient temperature of 773 K. Excellent agreement was observed between the center, surface, and

(11)

To systematically explore various approaches for defining a characteristic length scale to use in our 1-D simulations, we considered a range of proposed equivalent diameters for nonspherical particles as listed in Table 2.55−57 Temperatures at various locations in the 3-D particle model were compared to temperature profiles (center, surface, and average) obtained from the 1-D solution. Figure 5 (left) displays the locations on the biomass particle model where temperatures were logged during the Comsol simulation. Several different characteristic lengths (R = D/2) were evaluated based on the list summarized in Table 2. D

DOI: 10.1021/acs.energyfuels.6b00554 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

Figure 5. Locations of 3-D temperature profiles (left) and dimensions evaluated as a characteristic length or equivalent spherical diameter (right).

changing b alone includes the implicit assumption that the simulated particles have infinite size in the other dimensions. To illustrate this point further, average temperature profiles generated from 3-D simulations of spherical and cubic particles having the same volume and the same constant properties used above are compared in the right frame of Figure 6. In this case, we note that a spherical particle heats more slowly, because the external surface area to volume (or mass) ratio is smaller for the sphere than the cube. The higher external surface area is of course characteristic of all nonspherical particles and is why cryogenic storage tanks are designed to be spherical to minimize heat transfer with the surroundings.58 This has important implications for biomass pyrolysis because the particles are typically nonspherical. Since the degree of nonsphericity in biomass particles can be affected by both biomass type and the milling method, this makes particle shape an important factor to be considered in process simulations. The variation of heating rate with surface to volume ratio suggests that the surface area to volume diameter (DSV) might be the appropriate characteristic length scale to use for approximating particle heat transfer in 1-D. In the following section, we illustrate the results of using this and other characteristic length scales with 1-D approximations to evaluate how well different choices can account for nonspherical particle shapes. 3.3. Comparison of 1-D and 3-D Transient Temperature Profiles. To compare 1-D temperature profiles for nonspherical particles with 3-D profiles, we used results from a previous study of 3-D heat conduction in milled white oak particles. The 3-D study included heat conduction simulations with explicit intraparticle microstructure and demonstrated that microstructural impacts could be accounted for in heat transfer simulations using isotropic bulk properties as long as appropriate particle shapes were used.18 The specific properties and physical dimensions assumed for heat conduction simulations here are summarized in Tables 3 and 4. Seven of the eight different characteristic length scales (as defined in

Table 3. Bulk Thermal Properties and Model Parameters used for 1-D and 3-D Simulations of a Representative Hardwood (White Oak) and Softwood (Loblolly Pine)a property

Loblolly Pine

White Oak

ρ (kg/m ) k (W/m·K) h (W/m2·K) Cp (J/kg·K) To (K) T∞ (K)

540 0.12 350 103.1 + 3.867T 293 773

720 0.16 350 103.1 + 3.867T 293 773

3

Values for ρ and k along with equation for Cp were obtained from the Wood Handbook.34 a

volume-averaged temperatures and demonstrates that the simulation methods converge as expected in the limiting spherical case. 3.2. Transient Heat Conduction in Non-Spherical Particles. A 1-D simulation of nonspherical particles requires application of the shape factor (b in eq 5); however, as we explain below, this must be done with care to properly account for the controlling heat transfer factors, especially the particle volume (mass) and external surface area. It is tempting to simply adjust the shape factor b to simulate the heating of rectangular and cylindrical particles.53 Unfortunately, changing the shape factor alone is not sufficient, because the resulting description no longer takes into account the roles of volume (mass) and external surface area. We believe this has led to some confusion in the literature, as illustrated in the left frame of Figure 6, which replicates comparisons made by Babu et al.30 of temperature profiles with 1-D heat conduction in rectangular, cylindrical, and spherical particles produced by adjusting only the shape factor. A casual comparison of these profiles might lead to the conclusion that spherical biomass particles should heat faster than cylindrical or rectangular particles. This is actually an artifact resulting from the fact that

Table 4. Physical Dimensions of CAD Geometry Used for 1-D and 3-D Simulations of Loblolly Pine and White Oak Particles Feret Diameter 0.2 mm 0.4 mm 0.7 mm 1.4 mm 2.8 mm 5.4 mm 10 mm 20 mm

volume (m3) −13

8.895 × 10 5.553 × 10−12 2.11 × 10−11 8.442 × 10−11 4.011 × 10−10 2.877 × 10−9 1.827 × 10−8 1.462 × 10−7

surface area (m2) 5.355 1.879 4.836 1.394 4.614 1.716 5.885 2.354

× × × × × × × ×

10−8 10−7 10−7 10−6 10−6 10−5 10−5 10−4 E

height (m)

width (m)

1.8914 × 10−4 3.8505 × 10−4 6.8678 × 10−4 0.001409 0.002839 0.005475 0.01014 0.02027

6.9766 × 10−5 1.1694 × 10−4 1.6084 × 10−4 2.0012 × 10−4 2.6284 × 10−4 5.069 × 10−4 9.387 × 10−4 0.001877

length (m) 8.8923 × 1.6267 × 2.5184 × 3.9492 × 7.0966 × 0.001369 0.002534 0.005069

10−5 10−4 10−4 10−4 10−4

DOI: 10.1021/acs.energyfuels.6b00554 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

Figure 6. Temperature profiles for sphere, cylinder, slab, and cube shapes. (left) Average temperatures produced by the 1-D numerical model by only changing the shape factor for a characteristic length of 1 mm. (right) Surface and center temperatures from 3-D simulations of spherical and cubic particles of the same volume but different surface area to volume ratios.

Figure 7. Simulated temperature profiles for white oak particles at two different particle sizes with Feret diameters of 200 μm (left) and 20 mm (right). The symbols denote 3-D results for the specific particle locations indicated in Figure 5. The solid lines represent volume average temperatures from 1-D simulations assuming different characteristic lengths. Time step remained constant in 3-D simulation but not all data points are displayed in the figure.

Figure 8 compares the 1-D predictions for the surface, centerline, and volume averaged temperatures in a loblolly pine particle with corresponding 3-D simulations that accounted for realistic geometry. The 1-D model utilizing the surface area to volume diameter as the characteristic length provided the best agreement with the center, surface, and volume average temperature profiles exhibited by the 3-D Comsol simulation for a particle size of DF = 5.4 mm. To further assess the impact of particle size, we expanded the low-order DSV and 3-D comparisons to include a range of white oak particle sizes: DF = 200 μm, 400 μm, 700 μm, 1.4 mm, 2.8 mm, 5.4 mm, 10 mm, and 20 mm with the parameters listed in Tables 3 and 4. Figure 9 depicts the results for the volume average temperature profiles up to DF = 20 mm. We observe that the volume average (TV) profiles appear to agree well up to DF = 5.4 mm. Beginning at DF = 10 mm, the 1-D TV profile starts to noticeably deviate from the 3-D profile but remains in good agreement with the 3-D simulations. 3.4. Linking Particle Heat Conduction with Pyrolysis Reaction Kinetics. As originally noted by Pyle and Zaror45 and more recently by Paulsen et al., 59 the rate of devolatilization in pyrolysis involves a balance between internal and external heat transfer and between heat transport and

Table 2) were used for the 1-D simulations to determine which gave the best agreement with the 3-D temperature profiles. The Feret diameter (DF) was not included because of its similarity to DH. 1-D simulations for small (DF = 200 μm) and large (DF = 20 mm) oak particles based on the surface area to volume diameter (DSV) were found to give the best agreement with the 3-D simulations for volume average temperatures as perceived in Figure 7. We saw similar good agreement between the 1-D DSV model and 3-D predictions for volume average temperatures in the case of a loblolly pine particle with the properties and dimensions listed in Tables 3 and 4. For the large particle size (Figure 7 right), there was more disparity between the temperatures (TV, TST, TC, TL, TW, TS) recorded at various locations on the 3-D particle. We conjecture that the later offset is a consequence of the increased temperature gradients experienced by larger particles, as heat transfer becomes an even more rate limiting step than it is for smaller particles. Even with the difference in intraparticle temperature profiles for larger particle sizes, the volume-average temperature TV from the 3-D simulation can be reasonably reproduced with the loworder model when DSV is used as the equivalent spherical diameter. F

DOI: 10.1021/acs.energyfuels.6b00554 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

an appropriate temperature for evaluating the rate constant and thermal properties, it is also necessary to select the length scale in the above dimensionless groups. The effects of assuming different possible length scales are illustrated by the symbols in Figure 10 for an example case where the thermal properties of white oak in Table 3 are assumed, along with a reaction rate coefficient of K = 0.0622 s−1 based on the work of Sadhukhan et al.47 Note that different assumptions for the characteristic length place the example case in very different pyrolysis regimes on the map. Based on the results presented above, we suggest that R = DSV/2 would probably be the best choice for the characteristic length since it reflects the proper external surface to volume properties of the particle. 3.5. Adding Heat of Reaction. In some instances, the heat of reaction is sufficiently large that it can also affect the particle heat up. An example of this is illustrated in Figure 11, which is based on the experimental measurements of Sadhukhan et al.47 For simulating this case, we coupled the 3-step kinetic scheme from Sadhukhan et al. with the 1-D heat transfer model based on DSV and the parameter values summarized below in Table 5. As indicated by Figure 11, reaction heat causes the center particle temperature to overshoot the surrounding temperature. Given the importance of reaction heat (which might be negative or positive in some cases at certain temperatures), two other dimensionless numbers appear to be relevant to consider: PyIII = k(T∞ − To)/(ρR2KΔH) and PyIV = h(T∞ − To)/ (ρRKΔH) where ΔH is the heat of reaction released by pyrolysis. These groups represent, respectively, the ratio of conductive heat transport to reaction heat generation and convective heat transport to reaction heat generation. As above, we propose that the most appropriate characteristic length scale to use in these ratios is R = DSV/2.

Figure 8. Comparison of 3-D temperature profile results for a loblolly pine particle with 1-D results using the particle surface area to volume diameter (DSV) as the characteristic dimension. Dimensions and properties for a DF = 5.4 mm particle size are listed in Tables 3 and 4. Temperatures are evaluated at the surface Ts, center Tc, and volume average Tv.

reaction. These balances can be expressed in terms of two types of dimensionless numbers: a Biot number (Bi = hR/k), which reflects the ratio of convective vs conductive heat transport, and two pyrolysis numbers represented by PyI = k/(ρCpR2K) and PyII = h/(ρCpRK). The former reflects the ratio of heat transfer by conduction to the rate of heat loss in the evolved products leaving the particle. The latter is a similar ratio that reflects the effect of convective heat transfer to the external particle surface relative to the rate of heat loss in the reaction products. In this case, h, k, Cp, and ρ are defined in Table 3, while R is a characteristic particle length scale and K is a reaction rate constant that must be evaluated at some selected temperature.45,59,60 Using the above dimensionless groups, it is possible to construct a regime map that identifies four distinct regions of behavior, depending on which process is rate controlling, as depicted in Figure 10. Such maps can be useful for recognizing cases in which one particular heat transfer mode is dominant, because it is then possible to focus on the controlling process and thereby reduce computational complexity. Besides selecting

4. CONCLUSIONS The results presented here underscore the importance of using a characteristic length scale based on the external surface area to volume ratio of biomass particles when simulating heating during pyrolysis. This can be challenging because of the range of shapes in realistic biomass feedstocks and the lack of information about the different particle dimensions involved in defining shape that are not available from sieve analysis alone. In many instances, it may be important to account for a range

Figure 9. Comparison of low-order (DSV) temperature profiles (solid lines) to 3-D temperature profiles (TV volume average) for white oak particle. Models evaluated at Feret diameters of 200 μm−2.8 mm (left) and 5.4−20 mm (right). Time step remained constant in 3-D simulation but not all data points are displayed in the figure. G

DOI: 10.1021/acs.energyfuels.6b00554 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

Figure 10. Pyrolysis regimes according to dimensionless Biot (Bi) and pyrolysis numbers (PyI and PyII) for white oak particles at T = 773 K. (left) Effect of characteristic length for particles with DF of 200 μm (▲) and 20 mm (●). (right) Particles characterized by DSV where DF = 200 μm to 20 mm.

Table 5. Particle Properties and Kinetic Parameters Used to Model Pyrolysis of Wood Cylinder As Demonstrated by Sadhukhan et al.47a

Figure 11. Center temperature profiles for a 20 mm × 100 mm wood cylinder at 683 K. Symbols represent experimental data provided by Sadhukhan et al.47 The solid blue line (Tc_dsv) denotes results from the 1-D model with ΔH = −240 000 J/kg; the solid red line (Tc_dsv, H = 0) represents the 1-D model with no heat of reaction, ΔH = 0.

parameters

value

diameter (mm) height (mm) ρ (kg/m3) k wood (W/m·K) k char (W/m·K) Cp wood (J/kg·K) Cp char (J/kg·K) h (W/m2·K) To (K) T∞ (K) A1, A2, A3 (1/s) E1, E2, E3 (J/mol) ΔH (J/kg) δ (−)

20 100 682 0.13 + 0.0003(T − 273) 0.08−0.0001(T − 273) 1112 + 4.85(T − 273) 1003 + 2.09(T − 273) 48 285 683 168.4, 13.2, 5.7 × 106 51965, 45960, 92400 −240000 1.38

a

Pre-factor is represented as A, activation energy, as E, overall heat of reaction, as ΔH, and deposition coefficient, as δ.

On the basis of the results summarized in this paper, we propose the following guidelines for low-order simulation of biomass particle heat up under fast pyrolysis conditions: • The surface area to volume equivalent spherical diameter, represented as DSV, offers a good approximation of the characteristic length scale needed to approximate heat conduction in realistic woody biomass particles with onedimensional simulations. • Effective (bulk averaged) values for thermal conductivity and heat capacity, such as those provided by the Wood Handbook,34 can be utilized in low-order models to approximate heat conduction in realistic woody biomass particles, even though the particles can be highly anisotropic. • The dimensionless Biot and pyrolysis numbers are useful for assessing the relative effects of convective and conductive heat transport and reaction kinetics and for identifying regions of behavior where each of these effects is dominant; however, they do not account for heat of reaction.

of particle shapes in widely distributed feedstocks. Such factors are probably most important when different types of biomass and different milling techniques are involved, which is a likely scenario that will be encountered in commoditized biomass feedstocks pyrolyzed at the industrial scale. Simulations of conductive heat transfer based on onedimensional spherical approximations can yield relatively accurate predictions of the transient intraparticle temperatures inside realistic, three-dimensional, nonspherical particles if the appropriate geometry and bulk thermal properties are used. The external surface area to volume ratio of particles is a key geometric factor, since it determines the available area per unit mass through which heat can enter the particle. Even though the material properties of biomass particles are inherently anisotropic, our findings suggest that bulk average properties can provide acceptable results in one-dimensional, particle-scale simulations because they already account for directional averaging in realistic particles. H

DOI: 10.1021/acs.energyfuels.6b00554 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels • When heat of reaction is significant, modified dimensionless pyrolysis numbers that include the heat of reaction should be utilized. We expect that the above approach for approximating particle heat transfer in one-dimension will be useful in reducing the computational overhead involved in simulating biomass fast pyrolysis in multiphase flow reactors. This is important because the computational demands for modeling other processes such as hydrodynamics and chemical reaction kinetics can be extremely high. While the results presented here are focused primarily on intraparticle conductive heat transfer, it seems clear that implementation of these simplified heat conduction models can be a starting point for more sophisticated particle-scale models of pyrolysis that include the effects of reaction kinetics, particle size distributions, and extra-particle transport.

■

(3) Bridgwater, A. V. Biomass Bioenergy 2012, 38, 68−94. (4) Di Blasi, C. Chem. Eng. Sci. 1996, 51 (7), 1121−1132. (5) Di Blasi, C. Chem. Eng. Sci. 2000, 55, 5999−6013. (6) Okekunle, P. O.; Pattanotai, T.; Watanabe, H.; Okazaki, K. J. Therm. Sci. Technol. 2011, 6 (3), 360−375. (7) Gaston, K. R.; Jarvis, M. W.; Pepiot, P.; Smith, K. M.; Frederick, W. J.; Nimlos, M. R. Energy Fuels 2011, 25 (8), 3747−3757. (8) Kersten, S. R. a.; Wang, X.; Prins, W.; van Swaaij, W. P. M. Ind. Eng. Chem. Res. 2005, 44 (23), 8773−8785. (9) Mellin, P.; Kantarelis, E.; Yang, W. Fuel 2014, 117, 704−715. (10) Kersten, S.; Garcia-Perez, M. Curr. Opin. Biotechnol. 2013, 24 (3), 414−420. (11) Alves, S. S.; Figueiredo, J. L. Chem. Eng. Sci. 1989, 44 (12), 2861−2869. (12) Babu, B. V.; Chaurasia, A. S. Energy Convers. Manage. 2004, 45 (9−10), 1297−1327. (13) Bharadwaj, A.; Baxter, L. L.; Robinson, A. L. Energy Fuels 2004, 18 (4), 1021−1031. (14) Bryden, K. M.; Ragland, K. W.; Rutland, C. J. Biomass Bioenergy 2002, 22 (1), 41−53. (15) Chan, W.-C. R.; Kelbon, M.; Krieger, B. B. Fuel 1985, 64 (11), 1505−1513. (16) Lu, H.; Ip, E.; Scott, J.; Foster, P.; Vickers, M.; Baxter, L. L. Fuel 2010, 89 (5), 1156−1168. (17) Chaurasia, A. S.; Babu, B. V. Influence of Product Yield, Density, Heating Conditions and Conversion on Pyrolysis of Biomass; Proceedings of International Conference on Desert Technology-7 (DT-7), November 2013; pp 9−14. (18) Ciesielski, P. N.; Crowley, M. F.; Nimlos, M. R.; Sanders, A. W.; Wiggins, G. M.; Robichaud, D.; Donohoe, B. S.; Foust, T. D. Energy Fuels 2015, 29 (1), 242−254. (19) Di Blasi, C. AIChE J. 2002, 48 (10), 2386−2397. (20) Galgano, A.; Di Blasi, C. Combust. Flame 2004, 139 (1−2), 16− 27. (21) Greenhalf, C. E.; Nowakowski, D. J.; Harms, A. B.; Titiloye, J. O.; Bridgwater, A. V. Fuel 2013, 108, 216−230. (22) Haseli, Y.; van Oijen, J. a.; de Goey, L. P. H. Chem. Eng. J. 2011, 169 (1−3), 299−312. (23) Janse, A. M. C.; Westerhout, R. W. J.; Prins, W. Chem. Eng. Process. 2000, 39 (3), 239−252. (24) Yuen, R. K. K.; Yeoh, G. H.; de Vahl Davis, G.; Leonardi, E. Int. J. Heat Mass Transfer 2007, 50 (21−22), 4371−4386. (25) Koufopanos, C. A.; Papayannakos, N.; Maschio, G.; Lucchesi, A. Can. J. Chem. Eng. 1991, 69, 907−915. (26) Papadikis, K.; Gu, S.; Bridgwater, A. V. Fuel Process. Technol. 2010, 91 (1), 68−79. (27) Yashwanth, B. L.; Gallacher, J.; Shotorban, B.; Mahalingam, S.; Fletcher, T. H.; Weise, D. R. Experimental and numerical investigation of the effect of heating modes and moisture content on pyrolysis and ignition of live fuels. In U.S. National Combustion Meeting, Cincinnati, Ohio, May 17−20, 2015; pp 1−14. (28) Wang, X.; Kersten, S. R. A.; Prins, W.; van Swaaij, W. P. M. Ind. Eng. Chem. Res. 2005, 44 (23), 8786−8795. (29) Himmel, M.; Tucker, M.; Baker, J.; Rivard, C.; Oh, K.; Grohmann, K. Biotechnol. Bioeng. Symp. 1985, No. 15, 39−58. (30) Babu, B.; Chaurasia, A. Energy Convers. Manage. 2003, 44 (14), 2251−2275. (31) Papadikis, K.; Gu, S.; Bridgwater, A. V. Fuel Process. Technol. 2010, 91 (7), 749−758. (32) Di Blasi, C. Int. J. Heat Mass Transfer 1998, 41 (24), 4139− 4150. (33) Bonnefoy, F.; Gilot, P.; Prado, G. J. Anal. Appl. Pyrolysis 1993, 25, 387−394. (34) Glass, S. V.; Zelinka, S. L. In Wood Handbook: Wood as an Engineering Material; Forest Products Laboratory: Madison, WI, 2010; pp 4.1−4.19. (35) Ragland, K.; Aerts, D.; Baker, A. Bioresour. Technol. 1991, 37, 161−168. (36) Prakash, N.; Karunanithi, T. Asian J. Sci. Res. 2009, 2 (1), 1−27.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.energyfuels.6b00554. Table of thermal properties found in biomass pyrolysis literature, table comparing volume equivalent shapes, and a figure plotting conversion of a 6 mm wood cylinder (PDF) An overview of the Computational Pyrolysis Consortium (CPC) can be found at http://cpcbiomass.org/. The Python code developed for the low-order intraparticle model is available on GitHub at https://github.com/pyrolysis in the “low-order-particle” repository.

■

AUTHOR INFORMATION

Corresponding Author

*Phone number: 865-946-1340. E-mail address: wigginsg@ ornl.gov. Notes

Notice: This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC0500OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for the United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doepublic-access-plan). The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Computational Pyrolysis Consortium (CPC) funded by the U.S. Department of Energy, Bioenergy Technologies Office (BETO).

■

REFERENCES

(1) Basu, P. In Biomass Gasification, Pyrolysis and Torrefaction: Practical Design and Theory; Academic Press: San Diego, CA, 2013; pp 47−84, 147−176. (2) Butler, E.; Devlin, G.; Meier, D.; McDonnell, K. Renewable Sustainable Energy Rev. 2011, 15 (8), 4171−4186. I

DOI: 10.1021/acs.energyfuels.6b00554 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels (37) Di Blasi, C. Combust. Sci. Technol. 1993, 90, 315−340. (38) Galgano, A.; Di Blasi, C. Ind. Eng. Chem. Res. 2003, 42, 2101− 2111. (39) Galgano, A.; Di Blasi, C.; Horvat, A.; Sinai, Y. Energy Fuels 2006, 20, 2223−2232. (40) Gronli, M. G.; Melaaen, M. C. Energy Fuels 2000, 14, 791−800. (41) Haseli, Y.; van Oijen, J. a.; de Goey, L. P. H. J. Anal. Appl. Pyrolysis 2011, 90 (2), 140−154. (42) Larfeldt, J.; Leckner, B.; Melaaen, M. C. Fuel 2000, 79, 1637− 1643. (43) Melaaen, M. C. Numer. Heat Transfer, Part A 1996, 29, 331− 355. (44) Papadikis, K.; Gu, S.; Bridgwater, A. V.; Gerhauser, H. Fuel Process. Technol. 2009, 90 (4), 504−512. (45) Pyle, D. L.; Zaror, C. A. Chem. Eng. Sci. 1984, 39 (1), 147−158. (46) Rath, J.; Steiner, G.; Wolfinger, M.; Staudinger, G. J. Anal. Appl. Pyrolysis 2002, 62 (1), 83−92. (47) Sadhukhan, A. K.; Gupta, P.; Saha, R. K. Bioresour. Technol. 2009, 100 (12), 3134−3139. (48) Kung, H.-C. Combust. Flame 1972, 18 (2), 185−195. (49) Miller, R. S.; Bellan, J. Combust. Sci. Technol. 1997, 126, 97−137. (50) Xiong, Q.; Kong, S.-C.; Passalacqua, A. Chem. Eng. Sci. 2013, 99, 305−313. (51) Xiong, Q.; Aramideh, S.; Kong, S. Energy Fuels 2013, 27, 5948− 5956. (52) Kretschmann, D. E. In Wood Handbook: Wood as an Engineering Material; Forest Products Laboratory: Madison, WI, 2010; pp 1−46. (53) Bergman, T. L.; Lavine, A. S.; Incropera, F. P.; Dewitt, D. P. Fundamentals of Heat and Mass Transfer, 7th ed.; John Wiley and Sons, Inc.: Hoboken, NJ, 2011. (54) Ozisik, M. N. In Finite Difference Methods in Heat Transfer; CRC Press: Boca Raton, FL, 1994; pp 249−257. (55) Jennings, B. R.; Parslow, K. Proc. R. Soc. London, Ser. A 1988, 419 (1856), 137−149. (56) Allen, T. In Powder Sampling and Particle Size Determination; Elsevier, 2003; pp 56−136. (57) Trottier, R.; Dhodapkar, S. CEP Mag. 2014, No. July, 36−46. (58) Mital, S. K.; Gyekenyesi, J. Z.; Arnold, S. M.; Sullivan, R. M.; Manderscheid, J. M.; Murthy, P. L. N. Review of Current State of the Art and Key Design Issues With Potential Solutions for Liquid Hydrogen Cryogenic Storage Tank Structures for Aircraft Applications, NASA/TM2006-214346; NASA Glenn Research Center: Cleveland, OH, 2006 (59) Paulsen, A. D.; Mettler, M. S.; Dauenhauer, P. J. Energy Fuels 2013, 27 (4), 2126−2134. (60) Mettler, M. S.; Vlachos, D. G.; Dauenhauer, P. J. Energy Environ. Sci. 2012, 5 (7), 7797−7809.

J

DOI: 10.1021/acs.energyfuels.6b00554 Energy Fuels XXXX, XXX, XXX−XXX