Release Profile of Volatiles in Fluidised Bed Combustion of Biomass

There is an increasing interest for power production using biomass fuels. These fuels are considered as CO2 neural, since CO2 released during burning is again captured by growing plants. When applying CCS technology along with biomass combustion a net removal of CO2 from atmosphere takes place [1]. Fluidised bed combustion has long been applied as a flexible means to burn biomass. Oxy-fuel circulating fluidized bed combustion technology has been under development [2]. The combustion is achieved in a mixture of oxygen and recirculated flue gases instead of air and the flue gases consist of carbon dioxide and water vapour allowing CO2 capture. Especially the use of waste products from wood and food industry is beneficial, since no extra land area is required. Such fuels are plentiful available, but they are dispersed in wide area. There are two ways to utilize this biomass resource. It can be co-combusted in a large boiler plant using coal as the main fuel [36]. The other option is to burn the fuel in smaller plants producing electric energy. The efficiency for generation of electricity in smaller plants is lower, but in cold climate the waste heat can be utilized in the local heat production. The combined heat and power is widely applied in Finland using local wood and peat as auxiliary fuel. In hot climate the waste heat could be used for producing cooling by the absorption heat pump system. In fact, this technology is applied even in cold climate in Helsinki in Finland for district cooling of office buildings.


Introduction
There is an increasing interest for power production using biomass fuels. These fuels are considered as CO 2 neural, since CO 2 released during burning is again captured by growing plants. When applying CCS technology along with biomass combustion a net removal of CO 2 from atmosphere takes place [1]. Fluidised bed combustion has long been applied as a flexible means to burn biomass. Oxy-fuel circulating fluidized bed combustion technology has been under development [2]. The combustion is achieved in a mixture of oxygen and recirculated flue gases instead of air and the flue gases consist of carbon dioxide and water vapour allowing CO 2 capture. Especially the use of waste products from wood and food industry is beneficial, since no extra land area is required. Such fuels are plentiful available, but they are dispersed in wide area. There are two ways to utilize this biomass resource. It can be co-combusted in a large boiler plant using coal as the main fuel [3][4][5][6]. The other option is to burn the fuel in smaller plants producing electric energy. The efficiency for generation of electricity in smaller plants is lower, but in cold climate the waste heat can be utilized in the local heat production. The combined heat and power is widely applied in Finland using local wood and peat as auxiliary fuel. In hot climate the waste heat could be used for producing cooling by the absorption heat pump system. In fact, this technology is applied even in cold climate in Helsinki in Finland for district cooling of office buildings.
Fluidised bed combustion technology is advantageous for various fuels [3][4][5][6][7][8] since it is not always necessary to process the waste to more refined fuels but it can usually be burned directly. However, processing the fuel for example by pelletizing before combustion could be beneficial [9]. Drying of fuel with low temperature waste or solar heat is means to increase the availability of energy from biomass by increasing the heating value. Biomass fuels contain a large faction of volatile matter, which is even higher than the proximate value due to high heating rate in the furnace, which increases the part of fuel released as volatiles and reduces the amount of char. In this paper a simplified method to predict release profile of volatiles in fluidized bed (FB) and circulating fluidized bed (CFB) boilers is presented.
From boiler manufacturers' point of view it is important to know burning profiles for different fuels. This enables an optimal dimensioning and placement of heat exchange surfaces in a furnace. Volatile release from fuel is divided between the dense bed and the riser. The burning of the volatiles and heat release take place only after the volatiles are released from the particles. Combustion and gas concentration profiles affect also the temperature profile and emissions of CO, NO x and SO x in fluidized and circulating fluidized beds [10][11][12]. The knowledge of the profile of the release of volatiles gives means to optimize primary, secondary and tertiary air locations and rates. CO emissions can become high if much of the volatiles are released high above the dense bed leaving too short time for oxidation. This could be avoided for specific fuels by the design of the combustor (gas velocity) and by the processing of the biomass to optimal size or by pelletizing to suitable size and density. There are several studies [7,8,9,[13][14][15][16][17][18] concerning combustion of high volatile fuel in fluidized or circulating fluidized beds. Circulation in CFB has been studied [19,20].
There is a wide literature concerning the computational methods for calculating the flow of particles in fluidized beds [21,22], but they require heavy computing and the accuracy is not so good due to the fact that flow phenomena are complex and biomass fuels are rather heterogeneous by size and shape [23]. A simplified method is developed here to predict volatile release profiles. The prediction is based on estimating the physical location and distribution of fuel particles and release of volatiles depending on particle size after feeding a fuel batch into the boiler. The final release profile is found by summing up the contributions from different times and from different size fractions.

Heating, Drying and De-volatilization of a Biomass Particle
The studies on de-volatilization times of coal particles in fluidizedbeds have been reviewed [24]. Effects of moisture u, particle shape Page 2 of 12 has been presented [29]. In some models the biomass is divided into different compounds with different reaction kinetics. Some models developed for the kinetics of solid fuels can predict the composition of the gases. The single step Arrhenius-type reaction is generally applied to describe the rate of pyrolysis with a constant final density ρ f . and Arrhenius type rate coefficient . Even this model is frequently used; it is physically not valid, because ρ f depends on some final temperature of the particle. The particle does not reach this temperature immediately after injection to the furnace, but pyrolysis takes place during the heat-up period to a great extent. During heat-up the particle cannot "predict" its future final temperature, which depends on the conditions of the ambient (which can even be changed with time). For coals a correlation has been presented [30,31] for the final density or yield of volatiles depending on final temperature. There is wide variation between the different values of pre-exponential coefficient and activation energy E measured by different researchers [32][33][34][35][36] with different heating rates and conditions ( Figure 1). The higher activation energy for coal [37] shows that pyrolysis takes place at higher temperatures.
The distributed activation energy model accounts for reactions taking place on a broader temperature range for coals. Parameters for this model for wood have been reported [38]. The model of Kobayashi [39,40] for coal consists of two competing reactions. This model has the advantage that it can predict char yield depending on the heating rate.
A model, in which the driving force is the difference between the instantaneous density and the final density at that temperature, ρ-ρ f (T p ), instead of commonly applied difference between instantaneous density and the final density at some indefinite future temperature T f , ρ-ρ f (T p ), better describes the real situation [41]. In this case the rate of pyrolysis is modelled by Eq. (1) but with temperature dependent final density. The temperature changes with time depending on the heating conditions. The temperature dependent final density [42] can also be approximated by a logistic function as shown in Figure 2 with values C≈0.03/K and T ref ≈616K. The fraction released as volatiles is v 0 ≈0.75…0.8 for wood at relative low heating rate. The amount of volatiles v 0 depends on the heat-up rate of the particle and with high heating rate v 0 is higher even 0.9-0.95 [43]. This heating rate depends on particle size and ambient temperature (gas, walls). The reaction rate coefficient for Eq. (3) with temperature dependent ρ f has been reported [42] and its dependence on temperature is milder than that by the usual exponential Arrhenius rate equation, because the strong temperature dependency is now accounted for by the final density, which is not a constant but based on the prevailing instantaneous temperature. Different models have been developed for de-volatilization of single particles (models with spatially uniform temperature and density, shrinking core models, models for particles with non-uniform temperature). In the model assuming uniform temperature and density the particle temperature is calculated by the equation If the effect of heat of pyrolysis is assumed insignificant, Eq. (4) is becomes (sphericity) φ s , size d and ambient temperature on devolatilization and burning of wood particles have been studied [25,26]. Equations of the form have been presented for the de-volatilization time of a wood particle. f(T) is an experimental function and n and C are experimental constants. The knowledge of the total time of de-volatilization is not enough, when considering the location, where the volatiles are released. The de-volatilization as function of time can be presented by experimental correlation for example by suitable cumulative distribution. The Rosin-Rammler type has been applied [27]. Pyrolysis does not start immediately after feeding to the boiler and here a time lag t d due to drying and heating up (to about 473 K) is added to shift the zero point of the time scale. The parameters t d , τ v and n can be found from batch combustion or pyrolysis tests.
More fundamental calculations require the knowledge of the kinetics of pyrolysis. Several models for the kinetics of pyrolysis of biomass have been presented [28]. Kinetics data for different biomass    = + Γ is assumed constant, the particle temperature T p can be solved In reality the time constant depends on temperature and particle density, which change with time. By eliminating time t Eqs. (1) and (5) This equation can be integrated [44]. The result in slightly different form is where Θ= RT/E and E 1 is exponential integral. This is illustrated in Figure 3.
In reality the particle density and specific heat capacity and heat transfer coefficient change with time. There may also be a flame around the particle increasing the effective temperature of the ambient [45][46][47]. Then a numerical method can be applied to calculate the particle temperature and mass as function of time [45]. In practice the temperature of the particle is also not uniform. The internal heat transfer resistance of the particle can approximately be taken into account by the coefficient φ r [45,47] in the effective heat transfer coefficient For complicated particle shapes it can be found numerically [48]. The non-uniform temperature affects also on the average pyrolysis kinetics and a correction coefficient has been presented [47].
In reality there will be temperature differences even inside small pyrolysing particles. Then more accurate approach is to account for non-uniform temperature distribution. The energy equation for the particle including storage, conduction and convection of drying and pyrolysis products is This equation is coupled with Eq. (3) for the local pyrolysis inside the particle. Numerical methods to solve the equations have been presented (see e.g. [41]).
Wood fuels usually have higher moisture content that coal. The drying delays devolatilization and the processes become overlapping [41] for large particles. This is illustrated in Figure 4.
The application of the detailed model, Eq. (10), is cumbersome as a sub-model. Drying and pyrolysis take place in narrow temperature ranges inside a biomass particle leading to steep moisture and density distributions even for relatively small particles [46]. Simplified shrinking core models for drying [49,50] and simultaneous drying and pyrolysis [51] have been presented.

Division of Combustion of a Particle in the Dense bed and above it
Drying may affect the release profile of volatiles inside the boiler in some cases. If the particle is completely dried before its weight reduces to the limit where its free settling velocity becomes below the gas velocity, moisture has no effect on the release profile, because the particle stays in the dense bed. A smaller particle may reach this limit or can even initially be below the limit. Then drying may have a great influence on the release profile. The heat up of the particle also affect the release profile of the volatiles if the particle mass is below the critical mass. For small particle heating up and drying take place before significant pyrolysis takes place and much of the volatiles can be released above the bed. For a large particle drying and pyrolysis can be overlapping [41]. Then more volatiles can be released in the bed, since moister particle is also heavier. However, particles may be completely dried before pyrolysis and before reaching the critical weight and then drying will have no effect on the release profile. Models can be used to estimate, if drying and pyrolysis are overlapping [41,50].
The mass of particle decreases during drying and devolatilization and its size reduces during char combustion. A model for division of the release profiles of volatiles and char combustion in the bed and above it [17] is further developed here. It was based on the free settling velocity of a single particle. If the free settling velocity is lower than the fluidization velocity, the particle remains in the dense bed. If it is higher, then the particle may escape the bed, but does not necessarily do so.
The movement of particles in the dense bed is rather chaotic. Above the bed the concentration of solids decreases [52,53]. Then it is possible to estimate the relation between diameter and density of particles that will be entrained with the flow, if the gas above the dense bed is assumed to be dilute enough. Gravity will make the particles with high density to go back to bed, but small light particles escape the bed. The propensity of the particles to burn in the bed or in the riser is estimated by using the equation of motion of a single particle. A border between  Figure 3: Calculated pyrolysis conversion of particles during devolatilization as function of (a) particle temperature (b) time. d=0.1, 1 and 10 mm (τ p =0.0087, 0.41 and 6.40 s), these regions corresponds to situation, in which the free settling velocity of the particle is equal to the gas velocity ( Figure 5). Thus a relation between the critical particle (volume average) density and size is obtained between the regions. The location of the border line depends on the gas velocity in the riser. So it is possible to estimate the proportions of drying, pyrolysis and char burning that take place in the dense bed at minimum, but in reality the proportions can be higher, since the particles may stay in dense bed for some time before the escape even the free settling velocity is lower that the fluidization velocity.

Burning Regions (Dense Bed and above Dense in Riser)
A particle first dries (blue line) and releases volatiles (red line) with decreasing density shown by the horizontal lines ( Figure 5). It is possible to see from such a figure the part of de-volatilization that will take place in the bed region at minimum, since heavy particles cannot easily leave the bed. However, small and light particles do not necessary leave the bed, but their burning can take place also in the dense bed. Then some analysis about their residence time in the bed is required. The borderline is determined by the free settling velocity of the particles and it is in reality only a rough approximation for the border of the dense bed and dilute region, since the shape of the fuel is not spherical and the particles do not behave as single particles. The border is also not abrupt, but its shape and location depend on conditions. The devolatilization and char burning may also be overlapping especially for small particles [45]. It is seen (Figure 5a) that most of the burning of a large heavy particle takes in the bed and only after the mass is decreased they can escape the bed.
For large particles (Figure 5a) drying and de-volatilization take place completely in the bed regions, since the burning (red) line crosses the (green) free settling velocity curve only in the char combustion stage. Much of char burning also takes place in the dense bed (because char mass ~ d 3 ). The char after crossing the critical intersection point (in the region "burning in circulation") can partly be burned in the dense bed and partly in the riser. Part of volatiles is released above the dense bed for lighter large particles and all volatiles from small particle can be released above the bed (Figure 5b).
After fuel is fed into a boiler drying de-volatilization and char combustion take place. Whether a particle escapes the bed or not depends on its free settling velocity, which depends on size and mass of the fuel particle, which change due to drying, pyrolysis and char oxidation. Heavy particles fall back into the dense bed while the lighter ones are entrained by the gas flow. However, even the free settling velocity of the particles is lower than the gas velocity statistically it takes some time that they reach the boundary region between the dense bed and the riser and eventually escape the dense bed. This can be described by the elutriation rate constant or its inverse, the elutriation time constant. Most of the volatiles are released in the bed for particles residing long in the bed, but for smaller particles with shorter residence time, much of the de-volatilization can take place in the riser.
The particle velocity in an upward flow is described by the equation There is a great diversity of shapes of biomass particles. Coal particles are more regular in shape. There exists a vast literature on drag coefficient of particles of different shapes in different conditions [54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71]. In the present paper the drag coefficient is calculated from a correlation for a spherical particle [72]  when Re≤1000 and C D =0.44 when Re>1000, but the calculation method can be applied to other correlations for C D as well as distributed values of C D . The particle does not move, if the free settling velocity is equal to the gas velocity. In this case dU p /dt=0 and U p =0. Equation (11) gives the explicit relation between critical density and critical particle size at different gas velocities 2 , 3 Re Thus, it is possible to get information about the minimum fraction of fuel released from particles (as volatiles or char combustion products) in the dense bed. In reality the fraction can be larger, since the particles are not immediately escaping the dense bed when critical density and size are reached. If de-volatilization takes place with a constant size and decreasing density, the minimum fraction of volatiles released in bed is when v 0 ρ p0 >ρ p , c and d 0 >d c . When v 0 ρ p0 <ρ p , c and d 0 >d c , X C,c =0 char burning can take place even totally above the dense bed.
The above assumptions, de-volatilization with constant particle size and char combustion with constant density [17], are frequently applied in combustion modelling of solid fuels. In reality particles may shrink and oxidation takes place also inside the particle. Then the burning routes in reality are curved instead of the direct lines as shown in Figure  6 due to shrinkage during pyrolysis and internal combustion during char combustion. For some coals swelling also takes place affecting particle size and shape and its movement.
The volume change due to shrinkage during de-volatilization is approximated to depend linearly on degree of conversion where s is constant.
The mass of a particle is The particle mass decreases during de-volatilization and The PSD of the coal can be described by Rosin-Rammler type distribution for a size interval d min …d max . Then the PSD and cumulative size are The PSD becomes the normal Rosin-Rammler distribution, with d min =0 and d max =∞. Eqs. (15)-(17) give the relation for particle density during devolatilization In the right hand size term the numerator accounts for pyrolysis and the denominator for shrinkage or swelling. Both the size and density can change during char oxidation [73,74] Minimum 19.7 % of volatiles is released and 1.41 % of char is oxidized in the dense bed with the PSD given.

Mass Balance Model for Dense Bed
We consider the de-volatilization stage of particles. A fuel batch with particles of single size is fed in the CFB reactor. It may be considered as a differential batch in a continuous feed so that the batch itself does not affect the steady state environment in the reactor. The system under study is illustrated in Figure 8. The change of volatile mass in particles in the dense bed due to this batch is described by The term on the left hand side is the change of mass. The first term on the right hand side is the feed. The second term is the mass flow rate of volatiles in elutriated particles m K m e e = 0  , the third term is due to mass loss by de-volatilization , the fourth term is mass loss due to removal of bottom ash . The last term accounts for the mass flow of volatile mass in particles returned into the bed. Devolatilization can take place completely in the bed, if the free terminal velocity U t is higher than the gas velocity (i.e. the elutriation coefficient remains zero, K e =0).
The particle residence time in the riser part of the loop t r1 depends on the initial density and particle size. Thus the mass flow rate in the riser changes due to de-volatilization and it is described by can be removed also from the loop seal in addition to the bottom ash removal. It is assumed that no reactions take place in the flow down, but they could also be included by using longer residence time instead of t r1 . The mass flow rate of partly pyrolysed particle back to the reactor is The fact that the particles are retuned with an average time delay depending on the prevailing past bed mass is shown by the time difference r t t − in the bed mass.
We use notations . The time constant for elutriation depends on the particle size and density. Particle density decreases during de-volatilization. Particle size may also change due to shrinkage (or swelling for some fuels).
We follow the subsequent behaviour, when a batch with initial mass m 0 is fed in the boiler at time t=0. The bed mass can be solved numerically as function of time The mass volatile matter fed in the batch of particles is divided in different parts with different fates. Volatile matter released as volatiles inside the bed (i=v), volatile matter in particles removed as bottom ash (i=b), volatile matter (inside particles) removed as fly ash (i=f) and volatile matter in particle removed in the loop seal (i=s) are

Elutriation of Particles from the Dense Bed
Elutriation is described by the equation [52] There is a great number of correlations for elutriation rate constant in the literature. Some of these correlations are summarized [75][76][77][78][79][80][81]. As an example some of the equations for the elutriation coefficient can be presented in the form where c i are constants presented in Table 1. Other more suitable correlations based on measurements or results from comprehensive flow models could also be used.

Particle Motion in the Riser
The residence time of the particles is obtained by solving the equation of motion, Eq. (11) for a single particle. Eq. (11) can also be presented as using the time constant If the time constant is assumed constant or average value is applied, the solution is for large times t>>τ w t g p U U U − ≈ (38) where the terminal velocity can be calculated by iteration. The particle size decreases due to the shrinkage according to Eq. (15) giving . Density decreases due to pyrolysis but it is also affected by the shrinkage as described by Eq. (18).
The terminal velocity will change during pyrolysis, since particle size and density are changing C D Re also depends on the changing particle size (especially for large particle). The velocity of the particle can be solved numerically using the equation However, if the time constant or relaxation time τ w is small, the particle reaches its equilibrium velocity t g p U U U − = fast, which can be directly applied. In addition, the starting velocity of particles from the dense bed varies. Some particles are accelerated after leaving the bed and some particles have already excess velocity and are decelerated. This distribution of initial velocities is not well-known and it gives some dispersion for the volatile release due to particle-to-particle differences. In this context it should be noted that particles of same size and density have not in reality the same residence time in the riser, but there is a residence time distribution (RTD) (see e.g. [82][83]). The gas does not flow in plug-flow manner, but there is a velocity profile affecting also to the particle velocities in the riser. Part of the particles may flow downwards in near the wall region. The particles are also not of the same shape and oriented similarly, which gives some particle-toparticle differences in values of C D . Thus, the release profile in the riser that is found here for one particle type is smoothened in reality due to this RDT. The analysis at different locations could be carried out with a measured or modelled RTD [84], but the weight and size reduction of the particles complicate the use of such methods.
The particle velocity U p =dx/dt gives the dependence between location of particle (distance from the dense bed) and time Solving the particle location and release of volatiles as function of time after it is removed from the bed by elutriation poses a difficulty in finding the steady state release profile of volatiles in the riser based on the unsteady state response of a batch. The particles are reducing in mass and size and their velocity is changing. This means that x-values calculated with the above equation will change. It would be better to find the particle velocity and release of volatiles directly as function of co-ordinate x. Using relation U p =dx/dt equation (36) can be presented in the form which can be presented in the form for numerical calculations , .
, , Then the corresponding time is found by When the fine and coarse particles were mixed, the carry-over rate was increased significantly. A simple equation and calculation procedure has been presented [85]. The use of the effective gas density have been recommended [14] (1 ) g g s ρ ερ ε ρ = + − for calculation motion of large biomass particles in flow of fine particles.

Discussion
The method is illustrated in the following for a pilot CFB reactor. Particles with size 1.5 mm are fed to the boiler. In this case the drying occurs completely in the bed (see Figure 6d). The gas velocity is assumed to be 3 m/s. The heat up time t d and de-volatilization time constant τ v depend on particle size (~d n , n≈ 1.5). In the illustration heating up

Page 8 of 12
and drying time are given. Here we assume that t d =0.5 s (which has no effect in this case on the release profile) and τ v =0.5 s. These values depend on particle size and can be estimated from literature or from bench scale de-volatilization tests. A linear decrease of moisture during this initial period is assumed. Particle density is described simple by equation during heat up and drying. The time lag t d accounts for heat up and drying time before pyrolysis starts. The particle diameter and density are changing simultaneously with the flight due to de-volatilization as described by Eqs. (15) and (18). In the illustration de-volatilization is assumed here to obey Eq. (2) with n=1, but more comprehensive models to calculate pyrolysis and mass as function of time could also be used. No fragmentation is assumed during pyrolysis but it the effects could be included provided that the sufficient information is available.
The values of Tanaka et al. [81] (Table 1) are applied for elutriation. In this illustration it is assumed that the total bed mass is 30 kg/m 2 . The release of volatiles in the bed is illustrated in Figure 9. The mass (volatiles in particles) starts to decrease at t=0.5 s after the particles have been dried and heated up. Elutriation starts at t=1.115 s, when the particles have lost weight enough so that U t <U g and some of the particles are entrained by the gas. The rates of devolatilization and elutriation are also shown in Figure 9. Elutriation rate reaches a maximum (at t=1.4 s) when the particles have lost enough weight but still contain volatile matter. It should be noted here that m is the mass of volatiles. Char and ash are not included but their effect on particle mass and elutriation is taken into account. It is seen that about 80 % of volatiles is released in the bed in this case and the rest is released in the riser.
The subsequent distance above from the dense bed for these particles (elutriated at t=1.115 s) is shown as function of time in Figure  10a. If the devolatilization is considered to be competed at 3 s, the riser height < 2 m (not including the bottom bed) is enough. Velocities of particles are shown in Figure 10b. The particle velocities were calculated using two methods (Eq. (36) black line, Eq. (38) red line), which give practically the same line (the black line is under the red line). It takes some time before the particle velocity reaches the velocity U g -U t in this case with relative large particles. Release profile for volatiles is shown in Figure 10c.
The situation in the riser (in this batch study) changes with time, since the elutriated particles become lighter when pyrolysing in the bed. Figure 11 shows the situation for particles entrained at t=1.4 s. In this case with large particles the curves are rather similar as earlier, since the particles have released much weight already before elutriation. Small particles can be elutriated already even before pyrolysis starts so that the profiles in the riser will change with time due to the increasing conversion of the inlet particles to the riser.
Thus, the release profile in the riser changes (for the batch process) as can be seen by comparing Figures 10 and 11. The conversion as function of location in the riser is shown in Figure 12 at times t=1.115, 1.415, 1.715, 2.015 and 2.315 s. The inlet conversion to the riser increases with time. The average (steady state) conversion in the riser and its derivative are shown in Figure 12b. The profile is steep close

Page 9 of 12
to bed because the velocity of particles close to bed is low, since their mass is still close to the critical one. It is practically the same as the conversion at 1.415 s at which the elutriation reaches maximum. This average conversion was calculated by

Conclusions
Analytical formulas for the (minimum) amounts of volatiles released and char oxidized in the dense bed are presented. The effect of shrinkage of particles was included. The effect of fragmentation can be considered in further development.
A model for the release profile of volatiles in FB and CFB boilers has been presented. It gives the division of the release of volatiles in the dense bed and in the riser. It is based on the fact that large and heavy particles do not escape the dense bed before complete de-volatilization. Smaller particles also reside some time in the dense bed releasing volatiles and are entrained by the gas from the dense bed the bed to the riser. This elutriation depends on the terminal velocity of the particle. Particle size and density decrease also due to de-volatilization. Method to calculate the release profile of the volatiles in the riser is also presented. The propensity for elutriation increases with decreasing particle size and density. The analysis is carried out for all size fractions and the final profile is found by summing up each contribution. The model gives also the possibility to get information of the release distributions of fuel nitrogen and sulphur from particles during de-volatilization stage. The extension of the method to include char combustion stage is possible.

Acknowledgment
This work was performed with funding from the project ERA-NET Biomodelling-Advanced biomass combustion modelling for clean energy production within the framework of ERA-NET Bioenergy Joint Call on Clean Biomass Combustion.

Appendix: Analytical Solutions for Bed Mass
Here we consider simplified solutions with constant values for the coefficients K. In reality especially K e is not constant but depends on the changing size and density of the particle. K e is, however, constant for an inert particle. This inert case can be useful in evaluating elutriation time constant using cold reactor tests without devolatization or using tests at higher temperatures with inert material. By applying the Laplace transform