Heat Transfer within Direct Heating Integrated Moving Bed Torrefier by Method of Lines

Bosikawa Mosi Fabrice^*, Sumuna Temo Vertomene, Mbuyi Katshiatshia Haddy

Laboratory of New and Renewable Energy, Université de Kinshasa, Kinshasa, Democratic Republic of Congo.

DOI: 10.4236/jsbs.2021.111001                

Abstract

This work deals with the determination of the temperature profile within a direct heating moving bed torrefier in order to determine its minimum column height. A thermal model based on eulerian-eulerian two-phase solid-gas theory was performed and solved with the method of lines. In addition, this study allows to investigate the effect of the biomass particle size on the minimum torrefier column height. This investigation was performed by changing, simultaneously, the diameter of particles and the minimum fluidization velocity of the bed. Then, the calculations were made for a counter-current torrefaction reactor of 30 cm in diameter and for 5 kg/h of the feed rate of raw sugarcane bagasse. Results show that the height of the reactor column must be at least 30 cm for that are 1 mm in diameter and 108 cm for particles that are 2 mm in diameter.

Keywords

Torrefaction, Biomass, Sugarcane Bagasse, Moving Bed, Method of Lines

1. Introduction

Torrefaction is a thermochemical process in an inert (nitrogen gas is flowed into the reactor) or limited oxygen environment where biomass is slowly heated at a temperature between 200˚C and 300˚C [1] [2].

The torrefaction process occurs in two big steps. The first step consists of a drying in which biomass loses the majority of its moisture at 100˚C and, the second one in which hemicellulose, lignine and cellulose devolatization occur at 180˚C - 220˚C.

The importance of torrefaction lies in the fact that it leads to the production of fuel with high energy density, grindability and hydrophobicity.

Based on the reactor design, torrefaction technologies are subdivided into many categories such as moving bed, fixed bed and screw conveyor. In the moving bed, biomass is fed at the top of the reactor and it moves downward with respect to reactor wall. Furthermore, on the basis of the heat transfer mechanism, torrefaction technologies fall under two categories including direct and indirect heating. In the direct heating reactor, the heat is exchanged through direct contact between the heat carrier and the biomass. However, in the indirect heating reactor, heat is exchanged through a wall or through an electromagnetic radiation.

For a direct heating integrated torrefier, meaning that torrefaction reactor combines drying and torrefaction steps into one unit, as depicted in Figure 1, heat transfer analysis plays a critical role for successful initial sizing of the reactor. But the governing equations are very complex because they are nonlinear, coupled, etc. Therefore, they are very difficult to solve analytically. To alleviate this situation, numerical methods can be used (finite differences, finite elements, etc.). However, finite differences and finite elements methods have the disadvantage of being computationally more expensive than the method of lines. This is due to the fact that in the method of lines, independent variables are not discretized at the same time such with finite differences or finite element methods. This allows to reduce the time as well as the cost of computations [3].

The present work proposes the utilization of the method of lines to solve the equations describing the heat transfer between gas phase and solid phase within a direct heating integrated moving bed torrefier.

2. Methodology

The methodology adopted in this work and described below is a combination of three main parts such as the mathematical modeling, the collection of data and the numerical calculation. First, in the mathematical modeling part, a model allowing the estimation of a direct heating integrated moving bed torrefier was performed and used. Then, in the second part, the data collection was proposed due to the literature review. This allowed to describe the solid (sugarcane bagasse) and gas (nitrogen) phase physico-chemical properties. Finally, in the

Figure 1. Different temperature zones in a direct heating integrated torrefaction reactor [4].

numerical calculation part, the method of lines based on finite difference approach to solve a system of coupled partial differential equations obtained from the first part was performed.

3. Mathematical Modeling

3.1. Governing Equations

While using the eulerian-eulerian two-phase approach, the mathematical equations describing the heat exchange between the downward moving solid phase and the upward moving hot gas phase are based on the equations of energy conservation written for each phase. In this approach, different phases are treated as interpenetrating continuous fluids and the concept of phase volume fraction is need to be introduced.

Model assumptions

· The solid particles conduction effects are neglected as Biot number: $Bi=\frac{h{d}_p}{k_s}<0.1$

· The physico-chemical properties of solid and gas phase are constants;

· The solid particles are homogeneous spheres of diameter d_p;

· The gas phase flows upward at constant speed and the solid phase moves counter-currently to it at a constant speed too as depicted Figure 2;

· The heat exchange takes place only in an axial direction (1D problem).

Considering that $\alpha$ is the volume fraction of one phase, $\rho$ the bulk density (kg/m³), C_p the specific heat (J/kg˚C), T the temperature (˚C), $v$ the velocity (m/s) and k the heat conductivity coefficient (W/m˚C). The rate of the convection heat transfer between gas and solid particles is ${\varphi }_{sg}=h{A}_s\left(T_g-T_s\right)$, in which h is the convection heat transfer coefficient (W/m²˚C), A_s, the specific surface area (m²/m³) given by $6\frac{1-\epsilon }{d_p}$, $\epsilon$, the fractional voidage and d_p, the particles diameter.

Energy conservation equations

Figure 2. Simplified drawing the reactor.

For the case where the heat exchange by radiation and by friction is neglected, the energy conservation equations for the control volume dxdydz (Figure 3) and in the time dt can be written as follows:

$E_{in}-E_{out}=\text{d}E$ (1)

For solid phase, we have:

$E_{in}=h{A}_s\left(T_g-T_s\right)\text{d}x\text{d}y\text{d}z\text{d}t-k_s{\frac{\partial T_s}{\partial x}|}_x{\left(\text{d}y\text{d}z\right)}_{solid}\text{d}t$ (2)

$E_{out}=-k_s{\frac{\partial T_s}{\partial x}|}_{x+\text{d}x}{\left(\text{d}y\text{d}z\right)}_{solid}\text{d}t$ (3)

By substituting Equation (2) and Equation (3) into Equation (1), we have:

$\text{d}E=h{A}_s\left(T_g-T_s\right)\text{d}x\text{d}y\text{d}z\text{d}t+k_s\frac{{\partial }^2{T}_s}{\partial x^2}{\left(\text{d}x\text{d}y\text{d}z\right)}_{solid}\text{d}t$ (4)

Dividing by dt on both sides, we obtain:

$\frac{\text{d}E}{\text{d}t}=h{A}_s\left(T_g-T_s\right)\text{d}x\text{d}y\text{d}z+k_s\frac{{\partial }^2{T}_s}{\partial x^2}{\left(\text{d}x\text{d}y\text{d}z\right)}_{solid}$ (5)

On the other hand, the total amount of heat in ${\left(\text{d}x\text{d}y\text{d}z\right)}_{solid}$ is:

$\text{d}E={\rho }_s{\left(\text{d}x\text{d}y\text{d}z\right)}_s{c}_{ps}\text{d}{T}_s$ (6)

Hence the time rate of change of dE is:

$\frac{\text{d}E}{\text{d}t}={\rho }_s{\left(\text{d}x\text{d}y\text{d}z\right)}_s{c}_{ps}\frac{\text{d}{T}_s}{\text{d}t}$ (7)

In the Eulerian point of view:

$\frac{\text{d}{T}_s}{\text{d}t}=\frac{\partial T_s}{\partial t}+\frac{\partial T_s}{\partial x}\frac{\text{d}x}{\text{d}t}$ (8)

Consequently, by inserting Equation (7) and Equation (8) into Equation (5), we obtain:

${\rho }_s{\left(\text{d}x\text{d}y\text{d}z\right)}_s{c}_{ps}\left(\frac{\partial T_s}{\partial t}+\frac{\partial T_s}{\partial x}\frac{\text{d}x}{\text{d}t}\right)=h{A}_s\left(T_g-T_s\right)\text{d}x\text{d}y\text{d}z+k_s\frac{{\partial }^2{T}_s}{\partial x^2}{\left(\text{d}x\text{d}y\text{d}z\right)}_{solid}$ (9)

Figure 3. Control volume.

Since ${\left(\text{d}x\text{d}y\text{d}z\right)}_s={\alpha }_s\text{d}x\text{d}y\text{d}z$ and $\frac{\text{d}x}{\text{d}t}=v_s$, we obtain the partial differential equation below:

${\alpha }_s{\rho }_s{c}_{ps}\left(\frac{\partial T_s}{\partial t}+v_s\frac{\partial T_s}{\partial x}\right)=h{A}_s\left(T_g-T_s\right)+{\alpha }_s{k}_s\frac{{\partial }^2{T}_s}{\partial x^2}$ (10)

Considering the gas phase and applying the same procedure, we find the equation below:

${\alpha }_g{\rho }_g{c}_{pg}\left(\frac{\partial T_g}{\partial t}+v_g\frac{\partial T_g}{\partial x}\right)=-h{A}_s\left(T_g-T_s\right)+{\alpha }_g{k}_g\frac{{\partial }^2{T}_g}{\partial x^2}$ (11)

Fractional voidage

The volume fraction of a phase, $\alpha$, is defined as the ratio of the volume occupied by this phase to the total volume of the bed:

${\alpha }_s=\frac{V_s}{V_{bed}}$ (12)

Since $V_{bed}=V_s+V_{void}$, so the volume fraction of the solid phase can be calculated as follows:

${\alpha }_s=1-\frac{V_{void}}{V_{bed}}$ (13)

In the Equation (13), the ratio $\frac{V_{void}}{V_{bed}}$ is called the fractional voidage or void fraction ( $\epsilon$ ). Knowing that the gas phase occupes the volume $V_{void}$, volume fractions of the solid and gas phases can be written in terms of the fractional voidage as depicted in the Equation (14).

${\alpha }_s=1-\epsilon$ et ${\alpha }_g=\epsilon$ (14)

3.2. Simplified Equations

Considering all the assumptions mentioned above and using the fractional voidage instead of the volume fractions, the energy conservation equations become:

$\{\begin{array}{l}\left(1-\epsilon \right){\rho }_s{C}_{ps}\frac{\partial T_s}{\partial t}=\left(1-\epsilon \right)k_s\frac{{\partial }^2{T}_s}{\partial x^2}+h{A}_s\left(T_g-T_s\right)+\left(1-\epsilon \right){\rho }_s{C}_{ps}{v}_s\frac{\partial T_s}{\partial x}\\ \epsilon {\rho }_g{C}_{pg}\frac{\partial T_g}{\partial t}=\epsilon k_g\frac{{\partial }^2{T}_g}{\partial x^2}-h{A}_s\left(T_g-T_s\right)-\epsilon {\rho }_g{C}_{pg}{v}_g\frac{\partial T_g}{\partial x}\end{array}$ (15)

For ease of notation of the above system of partial differential equations, let us write:

$\left(1-\epsilon \right){\rho }_s{C}_{ps}=a_s;\left(1-\epsilon \right)k_s=b_s;$

$h{A}_s=c_s;\left(1-\epsilon \right){\rho }_s{C}_{ps}{v}_s=d_s;$

$\epsilon {\rho }_g{C}_{pg}=a_g;\epsilon k_g=b_g;$

$-h{A}_s=c_g;-\epsilon {\rho }_g{C}_{pg}{v}_g=d_g$

$\{\begin{array}{l}{a}_s\frac{\partial T_s}{\partial t}=b_s\frac{{\partial }^2{T}_s}{\partial x^2}+c_s\left(T_g-T_s\right)+d_s\frac{\partial T_s}{\partial x}\\ a_g\frac{\partial T_g}{\partial t}=b_g\frac{{\partial }^2{T}_g}{\partial x^2}+c_g\left(T_g-T_s\right)+d_g\frac{\partial T_g}{\partial x}\end{array}$ (16)

3.3. Initial Conditions

For the system of partial differential equations obtained above (Equation (16)) to be well-posed, we must include some auxiliary conditions to complete the statement of this mathematical model. The number of those auxiliary conditions depends on the highest order derivative in each independent variable (spatial and temporal). Since we have a system of two mixed hyperbolic (first order in x)-parabolic (second order in x) type partial differential equations containing the first order time derivative, we need to introduce two initial conditions and four boundary conditions as shown in Table 1.

3.4. Numerical Calculation

Method

The method of lines will be applied to compute the one-dimensional temperature profile within the torrefier. In this method, all spatial derivatives are approximated using finite differences, finite elements or finite volume method, and the semi-discrete (discrete in space but continuous in time) system of ordinary differential equations obtained is solved using the time integration method [5].

In this work, the first spatial derivative is approximated with the centered finite differences schema and the second spatial derivative with the three points centered finite differences schema, all with uniform grid spacing Δx.

$\{\begin{array}{l}{\frac{\partial T}{\partial x})}_i=\frac{T_{i+1}-T_{i-1}}{2\Delta x}\\ {\frac{{\partial }^{2}T}{\partial x^2})}_i=\frac{T_{i+1}-2T_i+T_{i-1}}{\Delta x^2}\end{array}$ (17)

As we define the uniform grid of N nodes, then we obtain a system of 2N linear coupled ordinary differential equations as depicted in the Equation (18) and where T is the vector containing nodal values of temperatures.

$\frac{\partial T}{\partial t}={\left(\begin{array}{c}A\end{array}\right)}_{2N\times 2N}T+b$ (18)

Table 1. Initial and boundary conditions.

$T={\left[\begin{array}{cccc}{T}_1& T_2& \cdots & T_{2N}\end{array}\right]}^{\text{T}}={\left[\begin{array}{cccccc}{T}_{s1}& \cdots & T_{sN}& T_{g1}& \cdots & T_{gN}\end{array}\right]}^{\text{T}}$

Then, the entries of the $2N\times 2N$ matrix A and the $2N\times 1$ vector b for all internal grid points ( $i=2,\cdots ,N-1;j=N+2,\cdots ,2N-1$ ) are:

$\{\begin{array}{l}A\left(i,i+1\right)=\left(\frac{b_s}{a_s\Delta x^2}+\frac{d_s}{2a_s\Delta x}\right)\\ A\left(i,i\right)=-\left(\frac{2b_s}{a_s\Delta x^2}+\frac{c_s}{a_s}\right)\\ A\left(i,i-1\right)=\left(\frac{b_s}{a_s\Delta x^2}-\frac{d_s}{2a_s\Delta x}\right)\\ A\left(i,i+N\right)=\frac{c_s}{a_s}\\ b_i=0\\ A\left(j,j+1\right)=\left(\frac{b_g}{a_g\Delta x^2}+\frac{d_g}{2a_g\Delta x}\right)\\ A\left(j,j\right)=\left(\frac{-2b_g}{a_g\Delta x^2}+\frac{c_g}{a_g}\right)\\ A\left(j,j-1\right)=\left(\frac{b_g}{a_g\Delta x^2}-\frac{d_g}{2a_g\Delta x}\right)\\ A\left(j,j-N\right)=-\frac{c_g}{a_g}\\ A\left(j,j-N\right)=-\frac{c_g}{a_g}\end{array}$ (19)

While using boundary conditions defined above we see that the remaining entries (for external grid points) must be written as follows:

$\{\begin{array}{l}A\left(1,1\right)=-\left(\frac{2b_s}{a_s\Delta x^2}+\frac{c_s}{a_s}\right)\\ A\left(1,2\right)=\frac{2b_s}{a_s\Delta x^2}\\ A\left(1,N+1\right)=\frac{c_s}{a_s}\\ b_1=0\\ A\left(N+1,N+1\right)=\left(\frac{-2b_g}{a_g\Delta x^2}+\frac{c_g}{a_g}\right)\\ A\left(N+1,1\right)=-\frac{c_g}{a_g}\\ A\left(N+1,N+2\right)=\left(\frac{b_g}{a_g\Delta x^2}+\frac{d_g}{2a_g\Delta x}\right)\\ b_{N+1}=\left(\frac{b_g}{a_g\Delta x^2}-\frac{d_g}{2a_g\Delta x}\right)T_{entree\text{'}}\end{array}$

$\{\begin{array}{l}A\left(N,N\right)=-\left(\frac{2b_s}{a_s\Delta x^2}+\frac{c_s}{a_s}\right)\\ A\left(N,N-1\right)=\frac{2b_s}{a_s\Delta x^2}\\ A\left(N,2N\right)=\frac{c_s}{a_s}\\ b_N=0\\ A\left(2N,2N-1\right)=\frac{2b_g}{a_g\Delta x^2}\\ A\left(2N,2N\right)=\left(\frac{-2b_g}{a_g\Delta x^2}+\frac{c_g}{a_g}\right)\\ A\left(2N,N\right)=-\frac{c_g}{a_g}\\ b_{2N}=0\end{array}$ (20)

Model parameters

The physico-chemical properties of both phases are reported in Table 2 and Table 3.

The convective heat transfer coefficient between spherical particles of sugarcane bagasse and gas is calculated using the following Ranz-Marschall correlation [7] because it suits for packed bed of uniform spheres [8]:

$h=\frac{k_g}{d_p}\left(2+1.8\times R{e}^{0.5}\times P{r}^{1/3}\right)$ (21)

In which $R{e}_p=\frac{{\rho }_g\times v_g\times d_p}{{\mu }_g}$ is the Reynolds number and $Pr=\frac{{\mu }_g\times C_{pg}}{k_g}$ is the Prandtl number.

In addition, since in a moving bed we have to avoid fluidization of the bed, the

Table 2. Physico-chemical properties of sugarcane bagasse [6].

Table 3. Physico-chemical properties of N₂ gas [7].

relative velocity between solid particles and gas (v_r) must remain less than minimum fluidization velocity (v_mf). That means the upper extrem value of the relative velocity is v_mf.

While taking into account the above descriptions, the relative velocity can be written as follow:

$v_r=v_g-v_s=v_{mf}$ (22)

For sugarcane bagasse of mass flow rate of ${\stackrel{\dot{}}{m}}_b$, continuously flowing downward the reactor of diameter D, the axial speed of the bed (solid particles) can be calculated as depicted in the Equation (23):

$v_s=\frac{4\times {\stackrel{\dot{}}{m}}_b}{\pi \times D^2\times {\rho }_s}$ (23)

Hence, by combining the two above Equation (22) and Equation (23), the gas phase velocity can be written as:

$v_g=v_{mf}+\frac{4\times {\stackrel{\dot{}}{m}}_b}{\pi \times D^2\times {\rho }_s}$ (24)

Considering the mathematical model developed in this work, the calculations will be made for a pilot scale cylindrical moving bed reactor, rated to process 5 kg/h of sugarcane bagasse and having a diameter of 30 cm.

Table 4 shows the minimum fluidization velocity of sugarcane bagasse as function of particle size.

Calculations

Before to solve the $2N\times 2N$ initial value problem Equation (18) obtained from the spatial discretization of the system of partial differential equations Equation (16), we first need to analyze its stiffness for the choice of an effective numerical method. To do that, the eigenvalues of the matrix A was calculated.

We can notice that from around 250 lines, the stiffness ratio (the ratio of the magnitudes of the real parts of the largest and smallest eigenvalues) depends heavily on the number of lines. Thus, as we increase the number of lines, the accuracy of the numerical solution increases [9], but the problem become stiffer meaning that explicit methods work extremely slowly [10] due to the need of a very small Δt to avoid numerical instability (Figure 4).

Table 5 shows a comparison of performance of implicit and explicit Runge-Kutta algorithms as function of numbers of lines.

Thus, we used the built-in Matlab solver ode15s based on implicit Runge-Kutta algorithm instead of ode45 based on explicit Runge-Kutta algorithm with a residence time of 10 minutes.

Table 4. Minimum fluidization velocity [6].

Figure 4. Stiffness ratio as a function of number of lines.

Table 5. Accuracy and computational cost (for d_p = 2 mm).

The temperature profiles along the axial direction are given in Figure 5 and Figure 6 for particles that are respectively 1 mm and 2 mm in diameter. One can see that the height of the reactor column must be at least 30 cm for that are 1 mm in diameter and 108 cm for particles that are 2 mm in diameter.

Furthermore, it is noticed that the effect of particles diameter on the minimum column height is not negligible. These temperature profiles indicate that the torrefaction of larger biomass particles require a reactor with a larger minimum column height. This is due to the fact that the heat transfer through the bed depends on biomass particles size (diameter). Indeed, large biomass particles have a resistance to the thermal propagation superior to that of small ones [11] [12].

Figure 7 depicts a qualitative diagram for biomass temperature as it passes through the torrefier column. The wet biomass enters the torrefier from the top at ambient temperature (point O), and as it gradually moves down, its temperature rises from ambient temperature. At point A, biomass reaches the drying temperature of 100˚C. Between A and B, the temperature remains constant while biomass loses its moisture. Then, between B and C, the temperature of biomass continues to rise to the desired temperature of torrefaction at which biomass begins to lose a fraction of its weigh in the form of volatiles. Finally, this temperature is maintained.

Figure 5. Temperature profile for d_p = 1 mm.

Figure 6. Temperature profile for d_p = 2 mm.

Figure 7. Qualitative diagram for biomass temprerature distribution along the column of a vertical torrefaction reactor.

Comparing Figure 5, Figure 6, and Figure 7, the shape of results presents good agreement with the one’s obtained qualitatively because of the similarity of trends.

Therefore, it can be concluded that the method of lines, proposed in this work, is useful for the estimation of the minimum column height of an integrated fixed or moving bed biomass torrefaction reactor.

4. Conclusions

The numerical simulation of the heat transfers between gas and solid phase within a direct heating integrated moving bed torrefaction reactor by the method of lines formed the subject matter of this work. The goal was to determine its minimum column height in order to complete the first part of the initial sizing of a pilot scale reactor. The effect of particle size on the reactor height was also evaluated. Emphasis has been given on the implementation of the method of lines because it seems to be a simple and versatile approach to provide numerical solution of coupled and uncoupled partial differential equations.

In order to validate our model, we first compared the obtained curves with the qualitative diagram as the experimental results are not yet available.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Tumuluru, J.S., Sokhansanj, S. and Baardman, R.D. (2011) A Review on Biomass Torrefaction Process and Product Properties for Energy Applications. Industrial Biotechnology, 7, 384-491. https://doi.org/10.1089/ind.2011.7.384
[2]	Rousset, P., Davreux, F., Macedo, L. and Perré, P. (2011) Characterization of the Torrefaction of Beech Wood Using NIRS: Combined Effects of Temperature and Duration. Biomass and Bioenergy, 35, 1219-1226. https://doi.org/10.1016/j.biombioe.2010.12.012
[3]	Sadiku, M.N.O. and Obiozor, C.N. (2000) A Simple Introduction to the Method of Lines. International Journal of Electrical Engineering Education, 37, 282-296. https://doi.org/10.7227/IJEEE.37.3.8
[4]	Basu, P. (2013) Biomass Gasification, Pyrolysis, and Torrefaction: Pratical Design and Theory. Elsevier Inc., USA. https://doi.org/10.1016/B978-0-12-396488-5.00007-1
[5]	Samir, H., William, E.S. and Graham, W.G. (2009) Method of Lines. Scholarpedia, 2, 2859. https://doi.org/10.4249/scholarpedia.2859
[6]	Hugo, T.J. (2010) Pyrolysis of Sugarcane Bagasse. Master of Sciences in Engineering.
[7]	Yaws, C.L. (1995) Thermophysical Properties of Chemicals and Hydrocarbons. 830-838.
[8]	Ranz, W.E. (1952) Friction and Transfer Coefficients for Single Particles and Packed Beds. Chemical Engineering Progress, 48, 247-253.
[9]	Esence, T., Bruch, A., Molina, S., Stutz, B. and Fourmigué, J.F. (2017) A Review on Experience Feedback and Numerical Modeling of Packed-Bed Thermal Energy Storage Systems. Solar Energy, 153, 628-654. https://doi.org/10.1016/j.solener.2017.03.032
[10]	Hairer, E. and Wanner, G. (1991) Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems. Springer-Verlag, Berlin, 5-8. https://doi.org/10.1007/978-3-662-09947-6
[11]	Mikhail, M.N. (1987) On the Validity and Stability of the Method of Lines for the Solution of Partial Differential Equations. Applied Mathematics and Computation, 22, 89-98. https://doi.org/10.1016/0096-3003(87)90038-5
[12]	Pamuk, S. and Erdem, A. (2007) The Method of Lines for the Numerical Solution of a Mathematical Model for Capillary Formation: The Role of Endothelial Cells in the Capillary. Applied Mathematics and Computation, 186, 831-835. https://doi.org/10.1016/j.amc.2006.08.024