On the Mathematical Modelling of a Moving-Bed Counter-Current Gasifier Fuelled with Wood-Pellets

The subject of this work is the mathematical modelling of a counter-current moving-bed gasifier fuelled by wood-pellets. Two versions of the model have been developed: the one-dimensional (1D) version-solving a set of Ordinary Differential Equations along the gasifier height-and the three-dimensional (3D) version where the balanced equations are solved using Computational Fluid Dynamics. Unique procedures have been developed to provide unconditionally stable solutions and remove difficulties occurring by using conventional numerical methods for modelling counter-current reactors.The procedures reduce the uncertainties introduced by other mathematical approaches, and they open up the possibility of straightforward application to more complex software, including commercial CFD packages. Previous models of Hobbs et al., Di Blasi and Mandl et al. used a correction factor to tune calculated temperatures to measured values. In this work, the factor is not required. Using the 1D model, the Mandl et al. 16.6 kW gasifier was scaled to 9.5 MW input; the 89% cold-gas efficiency, observed at 16.6 kW input, decreases only slightly to 84% at the 9.5 MW scale.


Introduction and Objectives
Combustion and gasification of biomass has become considerably important in the current energy scenario due to both resource protection and climate precaution (CO 2 neutrality). Biomass importance strongly varies from country to country depending on the financial incentives instituted. For example, co-firing of coals with biomass [1][2][3] became commonplace in the Nordic countries, in the UK and in the Netherlands, while it is hardly used in Germany. Denmark is certainly the leading country in producing electricity in full-scale boilers fired with biomass [4][5][6][7].
The conversion of biomass in a number of thermal and chemical processes has also being promoted [8,9]; however, wide large-scale commercialization has not yet occurred. Biomass gasification is typically carried out in fluidized-beds, entrained-flow reactors or in fixed-or moving-beds. Gomez-Barea and Leckner [10] reviewed the fluidized-bed technology, Higman [11] reviewed the entrained-flow gasifiers, while Chopra and Jain [12] reviewed fixed-bed technology. The most recent developments in Germany include bioliq entrained-flow technology [13,14] and TU Freiberg entrained-flow gasifier [15].

Counter-Current Gasification of Solid Materials
In moving-bed counter-current gasification, a gasifiable solid material is fed into the top of a reactor, as shown in Figure 1. Typical feedstocks are low-quality coals, wood, wood-waste, Refuse Derived Fuel (RDF) or RDF with CaO (Ecoloop) [16]. A gasification agent, usually air or an air-water vapour mixture, enters at the reactor bottom and flows through the bed upwards. The process is typically divided into four zones: drying, pyrolysis, gasification and combustion, as shown in Figure 1. The moist input material first undergoes drying at approximately 100°C. The material then flows through the pyrolysis zone and is pyrolysed at 200-500°C temperatures. After the volatile components are given off, the solid-bed contains char and ash. Part of the char then reacts with the gasification agent (H 2 O, CO 2 and H 2 ) in the gasification zone. The remaining char burns in the lower part of the reactor (combustion zone) and releases the energy necessary for the zones above. Ash is discharged at the bottom of the gasifier. There is no sharp separation between these four regions, as there is inevitably some overlap between them [17,18]. The continuously fed feedstock moves slowly downwards so that the bed-height remains approximately constant. Such a gasifier we regard as a moving-bed reactor as opposed to fixed-bed reactors realising batch processes.

Objectives
Over the last two decades, mathematical modelling of biomass gasification in fixedbeds has become commonplace. Although a detailed literature review is beyond the scope of this paper, the models can be simplistically classified into thermodynamic equilibrium modelling [21,22], models based on balance equations written separately for gas and solid-phases with over-simplified fluid-flow [17,18,23,24], CFD-modelling based on porous media approach or neural networks [25]. It is perhaps fair to say that the CFD-modelling approach, initially developed for grate-fired stockers and waste-incinerators [26,27], has nowadays become widely applied to biomass gasification in fixed-beds.
Most of the mathematical modelling work has been concerned with performance predictions of co-current gasifiers where a biomass and a gasifying agent are supplied at the top (down-draft) and both move downwards. The most recent developments concerning co-current steady-state gasification include the home-made models of Sharma [28] and Jaojaruek [29] and the CFD-models (based on the ANSYS-FLUENT package) of Janejreh and Al Shrah [30]. There are numerous numerical simulations of batch processes of biomass combustion [31] and gasification [32] where a gasification agent or combustion air is blown (typically from the bottom) through a biomass bed that remains stationary. Biomass gasification with high-temperature air [33] has also being modelled [34], confirming an increased gasification rate and a higher producer-gas yield. In the above cited works, either the Newton-Raphson method (home-made software) was used, or the CFD solver was used. Although, in the latter case, under-relaxation was required, generally, there were no difficulties in obtaining converged solutions.
The situation is different in the counter-current configuration shown in Figure 1, where inlet conditions for gas-phase balance equations are specified at the bottom, which is also an outlet for the solid-phase. Similarly, inlet conditions for the solid-phase are specified at the top, which is an outlet for the gas-phase. Thus, an integration along the z-coordinate (Figure 1), from bottom to top, requires guessed values of the solid-phase at z = 0 so as to match the solid-phase inlet conditions at z = L. If CFD-solvers were to be used, the same problem surfaces; however, in a different way-the inlet (numerical) cells for the air, entering the gasifier at the bottom, cannot be used for outlet conditions of the solid-phase. The existence of this difficulty has been recognized for a long time [18,23], and ad hoc adjustments were proposed [17,24]; none of them were general enough to be successful. Not surprisingly, while mathematical modelling of co-current combustion and gasification has progressed well over the last two decades, modelling of counter-current reactors has become a stumbling block. The objective of this paper is to close this gap and equip a onedimensional mathematical model of Mandl et al. [24] with a new numerical procedure that unconditionally provides stable and accurate simulations for counter-current gasification. In order to remain consistent with Mandl et al.'s [24] modelling work, whose experiments are used as validation, we use here kinetics (Section 3), which is almost identical to the one used by Mandl et al. [24]. We also demonstrate how to perform three-dimensional, CFDbased, simulations of counter-current gasification without any need for ad hoc corrections.

Experiments
The experimental investigations, which we use for validation, were carried out at TU Graz (Mandl [24]) using a laboratory-scale gasifier shown in Figure 2. (In Mandl et al. [24], the gasifier is described as a fixed-bed reactor. In this paper, we call it a moving-bed reactor, see Section 1) The aim of the experiment was to investigate the thermal-chemical conversion of wood-pellets. The pellets (quality standards DIN plus or ÖNORM M 7135), see Table 1, were continuously fed into the gasifier from above at a rate of 3.5 kg/h so as to keep a 0.42 m bed height. Air was supplied from below as a gasification agent. The reactor cylindrical tube (height: 0.6 m; diameter: 0.125 m) was insulated to minimise heat loses. During the experiments, the temperature profile over the height was measured (T1-T7 in Figure 2); the producer-gas composition (CO 2 , CO, H 2 O, H 2 , CH 4 , tar) was recorded at the reactor head only [24].

Mathematical Modelling
This section contains a detailed description of the mathematical model. The equations are presented in a coordinate-independent form as in our previous publication [35]. Onedimensional (1D) and three-dimensional (3D) versions of the model and their predictions are going to be discussed in Section 5. Throughout this paper, the subscript g stands for the gas-phase while the subscript s for the solid-phase. Sinks/sources in the massbalance equations are designated using the letter G (kg/(m 3 s)), while sinks/sources in energy-balance equations are marked using the letter S (J/(m 3 s)). The values of the kinetic parameters and reaction enthalpies are listed in Table 2, which is presented before the model equations to facilitate easy reading. Table 2.
Kinetic data and heats of reaction for drying, pyrolysis, heterogeneous and homogeneous reactions.

Conservation Equations for the Gas-Phase
The energy balance equation for the gas-phase takes the form where is the moving-bed porosity; ρ g stands for the gas density; u g is the gas velocity; h Gas is the physical enthalpy of the gas; k g is the thermal conductivity; S konv is the source term representing the convective heat transfer between the gas and the solid-phase (see Equation (51)); S i is the source term due to chemical and physical processes (see Reactions (30)-(32), (40) and (43)); S f g represents the source term associated with the en-ergy exchange between the gas-and solid-phases, see Equation (50); and S V,g,W stands for the heat loses through the reactor walls. The momentum balance equation takes the form where p stands for the static pressure, g is the gravity, τ stands for a Newtonian fluid tensor for a laminar flow, and S solid is the pressure drop in the moving-bed. Equation (2) is accompanied by the continuity equation for the gas-phase: For each of the components (N 2 , O 2 , CO 2 , CO, H 2 O, H 2 , CH 4 and tar), a mass conservation equation is solved, which reads: where G g i stands for six source terms as follows where Y i is the mass fraction of the ith component, and D i,e f f stands for the effective diffusion coefficient.

Conservation Equations for the Solid-Phase
The energy balance equation for the solid-phase takes the form where S konv stands for the convective heat transfer between the two phases, and S drying , S pyrolysis , S gasi f ication and S combustion represent the sources due to drying, pyrolysis, gasification and combustion. S V,s,W stands for the heat loss through the reactor walls, and S f s is the source term due to the phase change (see Equation (49)). In the 1D model, the radiative heat transfer within the solid-bed is accounted for through an effective thermal conductivity, following the work of Goldman et al. [42], while in the 3D-model, it is accounted for through a Discrete Ordinates (DO) radiation model [43]. The continuity equation for the solid-phase reads: and is accompanied by the mass conservation equation for the i-th component: where G s i is an appropriate sink or source term (G drying , G pyrolysis , G gasi f ication , G combustion ). Ash is considered as inert, and its mass fraction is calculated as: is modelled according to the following Arrhenius expression [17,24]: R·Ts (11) with where R is the universal gas constant.

Rate of Pyrolysis
The dried pellets are pyrolysed into components using the following reaction dry feedstock k (13) −→  [17,24]. The pyrolysis rate is described by the following relationship [17,24] with After the bound water and the volatile components are given off, the feedstock consists of ash and char CH 0.2526 O 0.0237 only.

Heterogeneous Char Combustion Reaction
The char is oxidised following the reaction with the temperature-dependent χ = CO/CO 2 ratio calculated as [44]: with k c = 2500 · e (− −6420 Ts ) (18) and the mass transfer coefficient: The kinetic rate constant of Equation (16) is calculated as (see Table 2) The overall rate of reaction (16) takes into account both the mass transfer and kinetics so that where A p is the activated specific surface in square meter per cubic meter, and C O 2 stands for the concentration of oxygen in the gas-phase. The rate of consumption of oxygen in reaction (16) is: Heterogeneous Char Gasification Reactions Char gasification is described by the following scheme: with the overall reaction rates (24), (25) with the mass transfer coefficient: Overall Char Consumption Rate The following term is the sink/source term, which describes the overall char consumption rate due to both combustion and gasification:

Gas-Phase Combustion Reactions
The following three combustion reactions are considered with their rates calculated as follows [17] where M i is the Molecular Mass of (CO,H 2 , and CH 4 ) and C i are the concentrations of CO, O 2 , H 2 and CH 4 . The appropriate sinks in the gas-phase balances are given by Water Gas-Shift Reaction The rate of the homogeneous water gas-shift reaction is calculated as [17,41]: where the equilibrium constant K eq is [45]: The source is then Tar Cracking The tar, produced during pyrolysis, cracks in the gas-phase as follows [17,24]: with the rate where ρ tar stands for the tar density. Thus,

Sinks and Sources in Energy Balances
Sources and sinks occuring in the Energy Balances are calculated as follows: where c p,w,g stands for the mean specific heat of the gaseous water, and T 0 is the reference temperature of 298 K, where c p,Volatiles is the mean specific heat of the volatiles. The source/sink terms for the heterogeneous char combustion and gasification are while the source term (S f g ), which appears in Equation (1), is calculated as follows

Heat and Mass Transfer, Bed Structure and Bed Thermal Properties
The heat transfer between the gas-phase and the solid-phase is accounted for through the source-terms appearing in the energy balance equation of the solid-phase (Equation (6)) and the gas-phase (Equation (1)) where the convective heat transfer coefficient is calculated as In the above equation, a correction factor ξ occurs [17,24,46], and its value varies between 0.02 and 1 [46], see Section 5.1.
The heat losses occurring through the reactor walls (see (Equations (1) and (6)) are calculated as: where T w stands for the wall temperature of 300 K with i = g, s and In the present work, the assumption was made that the feedstock in the reactor neither disintegrates, accumulates nor rubs off, and thus, the original shape of the pellets is retained. Therefore, the moving-bed porosity ( = 0.5) remains constant for the entire reactor. Mandl et al. [24] demonstrated that the wood-pellets remained stable in all zones. Consequently, the specific surface area of the particle A p is calculated as [24,39] Thus, the initial equivalent diameter (d p ) of the fed material (standardized wood-pellets) is d p = 0.010259 m for pellets diameter d cylinder = 0.006 m and length l cylinder = 0.02 m; u s stands for the solid velocity, while u s,0 represents its initial value.
The thermal conductivity and the gas dynamic viscosity are calculated as [47] The effective thermal conductivity of the solid-bed is calculated as and consists of the thermal conductivities of moisture and char [17,18,23,24], where k s0 is the thermal conductivity of the moist feedstock [47]: while k char is the thermal conductivity of the char [42]: The influence of thermal radiation is accounted for through the two terms k radiation,g and k radiation,s [42,48] k radiation,g = 4 · σ · 0.05 · T 3 g (61) The mean specific heat-capacity of the gas-phase components is calculated using a 4th order polynomial. The temperature interval is divided into low (300-1000 K) and high (1000-5000 K) temperature regions. The mean specific heat capacities of the solids are listed in Table 3.

Differences in the Mathematical Modelling
It is perhaps fair to say that the basis for mathematical modelling of fixed-and movingbed reactors was formulated by the work of Brigham Young University [18,23], where emphasis was placed on coal gasification. Di Blasi [17] extended and adapted the mathematical framework to handle gasification of 5 mm wood particles in a 50-centimetre long and 10-centimetre diameter fixed-bed. Mandl's et al. [24] work is based on Di Blasi's one-dimensional equations [17], while our model presented above, although developed in one-dimensional and three-dimensional versions, is essentially based on the work of Di Blasi [17] and Mandl et al. [24]. Thus, all the three models considered here-namely Napoli's [17], Graz's [24] and our model-have the same mathematical framework. The models include very similar mathematical descriptions of drying, pyrolysis, char gasification and char oxidation. For example, for drying and pyrolysis, a shrinking-density sub-model is used, while for char combustion and gasification, a shrinking-core sub-model is used. The bed-structure is described by an identical set of relationships with the assumption of constant bed-porosity. Reaction schemes for the gas-phase are identical with the same tar cracking scheme, oxidation reactions of CO, H 2 and CH 4 , including the water gas-shift reaction. The essential differences are in: (a) the correction factor, see Equation (52), and (b) the numerical solver.

The Kinetics
In the experimental work of Di Blasi [17], beech wood was gasified, while in the work of Mandl et al. [24], wood-pellets (see Section 2) were used. Neither for the beach wood nor for the wood-pellets fuel-characterization experiments/measurements were carried out to determine the pyrolysis-products yields and the pyrolysis rate. Neither rates of char heterogeneous reactions were measured nor rates of gas-phase reactions, including tar cracking, were experimentally determined. Instead, in the mathematical modelling [17,24], the kinetic data were taken from the literature (see Table 1 in Di Blasi [17] and Table 2 in Mandl et al. [24]). For example, the stoichiometric coefficients that appear in Equation (13) were determined by Di Blasi [17] from the work of Roberts and Clough [50], where a different wood was examined. The coefficients were retained in Mandl et al. [24] and in this work. A detailed analysis of all the kinetic data used is beyond the scope of this paper. Table 4 and Figure 3 are produced to show the kinetic rates for wood pyrolysis and tar cracking only. In our work, the pre-exponential factor and the activation temperature for both pyrolysis and cracking are taken from Glaister [51], as cited in Grønli [36] (see Table 4.5 in Grønli [36]). Di Blasi [17] derived the pyrolysis rate from Roberts and Clough [50] and the cracking rate from Liden et al. [52]. Mandl et al. [24] derived them from Grønli [36] and Liden et al. [52], respectively (see the footnote to Table 4). Figure 3 shows that the rates are similar. Table 4. Kinetic rates expressions used in the modelling; the rates in 1/s.

Sub-Model
This Work (see Table 2  Most of the data concerning heterogeneous kinetics, listed in Table 2 and used in Di Blasi [17]and Mandl et al. [24], as well as in their work, were generated in the 1980s and 1990s of the last century using thermogravimetry. In this period, thermogravimetry was regarded as the (perfect) tool for heterogeneous kinetics, and the mass and heat transfer limitations of the technique were not fully realized. Furthermore, difficulties in obtaining reliable values of both the pre-exponential factor and the activation temperature (energy) were unknown [53]. In the light of the current knowledge [53][54][55][56], the accuracy of the data is questionable. In other words, modern thermogravimetry would likely produce different figures [57]. Currently, one observes a proliferation of publications, as exemplified by Refs. [58,59], containing more comprehensive kinetic data on biomass devolatilization, combustion and gasification.

Numerical Solution
For stationary one-dimensional (1D) problems, the equations presented in Section 3 simplify to a set of ordinary differential equations (ODE) describing temperature, the mass fraction of N 2 , O 2 , CO 2 , CO, H 2 O, H 2 and CH 4 and tar for the gas-phase and temperature, wood, wood-char and ash for the solid-phase. The spatial z-coordinate ( Figure 1) is then the independent variable. The gas-phase equations are coupled with the solid-phase equations by the source terms described in Section 3; additional coupling is provided by the relationships describing the bed-structure. Since the gaseous gasification agent enters the gasifier at the bottom, while the wood-pellets are supplied at the top, that gasifier operates in a counter-current mode. Thus, the inlet conditions for the gasification agent are known at the bottom, while the inlet conditions for the wood-pellets are specified at the top. When ODEs are required to satisfy boundary (inlet) conditions at more than one value of the independent variable (z-coordinate), the resulting mathematical problem is called a two-point boundary value problem. Hobbs et al. [23] named such a situation as "a split boundary value problem" and developed an elaborative iterative procedure to satisfy the temperature boundary conditions at the top and bottom of the gasifier. Di Blasi [17] used a first-order Euler method to solve the ODEs; no information was given on how the boundary conditions were handled. Mandl et al. [24] used a two-step iterative method; in the first step, guessed values were used for the unknown boundary conditions at the top where the integration begun and proceeded downwards using an explicit Runge-Kutta [60] procedure. In the second step, the unknown (guessed) boundary conditions were varied using the secant method [60] to match the boundary conditions at the bottom. The method, often called shooting, is inherently unstable, and good initial guesses are the secret to convergence.
In this paper, we present a novel method of handling the 1D two-point boundary value problem of counter-current gasification. The method is simple and eliminates mathematical difficulties associated with conventional solvers. Figure 4 shows an organigram of the method. In the first initialization step, the gas-phase temperature and mass fractions of N 2 , O 2 , CO 2 , CO, H 2 O, H 2 and CH 4 are given values (for all values of z-coordinate) determined by the inlet (bottom) conditions; the solid-phase temperature and mass fractions of watervapour, volatiles, char and ash are initialized with the top inlet values. After initialization, the gas-phase ODEs are solved upwards, beginning at the bottom, using the multi-step solver (MATLAB ® , ode15s) for stiff equations. During the integration of the gas-phase equations, the solid-phase variables remain frozen. Now all the source terms are evaluated (step 4 in Figure 4); no linearization is required. Then the solid-phase ODEs are solved downwards, from the top to the bottom, using the same solver (MATLAB ® , ode15s). During this integration, the gas-phase variables remain frozen. Then, again, the source-terms are evaluated (step 6 in Figure 4). The procedure is terminated when the following convergence criteria are satisfied: for the mass fraction of all gas-phase and solid-phase components for gas-phase and solid-phase temperatures where indices k and k + 1 indicate two consecutive iterations (an iteration includes a solution of the gas-phase ODEs upwards and subsequent solution of the solid-phase ODEs downwards). The above criteria have to be satisfied at all numerical divisions (grid points) in the z-coordinate. After a converged solution is found, both the overall mass and energy balances are checked. In all calculations presented here, the mass and energy balances are closed with an accuracy of 10 −2 -10 −4 of the inlet values.
For the 3D version of the model described in Section 3, the CFD solver ANSYS ® FLUENT (version 17.1) for porous-beds is used. The problem of how to impose appropriate boundary conditions, at the bottom and at the top, still remains. Figure 5 shows the organigram of the solution procedure; the essence is the use of two instances of the CFD solver run in parallel. The first instance solves the gas-phase equations only, after which the sources are evaluated and passed on to the second instance, which solves the solid-phase equations. The communication between the two instances, facilitated every 25 iterations, is provided by the User-Defined Functions. The same 3D numerical grid is used for both instances. The procedure is terminated when the above convergence criteria are satisfied.
We wish to finish this Section with a remark. The above-described 1D and 3D procedures for computing of the counter-current reactor are used for steady-state calculations. We are sure that the same procedures can be used in time-dependent computing.

Results
Predictions of both the 1D and 3D versions of the mathematical model described in Section 3 are considered here. The predictions correspond to the experimental run of Mandl et al. [24] so that the boundary and inlet conditions listed in Table 5 are identical  to Table 4 of Mandl et al. [24]. The gasifier's diameter and height are 0.125 m and 0.42 m, respectively. Wood-pellets of a 0.0103 m equivalent diameter are gasified in a countercurrent, steady-state operation. With a 3.5 kg/h wood-pellets supply rate and 17 MJ/kg Lower Calorific Value of the wood, the gasifier operates at 16.5 kW thermal input. The reactor was ignited introducing an ignition energy source in the middle of the reactor; the source ceased after 10 iterations. Since the energy balances for the solid-and gas-phases are second-order differential equations, two additional boundary conditions must be specified for the upper and lower part of the counter-current gasifier, which are as follows [17,39]: Extensive studies have been performed to obtain grid-independent numerical solutions. In 1D computing, 100, 1000, 2000 and 3000 divisions (along 0.42 m gasifier height) were tested, while in 3D computing, 63,360, 214,900 and 1,235,600 nodes were used. The 1D gridindependent results were obtained using 2000 nodes, while in the 3D grid-independent results, 214,900 nodes were used. Details of the sensitivity studies can be found in Schwabauer [20], see also Sosnowski et al. [61].

Comparison with the Measured Data
Figure 6 (Top) shows the measured and predicted temperatures of both the gas-phase and the solid-phase. The 1D and 3D predictions, both obtained with ξ = 1, match each other very well, demonstrating the correctness of the method used for gas-phase/solid-phase coupling. In the 3D simulations, the heat-loss through the reactor walls is the only element that may bring 3D effects. Since the heat-loss is small, amounting to 4% of the fuel input (=16.58 kW), all dependent variables vary marginally with the gasifier radius. Thus, it is not too surprising that 1D and 3D predictions match well.
Differences between the predicted gas-phase and solid-phase temperatures are typically not larger than 50 K. The predicted temperatures match the measurements well, and the differences are typically not larger than 75 K; slightly larger values than measured are predicted at elevations above 0.1 m. The four zones, namely drying, pyrolysis, gasification and combustion, are also shown in Figure 6 (Top). Arbitrarily chosen criteria for their boundaries are as follows. The drying zone, beginning at z = 0.42 m, ends up when the moisture content of pellets is smaller than 10 −4 , while the pyrolysis zone ends when the mass fraction of the volatiles remaining in pellets is smaller than 10 −5 . The thickness of the oxidation zone corresponds to the z-coordinate at which the oxygen mass fraction in the gas-phase is smaller than 10 −4 . Indeed, the gasification zone is the largest, as one would expect, with gasification reactions proceeding with appreciable rates in the region where temperatures are larger than around 1000 K. The model predicts a rather thin pyrolysis zone as a consequence of the very fast pyrolysis rate. Figure 6 (Bottom) shows TU Graz's predictions (see Figure 5 in Mandl et al. [24]), which generally are of the same quality. It is important to realize that TU Graz's predictions were obtained using different numerics and with a 0.5 value of the correction factor. Figure 7 shows gas composition predictions using both 1D and 3D model versions. Again, both versions provide very similar values with the largest deviations occurring in CO predictions. This is attributed to slight differences in the gas-phase temperature predictions; carbon monoxide concentration-being temperature-dependent through the water gas-shift reaction-is then slightly altered. Both versions provide good predictions of the gas composition at the gasifier outlet, as shown in Table 6. Cold-gas efficiency, listed in Table 6, is calculated as [19,62] η cold gas e f f iciency =ṁ producergas · LCV producergaṡ m wood−pellets · LCV wood−pellets (64) Table 6 also contains predictions of Mandl et al. [24] where under-predicted CO, H 2 and CH 4 concentrations at the gasifier outlet (see Figure 6 in Mandl et al. [24]) are to be noticed. To remedy this shortcoming, the yield of the pyrolysis reaction was modified [24], and the improved predictions are shown in the fifth column of Table 6. No such adjustments were made in our simulations.

The Correction Factor
The correction factor ξ, appearing in Equation (52), was already introduced into the modelling in the 1980s and 1990s [18], so it is around forty years old. As shown in Equation (52), the factor governs the energy transfer between the gas-and solid-phases; for ξ = 1, the rate is maximum, while for ξ = 0, there is no energy transfer between the phases. Thus, with ξ approaching 1, the differences between gas-phase and solid-phase temperatures are minimized. In all the previous works [17,18,24], the factor was used for tuning the temperature predictions to the measured values, and, due to imperfections of the numerical solvers (Section 4.2), obtaining stable and converged solutions was possible for some (selected) values of ξ only. On the contrary, the solvers presented in Section 4.2 provide unconditionally stable and converged solutions for any ξ values in the range zero to one. Figure 8 demonstrates importance of ξ, showing both gas-phase and solid-phase temperatures for ξ= 1, ξ = 0.5 (Mandl et al. [24]) and ξ = 0.2 (Di Blasi [17]). The onset of the pyrolysis zone is a strong function of ξ, as shown in Figure 8. Table 7 shows the dependence of ξ on the predicted composition and LCV of the producer-gas.  Table 7. The influence of the correction factor ξ on the predicted composition and LCV of the producer-gas as well as on cold-gas efficiency.

Producer Gas
Measurements

Scaling
The small-scale (16.6 kW) gasifier of Mandl et al. [24] produced a low-BTU gas of 6.7 kJ/kg (LCV) with 89% cold gas efficiency (Table 6), which was indeed a good performance. Thus, the question of scaling the design to industrial-size is relevant.
Combustion system scaling is generally a complex procedure, particularly when flames are to be scaled [63]. However, for moving-bed reactors, the procedure is straight forward; the key parameters are the residence time of the solid-phase and the air-to-fuel ratio. According to the 1D predictions, the solid-phase remained in the reactor for 0.88 h. Table 8 shows a series of gasifiers for which the solid-phase residence time is kept at 0.88 h, and the air-to-fuel ratio is 1.45. The smallest thermal input of 16.6 kW (standard case in Table 8) corresponds to Mandl et al. [24] experiments, while the largest input is 9.5 MW with four scales between. All the gasifiers remain in geometrical proportions since the height-to-diameter ratio remains at 3.6 so that the velocity of gasification air remains the same. At each scale, wood-pellets of the same composition and size are used. It is seen in Table 8 that the calorific value of the producer-gas tends to decrease with increasing input, while the cold-gas efficiency decreases from 87% for 16.6 kW input to 84% for the 9.5 MW scale [20]. Figure 9 shows a substantial increase of both solid-phase and gas-phase temperatures with the reactor size. This is attributed to the reduction of heat loss with the increased size. The heat loss amounts to around 4% of 16.6 kW thermal input, while at 9.5 MW, 0.5% heat loss is applicable. In other words, the effect of the surface/volume ratio of the reactor is clearly visible. Since the exit temperatures increase with thermal input, the cold-gas efficiency decreases accordingly. The predictions show that only ash leaves the reactor bottom so that the full conversion of the wood-pellets is observed at all scales. Table 9 and Figure 10 show the influence of scale on the gas composition. At all scales, CO, CO 2 and H 2 O are formed in the combustion zone; CO 2 , H 2 O and H 2 are then consumed in the gasification reactions. The 6.5% CO 2 concentration in the producer-gas, observed at the 16.6 kW scale, drops to approximately 4.5% at 9.5 MW, while CO increases from 21.8% to 25.3%. The tar concentration is reduced from 4.4% to approximately 3.7% with increasing the scale, which is desirable. As described above, the calorific value of the producer-gas decreases slightly with scale [20]. In the above-described gasifier scaling, we have computed both the gas-phase and solid-phase temperatures as well as gas-phase composition along with the gasifier height. Table 8 lists the predicted gas-phase temperature and the calorific value of the producergas at the gasifier outlet for the six scales considered. A question is which of these two gasifier exit parameters (gas temperature or calorific value) has been predicted with a larger confidence? In light of the relatively close agreement (shown in Figure 6) between the measured and predicted temperatures, one may expect that the local heat-release rates have been correctly predicted for all the scales considered. Since the heat losses constitute a rather small fraction (see above) of the thermal input, we place a larger confidence on the temperature predictions than on the calorific value predictions.

Conclusions
The essence of this paper is the steady-state mathematical model of a moving-bed reactor for gasification of wood-pellets in air. Two versions of the model have been developed to simulate counter-current operation: the one-dimensional (1D) version-solving a set of Ordinary Differential Equations along the gasifier height-and the three-dimensional (3D) version where the balanced equations are solved using Computational Fluid Dynamics used separately for the gas-phase and solid-phase equations. Unique procedures were developed to provide unconditionally stable solutions and remove difficulties occurring from conventional numerical methods [17,18,23,24] in modelling counter-current reactors. The procedures reduce the uncertainties introduced by other mathematical approaches, and they open up the possibility of a straightforward application to more complex software, including commercial CFD packages.
Mandl et al. [24] describe both experiments and 1D mathematical modelling of the small-scale (16.6 kW) gasifier of wood-pellets. Their mathematical model contains a correction factor, appearing also in previous publications [17,18], which was used to tune the temperature predictions to the measured values. It has been demonstrated in this paper that the factor is not required.
The 16.6 kW gasifier of Mandl et al. [24], producing a low-BTU gas of LCV = 6.6 MJ/kg, operates at 89% cold-gas efficiency, which is indeed a good performance. Using the 1D model, the performance of a series of scaled versions of the reactor was examined, with the largest size of 9.5 MW. Calculations indicate a slight decrease of the producer-gas LCV from 6.6 MJ/kg (at 16.6 kW) to 6.3 MJ/kg (at 9.5 MW). The cold-gas efficiency decreases slightly from 89% to 84%, confirming the attractiveness of the design.
Mandl et al. [24] and Di Blasi [17] models, as well as the model presented in this paper, contain a large number of kinetic data, all taken from the literature. Now, after the uncertainties in the numerical solver have been removed, there is a need to examine the sensitivity of the models to the kinetics used. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.