Stresses in a Cylindrical Wood Particle Undergoing Devolatilization in

Mar 18, 2008 - ... Indian Institute of Technology Madras, Chennai-600 036, India ... A mathematical model capable of estimating the temperature and st...
0 downloads 0 Views 2MB Size
Energy & Fuels 2008, 22, 1549–1559

1549

Stresses in a Cylindrical Wood Particle Undergoing Devolatilization in a Hot Bubbling Fluidized Bed M. Sreekanth,* B. V. S. S. S. Prasad, and Ajit Kumar Kolar Department of Mechanical Engineering, Indian Institute of Technology Madras, Chennai-600 036, India

H. Thunman and B. Leckner Department of Energy and EnVironment, Chalmers UniVersity of Technology, SE-412 96 Göteborg, Sweden ReceiVed NoVember 5, 2007. ReVised Manuscript ReceiVed January 13, 2008

A mathematical model capable of estimating the temperature and stress distribution in a devolatilizing cylindrical wood particle is proposed. The model assumes wood as orthotropic and comprises a thermal submodel and a stress submodel. It takes into account the contributions to stress by the temperature rise, shrinkage during drying and devolatilization, and mechanosorption. The thermal submodel estimates the temperature distribution and shrinkage as a function of space and time to be used as inputs to the stress submodel. Calculations have been carried out for cylindrical wood particles of aspect ratio l/d ) 1 and with dimensions of 10, 20, and 30 mm. The model predictions compared with the results of an available drying model are found to be satisfactory. During devolatilization of the wood particle, the model predicts high compressive stresses at the center while those at the surface are tensile in nature and moderate in magnitude. A modified von Mises criterion predicts failure at the axial center of the cylinder, with the size of the crack formed by the failure increasing from 0.4 to 1.4 mm in diameter as the size of the wood particle increases from 10 to 30 mm. The predicted crack radius and the time of crack initiation are compared with the measurements carried out on wood particles undergoing devolatilization in a fluidized bed combustor. The predicted crack radius and time of crack initiation are found to be sensitive to the magnitude of the mechanical properties.

Introduction The estimation of stresses in a devolatilizing wood particle in a combustor is essential in assessing its structural integrity. Highly stressed wood particles could fragment into smaller particles whose reduced size would be the effective size for subsequent char combustion. A fragmented wood particle would have its surface area increased manifold, thereby accelerating all the processes dependent on the surface area while also contributing to elutriation out of a fluidized bed or dropping through the grate openings in a fixed bed, resulting in fixed carbon loss. Stress models for wood have been reported for applications in wood drying1–3 or for studying the structural stability of construction timber.4 Although drying is an important precursor to devolatilization, in a hot fluidized bed combustor it lasts for a short while and the magnitudes of temperature and shrinkage involved are small in comparison to those of devolatilization. Therefore, the stresses during devolatilization could be higher than those during drying. There are very limited stress models for drying of cylindrically shaped wood.3 Stress models are most commonly reported for rectangular boards. However, there is no model that predicts the stresses developed in a devolatilizing cylindrical wood * To whom correspondence should be addressed. Telephone: 91-44-2257 5664. Fax: 91-44-2257 4652. E-mail: [email protected]. (1) Tauchert, T. R.; Hsu, N. N. Wood Sci. Technol. 1977, 11, 51–58. (2) Pang, S. Drying Technol. 2000, 18 (8), 1677–1696. (3) Kang, W.; Lee, N. H. Wood Sci. Technol. 2002, 36, 463–476. (4) Fu, K. C.; Harb, A. Comput. Struct. 1987, 27 (2), 225–235.

particle. Though some work has been done on stress development in devolatilizing coal,5–7 it cannot be directly extended to wood due to the structural differences (wood is orthotropic while coal is considered isotropic in thermal stress models) and the differences in the structural changes that take place during the process (wood shrinks while some coals swell, though it was not considered in the models referred to5–7). Moreover, coal is generally considered to be spherical in shape, which rarely can be assumed for wood. The present model aims at estimating the stresses in a devolatilizing cylindrical wood particle by extending the available stress models for wood drying. The model treats wood as an orthotropic material and includes the stresses due to temperature rise, shrinkage during drying and devolatilization, and mechanosorption. Knowledge of these stresses, along with the pressure of volatiles inside the wood particle and the impact of the bed material, is necessary for predicting its structural integrity in a hot fluidized bed. Experimental Section Experiments were conducted to understand the nature of the fracture that occurs in a devolatilizing wood particle in a fluidized bed combustor. Cylindrical Casuarina equisetifolia wood particles with sizes ranging from 10 to 30 mm were subjected to devolatil(5) No, S. Y.; Syred, N. J. Inst. Energy 1990, 63, 195–202. (6) Stanmore, B. R.; Brillard, A.; Gilot, P.; Delfosse, L. Proc. Combust. Inst. 1996, 26, 3269–3275. (7) Dacombe, P.; Pourkashanian, M.; Williams, A.; Yap, L. Fuel 1999, 78, 1847–1857.

10.1021/ef700658k CCC: $40.75  2008 American Chemical Society Published on Web 03/18/2008

1550 Energy & Fuels, Vol. 22, No. 3, 2008

Sreekanth et al.

Table 1. Proximate and Ultimate Analyses of C. equisetifolia quantity proximate analysis moisture volatiles fixed carbon ash ultimate analysis carbon hydrogen nitrogen oxygen

content (%) 10.8 72.5 16.4 0.3 42.5 6.1 0.16 51.24

ization under oxidizing conditions in a laboratory scale bubbling fluidized bed combustor at a bed temperature of 850 °C. The proximate and ultimate analyses of the wood are shown in Table 1. The reactor has a diameter of 130 mm and a height of 600 mm and is electrically heated from the outside. A wood particle of known dimensions and mass was introduced into the bubbling bed with the help of a basket made of stainless steel mesh. The particle was retrieved from the bed at the end of devolatilization, which was identified by the extinction of the flame. The char particles obtained were quenched by using sand at room temperature. After cooling, the char fragments were observed for their nature of fracture. A detailed discussion of experimental details and results can be found in the work of Sreekanth et al.8 Figure 1 shows the nature and progress of the fracture. It can be seen that primary fragmentation takes place during the second half of the devolatilization time. Also, fragmentation is frequent and extensive for larger particles. The following observations, relevant to the present work, have been made from the experiments. Most of the time, wood cracked at the radial center with the cracks radiating outward for a few millimeters, as shown in Figures 2–5. Cracks were also seen circumferentially and 2–3 mm inside the periphery. These cracks became increasingly noticeable as the devolatilization progressed, which was observed during the study on the timing of fragmentation. The crack radius varied from less than 1 to 5 mm from the smaller to the larger sizes, and its formation was noticed during the first half of the devolatilization time. Also, the diameter of the wood particle at the midpoint of the length became a minimum due to formation of a neck. At the end of devolatilization, an unfragmented core remains, with fragments falling off from near the periphery. Formulation and Development of the Model. The model involves estimation of the temperature profile in the wood particle and particle shrinkage, as functions of time and space. These temperatures and shrinkages are further used to compute the stress. The estimated stress is then evaluated according to a suitable criterion to predict fracture/cracking in the wood material. The model comprises a thermal submodel and a stress submodel. Cylinders having an aspect ratio l/d ) 1 are chosen for the present study. Such cylinders cannot be treated as infinite cylinders, and therefore, the thermal model is formulated in two dimensions. The stress model aims at predicting the failure, its timing, and its radial location. Hence, a 1-D stress model in the radial direction is chosen for simplicity. Thermal Submodel. The following assumptions are made in the formulation of the thermal submodel: 1. The wood particle is cylindrical in shape and free of bark. 2. The properties of the solid fuel vary linearly between those of wood and char, depending on the degree of conversion. 3. The physical structural integrity is retained; that is, the cylindrical shape is retained even though shrinkage takes place. 4. The volatile gases and water vapor escape instantaneously from the wood without accumulation. Therefore, the convective term in the conservation equations is neglected. (8) Sreekanth, M.; Kolar, A. K.; Leckner, B. Proceedings of the 19th International Conference on Fluidized Bed Combustion, Vienna, Austria; Winter, F., Ed.; Rutzky Druck: St. Poelten, Austria, 2006; Paper No. 97.

5. The wood particle is assumed to shrink by 20% and 10% in the radial and longitudinal directions, respectively, as measured.9 6. Effects due to the condensation of water vapor in the cooler regions of the wood particle and cracking of tar (secondary reaction) to form gaseous volatiles are negligible. 7. The external heat transfer coefficient is uniform around the particle surface. The reaction scheme for devolatilization, shown in Figure 6, consists of (i) moisture evaporating to water vapor and wood reacting to yield (ii) gaseous volatiles, (iii) tarry volatiles, and (iv) char. The reactions are endothermic in nature, and the heat consumed is considered in the model. Governing Equations in the Thermal Submodel. The mass conservation equations are as follows: (a) Conservation of wood ∂F ˜w ) w˙w ∂τ where w˙w is the rate of consumption of wood. (b) Conservation of char

(1)

˜c ∂F ) w˙c ∂τ where w˙c is the rate of generation of char. (c) Conservation of moisture

(2)

∂F ˜m ) w˙m (3) ∂τ where w˙m is the rate of moisture loss. (d) The unsteady energy conservation equation in 2-D cylindrical coordinates is given by

(F˜wCp

w

+F ˜cCpc + F ˜mCpm)

∂T 1 2 ∂ ∂T ∂ ∂T ) f f rk + fr2fz k + ∂τ r r z ∂r r ∂r ∂z z ∂z

(

)

( )

4

∑ (KF˜H)

i

(4)

i)1

where the term on the left-hand side of eq 4 represents energy storage. On the right-hand side, the first term represents the heat conduction in the radial direction, the second term represents the heat conduction in the longitudinal direction, and the last term represents the heat generation (consumption in this case) within the particle. The shrinkage factors fr and fz along the radial and longitudinal directions, respectively, are defined as the ratio of the instantaneous dimension to the initial dimension. The shrinkage factors are estimated based on the degree of conversion. Their values vary from 1 to 0.8 in the radial direction and from 1 to 0.9 in the longitudinal direction, as measured by Kumar et al.9 Shrinkage, defined as the ratio of reduction in size to the initial size, can be obtained by subtracting the shrinkage factors from 1. The radius, r, is the instantaneous radius of the particle and is computed at the beginning of every time step. Interpolation Factor. For the properties related to the solid structure, a linear variation between the virgin wood and the char is assumed as follows: φ ) ηφw + (1 - η)φc

(5)

where the interpolation factor is η ) Fw /Fw0

(6)

Effective Thermal Conductivity. The thermal conductivity of the solid undergoing conversion varies with the direction of the grain, temperature, and degree of conversion. The effective thermal conductivity is (9) Kumar, R. R.; Kolar, A. K.; Leckner, B. Biomass Bioenergy 2006, 30 (2), 153–165. (10) Grönli, M. G.; Melaaen, M. C. Energy Fuels 2000, 14, 791–800. (11) Chan, R.; Kelbon, M.; Krieger, B. Fuel 1985, 64, 1505–1513. (12) Blasi, C. D. J. Anal. Appl. Pyrolysis 1998, 47, 43–64. (13) Bryden, K. M.; Hagge, M. J. Fuel 2003, 82, 1633–1644.

Stresses in a Cylindrical Wood Particle k ) kcond + krad

Energy & Fuels, Vol. 22, No. 3, 2008 1551 (7)

composed by the total conductivity kcond ) ηkw + (1 - η)kc

(8)

and the radiant conductivity term10 krad )

4εg σεdporT 3 (1 - εg)

(9)

where dpor ) ηdpor,w + (1 - η)dpor,c

(10)

by Bryden and Hagge,13 so that the model predicts drying around 373 K. The heat of the drying reaction is taken to be the latent heat of vaporization of water. The kinetic constants are given in Table 3. Initial and Boundary Conditions. The energy and mass conservation equations are solved by using the following initial and boundary conditions: at t ) 0, T(r,t) ) T0; Fw(r,t) ) Fw0; Fc(r,t) ) 0; Fm(r,t) ) Fm0 (13) at r ) 0,

∂T ∂T ) 0 and at r ) R, kr ) h(Tbed ∂r ∂r Ts,r) for all z (14)

at z ) 0,

∂T ∂T ) 0 and at z ) Z, kz ) h(Tbed ∂z ∂z Ts,z) for all r (15)

Specific Heat Capacities. The specific heat capacity of the solid undergoing conversion is given by Cp,s ) ηCp,w + (1 - η)Cp,c

(11)

The thermophysical data used in the model10 are shown in Table 2. Kinetic Constants. The reactions are assumed to follow the Arrhenius expression, Ki ) Ai exp(-Ei/RT)

(12)

No data regarding the kinetic constants obtained for C. equisetifolia (the wood used in the referenced experiments) under fluidized bed combustion conditions are available in the literature. So the data obtained under TGA conditions11 are employed to predict the qualitative behavior of wood pyrolysis, as they are not associated with any one particular wood species.12 The pre-exponential factor for the drying reaction (reaction i)11 has been modified, as suggested

The total heat transfer coefficient is considered to be the sum of the convective and radiative components. The convective component depends on the fluidization conditions and on the sizes of the fuel and the inert bed material, and its magnitude increases as the fuel particle shrinks. The radiative component depends on the temperature difference between the bed and the fuel particle surface. The bed-particle surface heat transfer coefficient is calculated by using Palchonok’s correlation.14 The inert bed material in the fluidized bed is assumed to be sand having a mean size of 550 µm and at 850 °C, fluidized by air at a superficial velocity 5 times the minimum fluidizing velocity of 0.1 m/s. The convective heat transfer coefficient is calculated by using the following correlation:

Figure 1. Sequence and mode of fragmentation of wood during devolatilization: (a) photograph showing wood particles (20/20, 25/25, and 30/30 mm) before devolatilization and at the end of one-quarter, half, three-quarters, and complete devolatilization; (b) schematic showing the structural changes in a cylindrical wood particle during devolatilization.

1552 Energy & Fuels, Vol. 22, No. 3, 2008 Nuin ) [0.85Arin0.19 + 0.006Arin0.5Pr0.333] + [(6 +

Sreekanth et al.

[ ]

0.117Arin0.39Pr0.333) - (0.85Arin0.19 + 0.006Arin0.5Pr0.333)]

din deq

0.67

(16)

where din is the mean size of the inert bed material and deq is the equivalent sphere diameter of the wooden cylinder, based on the volume to surface area ratio. The convective heat transfer coefficient is reduced by 15%14 to account for the reduction in the heat transfer coefficient due to the floating of the wood particle on the bed. The radiative heat transfer coefficient is given by hr )

(

σ[Ts4 - Tb4] 1 1 + - 1 (Ts - Tb) εs εb

)

(17)

Figure 2. Crack at the center and weakening at the periphery of a 15/ 10 mm particle at the end of the first and second quarter of the devolatilization time.

Figure 6. Reaction scheme for the devolatilization of moist wood. Table 2. Thermophysical Properties10 property

value/expression

Cp,w Cp,c Cp,m kr,w kr,c kz,w kz,c dpor,w dpor,c Fw,dry (initial) Fm (initial) ε σ R

1.5 + (1 × 103)T 0.42 + (2.09 × 10-3)T + (6.85 × 10-7) T2 4.182 0.15 0.05 0.3 0.1 5 × 10-5 1 × 10-4 500 50 0.9 5.67 × 10-8 8.314

unit kJ/(kg K) kJ/(kg K) kJ/(kg K) W/(m K) W/(m K) W/(m K) W/(m K) m m kg/m3 kg/m3 W/(m2 K4) J/(g mol K)

Table 3. Kinetic Constants11

Figure 3. Center crack and fragmentation near the radial surface of a 20/20 mm particle at the end of the first and second quarter of the devolatilization time.

Figure 4. Center crack in a 20/20 mm particle at the end the second quarter of the devolatilization time.

Figure 5. Center crack in a 30/30 mm particle at the end of the second quarter of the devolatilization time.

reaction no. (Figure 6)

A (1/s)

E (kJ/mol)

H (kJ/kg)

i ii iii iv

5.13 × 1010 1.3 × 108 2.0 × 108 1.1 × 107

88 140.3 133.1 121.3

2244 150 150 150

The emissivities of the bed and fuel particle surfaces are taken as 0.8 and 0.9, respectively.14 The net heat transfer coefficient is typically of the magnitude of 350–550 W/(m2 K) for the particle sizes and the operating conditions of the present study. Stress Submodel. The stress model is formulated based on the principles used in the published drying models and extended to include devolatilization because of the similarity between drying and devolatilization. Both drying and devolatilization take place due to the heating up of wood, though at different temperatures. Heating results in the loss of the constituents of wood: moisture during drying and volatiles during devolatilization. Also, shrinkage takes place during drying and devolatilization. An important aspect of wood-drying is mechanosorption, which is an effect that occurs when wood under stress is subjected to a change in moisture content. Among the least understood manifestations of the behavior of wood is the interaction of moisture change with mechanical loading.15 The term mechanosorption has been coined to convey succinctly that mechanical influences combine with moisture sorption to produce a response that cannot be predicted from the response to each influence separately. For example, consider the case of a wooden block that is mechanically prevented from swelling due to moisture absorption; the final pressure required to keep it within its original dimensions is an order of magnitude lower than the pressure required to compress it after it has been allowed to swell freely.15 This means that the stress developed during the absorption is lower when under pressure (i.e., when stressed), and hence, (14) Leckner, B. Heat and Mass Transfer. In Multiphase Flow Handbook; Crowe, C., Ed.; CRC Press: Boca Raton, FL, 2006; Chapter 5.2. (15) Grossman, P. U. A. Wood Sci. Technol. 1976, 10, 163–168.

Stresses in a Cylindrical Wood Particle

Energy & Fuels, Vol. 22, No. 3, 2008 1553

Table 4. Mechanical Properties20 property value in direction property

radial

tangential

longitudinal

unit

Young’s modulus coefficient of thermal expansion Poisson’s ratio

1000 2.5 × 10-5

500 3.4 × 10–5

12 000 4 × 10–6

MPa K-1

total shrinkage mechanosorptive coefficient (k) µrt µtr

20

υzr ) 0.37; υzθ ) 0.5; υrθ ) 0.67; υθr ) 0.33; υrz ) 0.044; υθz ) 0.027 25 10 5/12

%

dσr σr - σθ + )0 (18) dr r The strain-stress relations that are developed from the relations3 defined for orthotropic materials, such as wood, in terms of moisture content contain elastic strain, shrinkage strain, and mechanosorptive strain. In the present work, these equations are modified to include thermal strain and are written in terms of shrinkages instead of moisture content in order to include shrinkage during drying and devolatilization. The modified equations are given as the following: Strain-Stress Equations. εr )

1 [σ - νrθσθ] + Rr∆T - sr - ksrσr + kµθrsθσθ Er r

(19)

εθ )

1 [σ - νθrσr] + Rθ∆T - sθ - ksθσθ + kµrθsrσr Eθ θ

(20)

1 0.5

mechanical loading and sorption need to be mathematically treated together. During wood drying, structural changes take place in the form of shrinkage, which could be caused by internal forces. Mechanosorption captures the interaction between desorption and these forces that try to keep the material intact. It would not be possible to dry wood without mechanosorptive behavior.16 From a modeling point of view, noninclusion of mechanosorption results in an overestimation of the drying stresses.3 In the present model, mechanosorption during devolatilization is taken into account because of the similarities between drying and devolatilization. The following assumptions are made in the formulation of the stress model: 1. Stresses vary along the radial direction only (i.e., the model is 1-D along the radius). This is supported by the fact that the temperature gradient and shrinkage along the longitudinal direction are low compared to those in the radial direction. The inputs from the thermal model to the stress model are taken from a plane that is midway of the length and parallel to the faces of the cylinder. 2. The state of stress is assumed to be that of plane stress. A stress model can be one of either a plane stress or a plane strain. The plane strain formulation is more complex while the differences between the results of plane stress and plane strain are negligible for drying stress applications.16 Moreover, the plane strain formulation needs negligible strain in the longitudinal direction while in reality the strain is about 0.1.9 Hence, the choice of plane stress is made in the present study. 3. The mechanosorptive strain is linearly proportional to the stress and is 5/12 times the shrinkage strain.3 Though this ratio is estimated for drying, it is assumed to also hold during devolatilization. This is motivated partly because the amount of volatiles lost (desorption) and the amount of shrinkage during devolatilization are greater than the corresponding moisture loss and shrinkage during drying, but approximately the same ratio is maintained. Another reason for making this assumption is the lack of measured data. 4. The creep strain is neglected since the times involved in the devolatilization in a fluidized bed are low, on the order of a few minutes. This is further supported by the fact that the ratio of strains due to shrinkage, mechanosorption, elasticity, and creep is roughly 12:5:1:0.2, according to Kang and Lee.3 5. The amount of shrinkage in the radial direction is available9 for cylindrical wood particles devolatilizing in a fluidized bed, whereas the amount of shrinkage in the tangential direction is not reported in the literature. Therefore, the measurements17 on birch wood (a hard wood) cubes in a single particle reactor are used to estimate the tangential shrinkage to be 1.25 times the radial shrinkage. Hence, in the present study, the tangential shrinkage is assumed to be 1.25 times that of the radial shrinkage. 6. The particle is being heated uniformly from all sides, and the shear stresses in the wood particle are negligible. 7. The mechanical properties remain constant during the process. They are given in Table 4. Governing Equations in the Stress Submodel. Equilibrium Equation. (16) Martensson, A.; Svensson, S. Holzforschung 1997, 51 (6), 472– 478. (17) Davidsson, K. O.; Petterson, J. B. C. Fuel 2002, 81, 263–270.

In both of the above equations, the first term on the right is the elastic strain based on Hooke’s law, the second term is the thermal strain, the third term is the shrinkage strain, and the last two terms are due to mechanosorption. In the plane stress formulation, the longitudinal stress is equal to zero. Strain-Displacement Equations. εr )

dur dr

(21)

ur (22) r The equilibrium (eq 18), strain-stress relations (eqs 19 and 20), and strain-displacement equations (eqs 21 and 22) are then combined to obtain the equilibrium equation in terms of radial stress only, eq 23, since the boundary conditions are in terms of radial stress. The final form of the equilibrium equation is εθ )

d2σr

a

2

dr

+b

dsθ dσr dT + cσr ) (Rr - Rθ)∆T - sr + sθ - Rθr +r dr dr dr (23)

where a) b) c)

1 2 r - ksθr2 Eθ

dsθ 3 - νθr νrθ + kµrθsrr + r and r - (3 + µθr)ksθr - kr2 Eθ dr Er

1 - νθr 1 - νrθ dsθ dsr + kµrθr + ksr(1 + µrθ) - kr Eθ Er dr dr ksθ(1 + µθr)

The expression for tangential stress is obtained from eq 18, as follows: Tangential Stress. σθ ) σr + r

dσr dr

(24)

Boundary Conditions. dσr )0 (25) dr Solution Procedure. The partial differential equations of the thermal model (eqs 1-4) along with the initial and boundary conditions (eqs 13-15) were solved numerically by using an explicit finite differencing scheme. The wood particle was initially discretized into uniformly sized control volumes. As conversion progressed, the control volumes shrank, and the new size of each control volume was calculated, depending on its degree of conversion. The temperature at a new time step was computed by using properties at the previous time step. The temperature thus obtained was used to calculate the properties at the new time step. After the At r ) R, σr ) 0, and at r ) 0,

1554 Energy & Fuels, Vol. 22, No. 3, 2008

Sreekanth et al.

Figure 7. Comparison of the predicted and measured center temperatures of a 20/4 cylinder at bed temperatures of 712 and 900 K.

Figure 9. Temperature distribution along the radius in a 30/30 wood particle as a function of time.

Figure 8. Comparison of predicted and measured center temperatures of a 20/4 cylinder at bed temperatures of 1005 and 1107 K.

Figure 10. Shrinkage along the radius in a 30/30 wood particle as a function of time.

temperature and shrinkage were obtained at the end of a time step, the stresses were calculated using eqs 23-25. From stability considerations, a time step of 0.005 s was used for the computation. Grid independence was attained when the radius and height of each initial control volume were chosen as 0.1 mm. The CPU time for the smallest wood particle (d ) l ) 10 mm) is 10 min while for the biggest one (d ) l ) 30 mm) it is 40 min on a Pentium IV, 512 MB RAM computer.

time, the temperature at the center is the lowest and that at the periphery is the highest. The reduction in the particle size due to shrinkage can be noted from the figure. The right end of every temperature profile at a given time marks the outer surface of the particle, which can be read on the abscissa. As the particle heats up, the surface recedes. The initial radius was 15 mm while the radius shrank to 12 mm by the end of 350 s, by which time the devolatilization was completed. Shrinkage in the Particle. The wood particle, on heating, undergoes shrinkage in all three directions. Shrinkage takes place due to the loss of moisture and volatiles. Moisture loss takes place when the local temperature reaches about 373 K. Shrinkage during devolatilization starts when the local temperature reaches about 650 K. In the meantime, between drying and devolatilization, the particle does not shrink. Figure 10 shows the shrinkage in the radial direction of a 30/30 mm particle at different times. It is seen that wood shrinks by 5% at the end of drying and by 20% by the end of devolatilization. The shrinkage front moves from the surface toward the center. Initially, the outer layers shrink due to drying while the inner layers remain unaffected by the heating at the surface. At an intermediate time, the outer layers shrink due to devolatilization while the inner layers shrink due to drying. At the end of 200 s, the center begins to shrink due to drying, and at the end of 350 s, the entire particle would have shrunk by 20%, in the radial direction. Density during Wood-Char Conversion. The solid is assumed to consist of wood and char. Before devolatilization commenced, it is entirely comprised of wood. As the conversion proceeds, the amount of wood decreases and that of char increases. At the end of devolatilization, the amount of wood becomes zero, and the solid is entirely char. The density of the solid varies with time, as the wood converts to char. Figure 11 shows the temporal variation of the wood density at the center

Results and Discussion Thermal Submodel. The solution of the thermal submodel gives the spatial and temporal temperature distribution, the amount of shrinkage, and the density of the converting solid material. Hereafter, a cylinder is represented by the ratio l/d. Hence, a 30/30 cylinder means l ) 30 mm and d ) 30 mm. In the following, the results of the thermal submodel are shown for a plane located midway across the length and parallel to the faces of the 30/30 cylinder. The temperature and shrinkage on this plane are input to the calculation of the stresses. Temperature Distribution in the Particle. Figures 7 and 8 show the measured and predicted center temperature in a 20/4 mm wood particle at different bed temperatures. The temperature measurements were carried out by Di Blasi and Branca18 in a bed of sand fluidized with nitrogen. The qualitative agreement between the measured and predicted temperatures was found to be good, and quantitatively, the maximum deviation was found to be 17%. This instilled confidence in the thermal model, and hence, further calculations were carried out for other sizes. Figure 9 shows the temperature distribution along the radius of a 30/30 mm wood particle at different times. At any instant of (18) Di Blasi, C.; Branca, C. Energy Fuels 2003, 17, 247–254. (19) Boresi, A. P.; Schmidt, R. J. AdVanced Mechanics of Materials, 6th ed.; John Wiley & Sons Inc.: Hoboken, NJ, 2003.

Stresses in a Cylindrical Wood Particle

Energy & Fuels, Vol. 22, No. 3, 2008 1555

Figure 11. Reduction of wood density with time at the center and surface. The time of intersection of the density at the center with the 5 kg/m3 line gives the devolatilization time (∼325 s in the figure).

Figure 13. Comparison of the present model with the predictions of Kang and Lee3 for their Model 2, Case 2 (decreasing spatial moisture gradient in the disk).

Figure 12. Comparison of the present model predictions with those of Kang and Lee3 for their Model 2, Case 1 (uniform moisture content in the disk). Also shown are radial stresses computed by the present model (by including contributions from temperature rise and shrinkage due to drying and devolatilization) at 10, 20, and 30 s in a 10/10 mm wood particle.

Figure 14. Comparison of the present model with the predictions of Kang and Lee3 for their Model 2, Case 3 (increasing moisture gradient in the disk).

and surface of a 30/30 mm wood particle. The devolatilization time, noted from the figure, is about 325 s, which is the time when the density of the dry wood at the center is reduced to 1% (5 kg/m3) of its initial value of 500 kg/m3. The center begins to devolatilize by the end of 250 s, which is marked by the beginning in the reduction of dry wood density. At the same time, the center begins to shrink due to devolatilization. Stress Submodel. Stress Model Comparison. The stress model is used to obtain the analytically computed results for a cylindrical wooden disk having predetermined moisture distributions.3 The analytical model3 considers stress due to shrinkage caused by the moisture gradient and mechanosorption. It does not consider the influence of temperature rise or drop. This makes it3 a special case of the present model. For the purpose of comparison, the first and fourth terms, which denote the influence of temperature rise and its gradient, on the right-hand side of eq 23, were removed and solved with the rest of the stress governing equations. In such a case, the present model would be able to consider only stresses due to shrinkage and mechanosorption. The shrinkage term is substituted with shrinkage only due to drying, for the sake of comparison. The present model was numerically solved, and comparisons with that of Kang and Lee3 were made for three cases having different moisture distributions. For the sake of comparison, the values from Kang and Lee3 were obtained by using the data digitization software Windig 2.5. Figure 12 compares the two models for the case where the disk has uniform moisture but varying shrinkage in the radial and tangential directions. It also shows the stresses (due to temperature rise and shrinkage due to drying and devolatilization) in a 10/10 mm particle at 10, 20, and 30 s.

Figure 15. Radial stress distribution in a 10/10 cylinder having 10% initial moisture, undergoing drying in a convective environment at 150 °C.

It can be seen that the stresses during devolatilization are very much different from those of Kang and Lee.3 This justifies the additional effort undertaken to estimate the stresses in a devolatilizing wood particle. Figure 13 compares the two models for the case where there is a decreasing moisture gradient, and Figure 14 compares them for an increasing moisture gradient. Figures 12-14 are meant to determine the predictive capability of the present model during drying. The trends match well, and the comparison of the magnitudes is satisfactory with a maximum error of 4%. Figure 15 shows the radial stress distribution in a 10/10 wood particle undergoing drying in a convective environment. The initial moisture content, ambient temperature, and external heat transfer coefficient are assumed to be 10%, 150 °C, and 50 W/(m2 K), respectively. The drying time required has been estimated as 232 s from the thermal model. The radial stresses

1556 Energy & Fuels, Vol. 22, No. 3, 2008

Figure 16. Radial stress distribution at different times in a 30/30 cylindrical wood particle undergoing devolatilization.

at the end of one-quarter, one-half, three-quarters, and the total drying times are shown. It can be seen that the radial stresses are compressive in nature, and the magnitudes are highest at the center. The magnitudes of the drying stresses obtained in the present work fall within the range of stresses predicted for a rectangular slab2 (100 × 40 × 590 mm3) dried at 150 °C. It can be concluded that the present model satisfactorily describes the drying stresses, and it can further be used to predict the stress in a wood particle during devolatilization. Stress in Wooden Cylinders. The stresses in 30/30 cylindrical wooden particles undergoing devolatilization are presented below. Figure 16 shows the radial stresses in a 30/30 wooden particle at different times. The times shown in the plots correspond to 1/ τ , 1/ τ , 3/ τ , and τ 4 dev 2 dev 4 dev dev for each size. The radial stress is compressive in nature, due to the dominance of shrinkage over thermal expansion. In all the figures, the radial stress at the center is lowest for 1/4τdev while it is maximum for 1/2τdev. The radial stresses are steep near the center, and the slope gradually decreases toward the surface. The radial stress away from the center decreases with time. This could be due to the reduction in the shrinkage rate and temperature gradients with time, away from the surface. As the devolatilization approaches completion, the stresses are reduced as the rate of shrinkage decreases and the temperature inside the particle becomes uniform, making the thermal stresses lower. At this juncture, it would be interesting to know the influence of volatiles pressure on the stresses. Volatiles pressure causes tensile stress in all directions, and their maximum magnitude can be estimated to be 1.5-2.5 bar (0.15–0.25 MPa) at the center. The pressure decreases and approaches atmospheric pressure at the surface of the particle. Compared to the maximum radial stress (9 MPa) at the center (from Figure 16), the magnitude of pressure is low. Therefore, the assumption of neglecting volatiles pressure is justified. Figure 17 shows the spatial distribution of tangential stress in a 30/30 cylinder at times corresponding to 1/4τdev, 1/2τdev, 3/ τ , and τ . At all the times shown, the tangential stress is 4 dev dev compressive at the center and tensile at the surface. This is because, as the particle shrinks, the outer layers move toward the core and exert force on it, causing compressive stresses in the core. As a reaction, the core tries to push the outer layers away, causing tensile stresses in it. Once the shrinkage ceases to occur, thermal expansion will dominate, and the nature of the stresses will be reversed to yield tensile stresses in the core while the surface has compressive stresses. The undulating nature of the curve during the initial stages is due to the occurrence of drying and devolatilization at different locations within the particle. As time progresses, this undulating nature dies out, as only devolatilization occurs at an instance. It can be noticed that the tensile stresses near the surface are well

Sreekanth et al.

Figure 17. Tangential stress distribution at different times in a 30/30 cylindrical wood particle undergoing devolatilization.

Figure 18. Stresses at the center and surface of a 30/30 cylindrical wood particle.

below the magnitudes of the compressive stresses near the center. Another observation that can be made from the figures is that the spatial point where there is a transition from compressive to tensile continuously moves from the outside toward the center as time progresses. As in the case of radial stress, the tangential stresses are steep near the center. The magnitudes of stress are almost the same for all the sizes, most likely because all sizes are assumed to undergo the same amount of shrinkage. It is observed from eqs 24 and 25 that the tangential and radial stresses are equal at the center, and the radial stress at the surface is equal to zero (eq 25). To show the temporal variation of stress, extreme locations on the particle, that is, the center and surface, are chosen. Figure 18 shows the variation of radial and tangential stresses at the center and surface in a 30/30 wood particle as a function of time. The radial stress curve shows two peaks on the compressive side. The first peak is approximately at about half-devolatilization time and corresponds to the shrinkage at the center due to drying. The stress then decreases in magnitude for a short time and again reaches a second peak, which corresponds to shrinkage during devolatilization. The maximum radial stress is about 9 MPa (compressive) for the 30/30 particle. Probable Locations of Failure and Failure Criterion. The spatial distribution of radial and tangential stresses (Figures 16 and 17) indicates a concentration of highly compressive stresses close to the center, making that area vulnerable to failure. Toward the surface, the radial stress is compressive and is continuously decreasing, whereas the tangential stress transforms from compressive to tensile in nature and is rising. This transition in the nature of stress could cause separation between the inner core and outer surface, as shown in Figure 19. Therefore, these two regions are prime candidates for failure. A suitable failure criterion is needed to predict the precise

Stresses in a Cylindrical Wood Particle

Energy & Fuels, Vol. 22, No. 3, 2008 1557

Figure 19. Nature of stresses that could cause the failure near the surface. The outward arrows show tensile stresses while the inward arrows show compressive stresses. Table 5. Ultimate Stress (in MPa) of Hardwood in Different Directions and Loadings20

Figure 20. Failure criterion function (eq 26) for a 10/10 cylinder at different locations near the center as a function of time.

nature of loading direction

tension

compression

longitudinal radial tangential (assumed same as in radial direction)

92 5 5

33 4 4

location and timing of failure. Ideally, such a failure criterion should consider the ultimate stress in different directions, the differences in the nature of loading, and also, the change in strength with time. Existing failure criteria (viz. maximum principal stress theory, maximum principal strain theory, and maximum shear stress theory19) are meant for isotropic materials and, therefore, not applicable to wood.19 The Hankinson’s formula20 relates the applied mechanical load with grain orientation and reduces to the maximum principal stress theory when the grain direction coincides with the direction of the applied load. However, it is defined for compressive loads only, making it unsuitable for the present situation. Wood has a different ultimate stress in different directions, and its strength depends on the nature of loading (tensile or compressive) too (Table 5). Hence, a modified von Mises criterion21 is used:

(

σθ σr σUr σUθ

) ( ) ( ) 2

+

σr σUr

2

+

σθ σUθ

2

-2g0

Figure 21. Failure criterion function (eq 26) for a 20/20 cylinder at different locations near the center as a function of time.

(26)

Failure occurs when the above expression is satisfied. The ultimate stress of hardwood having 12% moisture content and at room temperature in the radial direction is 5 MPa and in the tangential direction is 4 MPa.22 Figures 20-22 show the failure criterion function (left-hand side of eq 26) plotted as a function of time for locations starting from the center of the cylinder and moving outward. In Figure 20, the failure criterion is satisfied for the locations from the center to a radius of 0.2 mm. This means that there will be a crack of a diameter of 0.4 mm at the center. Beyond 0.4 mm, the failure criterion is not satisfied, and the failure propagation stops there. A similar trend can be seen for the 20/20 and 30/30 particles from Figures 21 and 22, respectively. The diameter up to which the failure occurs increases with size and reaches a value of 0.8 mm for the 20/ 20 particle and 1.4 mm for the 30/30 particle. The figures show that the criterion (and, hence, stresses) increases even after failure. In reality, the stresses should be relieved on failure. This (20) Bodig, J.; Jayne, B. A. Mechanics of Wood and Wood Composites; Van Nostrand Reinhold Company Inc.: New York, NY, 1982. (21) Pang, S. Proceedings of the 7th International IUFRO Wood Drying Conference, Tsukuba, Japan; Forest and Forest Product Research Institute: Ibaraki, Japan, 2001; pp 238–245. (22) Wood Handbook: Wood as an Engineering Material; Forest Products Laboratory, USDA Forest Service: Madison, WI, 2004.

Figure 22. Failure criterion function (eq 26) for a 30/30 cylinder at different locations near the center as a function of time.

discrepancy is due to the manner in which the model was formulated. The submodels were not coupled, and the criterion was calculated after the stresses had been computed. This has the advantage, however, to allow prediction of subsequent fractures/crack propagation, if any. For example, if the stress buildup stops when the center fails, the crack propagation from the center toward outer surface, which happens over a period of time, cannot be captured. Also, the model remains simple when not coupled. Failure time is the time from the introduction of the wood particle into the bed to the moment when failure is initiated at a radial location. From Figures 20-22, the time at which the failure occurs can also be seen, as pointed out by the arrow toward the time axis. From Figure 20, the curve of the criterion function at the center crosses the zero mark at about 11 s. The curve at 0.2 mm crosses the curve at 20 s. Hence, one could conclude that failure is initiated at 11 s and ends at 20 s. The corresponding values for the 20/20 and 30/30 particles are found in Figures 20 and 21 and are 35–60 s and 65–125 s, respectively.

1558 Energy & Fuels, Vol. 22, No. 3, 2008

Sreekanth et al.

Table 6. Crack Radius predicted crack radius (mm) size (mm)

measured radius (mm)

100% ultimate stress

75% ultimate stress

50% ultimate stress

10/10 20/20 30/30

0.7–1.3 1.6–2.9 1.5–4.6

0.2 0.4 0.7

0.4 1 1.5

0.9 1.9 2.65

Table 7. Crack Initiation Time predicted crack initiation time (s) size (mm)

measured time (mm)

100% ultimate stress

75% ultimate stress

50% ultimate stress

10/10 20/20 30/30

0–21 0–72 0–150

11 40 60

10 25 42

6 15 25

These results are subject to the assumption that the mechanical properties and ultimate stress of the wood are assumed to be constant during devolatilization. The wood particle becomes weaker during devolatilization, and the failure diameter could be greater than that predicted by the present model. Also, failure near the surface could occur as the local tangential stress changes from compressive to tensile. The location and time of the failure near the surface could be known if the dependency of strength on the time/temperature is known. To understand the dependence of crack size and time of formation on the ultimate stress, a parametric study was conducted by decreasing the ultimate stress of the wood. Table 6 shows the measured and predicted crack radius in the three wood sizes considered. These wood particles were devolatilizing in a bed at 850 °C. As the ultimate stress decreases, the crack radius increases and approaches the measured value. This could be due to the ability of the crack to propagate further away from the center due to the low ultimate stress. Table 7 shows the measured and predicted time at which a crack forms. In this case, as the ultimate stress decreased, crack formation took place earlier during the devolatilization process. These observations show that the assumption of constant mechanical properties results in conservative estimates related to the crack formation. The above discussion points to the possible locations of failure within the wood particle. Depending on the conversion process wood is subjected to, the failed wood particle could cling together or fall apart, resulting in primary fragmentation. To predict fragmentation, additional forces that come into effect, such as the impact of an inert bed on the fuel, should be taken into consideration. This along with the presented stress model will facilitate the prediction of fragmentation parameters such as the number of fragments and the size of fragments. Conclusions The following conclusions can be drawn from this study: 1. Devolatilizing wood cylinders experience compressive radial stresses throughout the particle while the tangential stresses are compressive at the center and tensile at the surface. 2. The magnitude of the stresses is almost the same for all of the sizes (-9 to 3 MPa) due to the equal amount of shrinkage undergone by each particle. 3. Radial stress shows two peaks, one due to shrinkage during drying and the other due to shrinkage during devolatilization. 4. The maximum radial stress is found to be about 9 MPa (compressive at the center) while the maximum tangential stress is 3 MPa (tensile at the surface).

5. The proposed failure criterion predicts failure in the vicinity of the center, and the diameter of the crack formed by the failure varies from 0.4 mm for the 10/10 mm particle to 1.4 mm for the 30/30 mm particle. 6. Failure time, the time from introduction of the particle in the bed to the failure, is 11 s for the 10/10 particle, 35 s for the 20/20 particle, and 75 s for the 30/30 particle. 7. The crack formation is sensitive to the mechanical properties. Hence, a model aimed at predicting primary fragmentation in a wood particle needs to consider varying mechanical properties for higher accuracy. Acknowledgment. The work presented in this paper has been financially supported by the Swedish International Development Agency (SIDA), and the authors gratefully acknowledge the same.

Nomenclature A ) pre-exponential factor (1/s) Ar ) Archimedes number Cp ) specific heat (J/(kg K)) d ) diameter of the cylinder (m) E ) activation energy (J/mol), also Young’s modulus (N/m2) f ) shrinkage factor ∆r/∆r0 ) 1 - s H ) enthalpy of reaction (kJ/kg) h ) heat transfer coefficient (W/(m2 K)) K ) rate of reaction (1/s) k ) effective thermal conductivity of a solid (W/(m K)), also the ratio of mechanosorptive strain to shrinkage strain in eqs 19 and 20 kcond ) thermal conductivity based on Fourier’s law (W/(m K)) krad ) radiant conductivity (W/(m K)) l ) length of the cylinder (m) Nu ) Nusselt number P ) pressure (N/m2) Pr ) Prandtl number r ) radial coordinate R ) outer Radius (m), also the universal gas constant (J/(mol K)) s ) shrinkage ()1 - f) T ) temperature (K), also temperature excess in the stress model u ) displacement (m) w˙ ) rate of generation (kg/(m3 s)) z ) longitudinal coordinate (m) Z ) half-height (m) Greek Letters R ) coefficient of thermal expansion (1/K) ε ) strain, also emissivity, porosity η ) interpolation factor υrθ ) Poisson’s ratio (ratio of strain in the θ direction to strain in the r direction due to stress applied in the r direction) F˜ ) uncorrected density (kg/m3), ) fr2fzF F ) corrected density (kg/m3) σ ) Stephan Boltzmann constant (5.67 × 10-8 W/(m2 K4)), also stress (N/m2) σU ) ultimate stress (MPa) τ ) time (s) µ ) dynamic viscosity of fluidizing gas ((N s)/m2), also mechanosorptive coupling coefficient in eqs 19 and 20 Subscripts 0 ) initial bed ) fluidized bed c ) char, also convective dev ) devolatilization eq ) equivalent g ) gas

Stresses in a Cylindrical Wood Particle I ) reaction number in ) inert m ) moisture por ) pore r ) radial, also radiative s ) surface of wood particle

Energy & Fuels, Vol. 22, No. 3, 2008 1559 U ) ultimate w ) wood z ) longitudinal direction θ ) tangential direction EF700658K