It is perhaps fair to say that the basis for mathematical modelling of fixed- and moving-bed 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.
4.1. 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.
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.
4.2. 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
2O, 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
2O, H
2and 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 water-vapour, 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 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 – 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.