On the Thermal Stability of a Counter-Current Fixed-Bed Gasiﬁer

: In recent years, gasiﬁcation gained attention again, both as an industrial application and as a research topic. This trend has led to the necessity to understand the process and optimize reactors for various materials and conﬁgurations. In this article, the thermal structure of a counter-current reactor is investigated to demonstrate that constraints on the temperature mainly determine the oxidation and the pyrolysis region. A non-dimensional set of equations is written and numerically solved using the method of lines (MOL) with spatial discretization based on a spectral algorithm. The results show that four thermal structures can be identiﬁed, two of which are the most common ones found in reactors of practical applications. Two stationary operation positions have been determined, one in the upper and one in the lower part of the reactor. Existence and stability conditions have been discussed based on non-dimensional parameters. The knowledge derived from this analysis was applied to two conﬁgurations, one typical of a biomass gasiﬁer and one proposed for waste gasiﬁcation.


Introduction
Gasification is a well-known technology that has been used in industrial plants for nearly 200 years.In the beginning, coal and peat were used to produce light gases for energy or heat production.In fact, conversion of heterogeneous materials into light gaseous fuels such as methane, carbon monoxide, and hydrogen made it possible to increase the efficiency because of the higher adiabatic temperature.Later on, other technologies replaced gasification in the energy sector, and only recently, gasification gained attention again, both as an industrial application and as a research topic.In particular, gasification of less calorific fuels (peat, biomass, plastics, and various kind of wastes) has become of considerable importance in the current energy scenario due to both resources protection and climate precaution, specifically CO 2 neutrality (see for the example [1][2][3][4][5][6][7]).Nevertheless, the recycling aspects of the thermal treatment of complex materials for the modern circular economy [8][9][10] are no less important.
The change of feedstock materials requires, in addition to the knowledge of chemical reactions, their products, and their dependence on thermodynamic variables, a thorough and accurate knowledge of the functioning of the reactors used in the various processes.This article proposes an original approach to the mathematical analysis of a particular type of reactor, namely the counter-current (updraft) reactor.
Several types of reactors are currently available for biomass gasification: countercurrent fixed bed, co-current fixed bed, fluidized bed, entrained flow, plasma, and free radical reactors.Studies conducted by Gomez-Barea and Leckner [11] reviewed the fluidizedbed technology, Higman [12] studied the entrained-flow gasifiers, while Chopra et al. [13,14] studied the fixed-bed technology.In the present day, there has been a resurgence of interest in this technology, as evidenced by the numerous publications available in the literature.In entrained flow gasifiers, a dry pulverized solid, an atomized liquid fuel, or a fuel slurry is gasified with an enriched oxygen agent in co-current flow.High temperatures and pressures cause a high throughput; however, thermal efficiency is somewhat low, as the gas must be cooled before it can be cleaned with existing technology [15,16].The fuel is fluidized in oxygen and steam or air in fluidized-bed reactors.Since temperatures are usually relatively low, the fuel must be highly reactive.Fuel throughput is higher than in the fixed bed technology but not as high as in the entrained flow gasifier [11,[17][18][19][20][21].In co-current fixed-bed gasifiers (also known as downdraft type), the gasification agent flows in the same direction as the fuel (downwards, hence the name).Heat needs to be added to the upper part of the reactor, either by oxidizing a small amount of fuel or from external heat sources.The produced gas leaves the gasifier at a high temperature, but since most of this heat can be often transferred to the gasification agent injected at the top of the bed, the resulting energy efficiency is similar to the one achievable with the counter-current reactor type.Since the gas must pass through the hot bed, the tar level is much lower than in the counter-current reactors [13,[22][23][24][25][26][27][28].In counter-current fixed-bed gasifiers (also known as updraft type), the gasification agent and the solid fuel flow in a counter-current configuration.The throughput for this type of reactor is relatively low, while the thermal efficiency is high as the temperatures in the gas exit are relatively low.However, as a consequence, tar and methane production is significant in typical operation conditions, so the product gas must be extensively cleaned before use [29][30][31][32][33][34][35][36].
Despite their advantages and decades of economic and scientific efforts, the fixed-bed gasification process, especially on industrial scales, has yet to establish itself in the market.The low efficiency of the co-current process, the high tar content in the gas in the countercurrent process, and the resulting high costs for gas and wastewater treatment often lead to the limits of economic viability.For example, the most recent developments in Germany are based on the entrained-flow technology instead and include the Bioliq reactor [5,15], and TU Freiberg entrained-flow gasifier [37].In combination with these intense research activities, efforts have also been made, in the recent past, to develop processes for clean thermal treatment of the most challenging fuels, such as biomass, substitute fuels (RDF), and waste with high proportions of chlorine, sulfur, but also heavy metals.In particular, among the various projects started recently, the Ecoloop process [38][39][40] was an attempt at a new innovative counter-flow gasifier that converts various carbonaceous materials into purified syngas.The Ecoloop process was based on a counter-current moving-bed gasifier with a high proportion of calcined mineral rock.This material, almost chemically inert, serves as a support structure and allows the primary cracking of tars, and, simultaneously, is responsible for binding halogens, sulfur, and heavy metals.Due to this innovative process, and in contrast to other standard thermal conversion processes, no harmful organic substances such as dioxin or furan are produced.The movement of the bed allows for a continuous reactor operation and recirculation of the solid material.Energy is not released into the environment but kept inside the process using the recirculation of the calcined rocks.A project has been supported by the German ministry [41] to develop the technology from a laboratory scale to a complete industrial application.
During the execution of the project, the use of mathematical models for the analysis of the behavior of the reactor under different conditions became indispensable but all the models presented in the literature lacked the necessary generality to be applied to this new type of gasifier.In fact, in the past and present literature, there is a multitude of onedimensional models for the simulation of fixed bed reactors, with a large part limited to the gasification of coals.Accurate and detailed investigations were carried out by the Advanced Combustion Engineering Research Center of the Brigham Young University in Provo (UT, USA) [42].Despite many simplifications, the model was successfully validated using experiments in a pilot plant.Similar models have been successively developed primarily based on this work (see [43][44][45][46][47]).
On the contrary to the large body of literature available for coals, models for the fixed-bed gasification of non-coals are very sparse.Di Blasi described an unsteady model for the gasification of wood in a counter-current fixed-bed configuration [48].In addition, these models were successfully validated using pilot-plant scale measurements.Similar models, using wood and jute sticks as fuel, can be found, for example, in [49,50].Costa et al. [22] presented a phenomenological model flexible to the feedstock composition.The model was developed with the aim of investigating the advantages and limits in the utilization of feedstock different from woodchips, such as hydro-char derived from the hydrothermal carbonization of green waste, or a mix of olive pomace and sawdust.Recently Schwabauer et al. [30] developed a 1D and 3D model for the gasifier configuration described by Mandl et al. [51] in which, together with experiments, a model is also presented.All the aforementioned models are comprehensive of the main chemical and physical processes, including moisture evaporation/condensation, finite-rate kinetics of wood devolatilization and tar degradation, heterogeneous gasification (steam, carbon dioxide, and hydrogen) and combustion of char, combustion of volatile species, finite-rate gas-phase water-gas shift, extraparticle mass transfer resistances, heat and mass transfer across the bed resulting from macroscopic (convection) and molecular (diffusion and conduction) exchanges, solidand gas-phase heat transfer with the reactor walls, and radiative heat transfer through the porous bed.In addition to the investigation presented by Di Blasi et al., Schwabauer et al. compared a simple 1D model with a more complex and detailed 3D model, demonstrating the validity of the simpler 1D procedure.Both authors presented also results of calculations with parameters varied in a small interval.However, the number of input parameters in both investigations is considerably high and the changed quantities were chosen either for their importance in the process or for their theoretical uncertainty.Therefore, neither of the two authors mentioned above had the interest of analyzing the results from a more purely mathematical and theoretical point of view.
In particular, the investigations mentioned above, both experimental and theoretical, have never investigated the thermal structure of the reactor and its thermal stability as a function of the relevant thermal parameters of the system.The results presented in the literature and highlighted above, are unable to justify why the various types of reactors work as they work, which kind of thermal characteristics can be achieved, and the relevant parameters influencing this structure.
Consequently, new industrial reactors are developed using an analogy principle, assuming that the thermal structure is similar to known, well-functioning reactors, despite the differences in many process parameters.The analogy principle has been the base for the Ecoloop reactor development, whose operational capabilities were based on the thermal structure of Limekiln, a reactor used to calcinate limestone (calcium carbonate) to produce calcium oxide.During the project mentioned above [41], this assumption has been rigorously investigated, and the limitation of the analogy principle has been demonstrated.
As a contribution in closing the gap, this article presents a new analysis of the oxidation front position and its stability in a counter-current continuous moving-bed configuration.The configuration described by Mandl et al. [51] (denoted as MA configuration) and the design suggested in [38][39][40] (denoted as EL configuration) are taken as a reference, but the analysis is far more general.The configurations above are considered as typical design for a standard gasifier (the one presented by Mandl et al.) and for an innovative design (the one shown in [38][39][40]).To circumvent the difficulties of interpreting the reactor phenomenology based on many parameters, three hypotheses have been postulated.Qualitatively and, with a good approximation, also qualitatively correct understanding of the reactor phenomenology can be obtained when: 1.
The equations are formulated in non-dimensional form; 2.
They are simplified to such an extent that they can be easily solved, in same cases also analytically; 3.
Non-dimensional parameters which arise in the non-dimensional equations are responsible alone for the thermal structure of the system.
The decision to use equations in non-dimensional form is based on several reasons.The normalization procedure results in the formation of a set of a few parameters in terms of the many thermodynamic quantities available in the system.These parameters and not the individual thermodynamic quantities are decisive for the thermal structure of the reactor.Some of these parameters acquire an important physical meaning which makes it possible to measure the relative importance of individual processes and, eventually, give a quantitative measure of the chance to discard some of them.
As a fundamental consequence of the three hypothesis formulated above, the thermal structure of counter-current reactors is primarily due to static factors, such as thermal couplings, mass flow rations, and fuel composition, and secondarily, it is influenced by the dynamic details of the chemical reactions in the process.Based on these assumptions, the analysis is started by solving the simplified non-dimensional equations using a flamesheet approach.In a successive step, a refined solution is numerically achieved using the method of line with a spectral differentiation scheme.In this step, several simplifications made before are relaxed.Existence and stability conditions are discussed based on nondimensional parameters.The knowledge derived from this analysis is applied to the two configurations mentioned above, the MA configuration typical of a biomass gasifier and the EL configuration proposed for waste gasification.

The Equation of the Model
The gas and the solid phase are modeled as two different sub-systems, and their description is based on equations balancing the usual thermodynamic quantities:

•
Total mass balance (for the gas and the solid system); • Enthalpy balance written as an equation for the temperature (for both systems); • Mass balances for each component in each phase.
The equations are formulated for a 1D domain in conservative form for an unsteady process (see, for example [30,48]).The continuity equations for the gaseous and the solid phase are, respectively: where ρ is the density, ε is the porosity of the bed, U is the velocity, and Ṁa is the mass transferred per unit volume from the solid to the gaseous phase during the process a.
The z-axis is aligned with the flow of the gas; therefore, the gas velocity is positive, while the velocity of the solid is negative.
The energy equations for the gaseous phase and the solid phase are: where Qg and Qs are the sources due to reactions in the gas and solid phase, respectively, h g , h s are specific enthalpies, k is the effective conductivity, S c (ε) is the specific contact surface, and h(ε) is the convective heat transfer coefficient.The specific physical enthalpy of the material transferred from the solid phase to the gaseous phase is denoted h r,a .
The effective conductivity coefficient k s incorporates heat transferred by conduction through the solid phase and the heat transferred by radiation.In fact, in packed bed reactors, the energy transfer due to the radiation is active mostly between neighboring particles.If the diameter of the particles is much smaller than the reactor dimensions (as in both processes analyzed in this paper), the radiative heat transfer becomes nearly a local process.The linearized formula for the heat transfer due to radiation can be written as [52]: According to this relation, the source term due to radiation can be approximated using the gradient of the solid temperature and k rad can be incorporated into an effective conductivity coefficient.
The species equations for the gas and the solid phase are: where f s,i and f g,i are the mass fractions while D s,i and D g,i the diffusion coefficients.
The mass transfer rate of component i from the solid phase to the gaseous phase is denoted with Ṁa,i .Evidently ∑ i Ṁi = Ṁa .The homogeneous reactions are indicated with the terms ṁg,i and ṁs,i and, since homogeneous reactions will not be further considered, these terms will be dropped from the equations in the rest of the article.
The following thermochemical processes are responsible for the phase coupling: 1.
Heat exchange among the phases due to convection and radiation.In Equations ( 3) and ( 4), the thermal coupling due to conduction is explicitly written Drying and pyrolysis of the fuel.

3.
Heterogeneous reactions of the fixed char with the gaseous phase (heterogeneous oxidation and gasification reactions); 4.
Homogeneous reactions in the gas phase (homogeneous oxidation and water gas shift reaction).
The momentum equation for the gas phases reduces to the equation for the pressure: The function F(p, T, ε, Re) on the right side describes the local pressure loss in the system and can be calculated by a correlation similar to the Ergun equation [53].
The equation system must be completed with boundary conditions (BCs) for both phases, whose inlet conditions are known at z = 0 for the gas phase and z = L for the solid phase.The pressure is controlled at z = L.For the gas phase, they read: for the solid phase: and for the pressure: If diffusion is incorporated into the model, the equations are of second order and extra boundary conditions must be taken into account, neglecting diffusion compared with the convection.The assumed zero-gradient BCs are not correct for the energy equation, where a loss must be considered, usually by radiation.The inclusion into the model of this loss function would have required the introduction of additional parameters.Tests were also conducted considering a more correct form of BC, and the results, which are not reported here, were not significantly altered.
If the steady state solution is of interest, a simplified form of the previous equations must be solved with ∂ ∂τ = 0.The equations can then be expressed using the mass flow rates per unit area of the reactor cross section (Φ g and Φ s ) defined as: The heat capacities can also be introduced here, to be used later in the equations:

Assumptions and Simplifications
Several assumptions and simplifications are introduced, which also help the simplification of the mathematical problem and the reduction in the dimension of the parameter space to be analyzed: 1.
Ideal gas with constant specific heat.The specific enthalpy can then be calculated as follows: Ideal solid with constant specific heat.From this assumption, it follows that: Constant porosity ε.
Negligible mass diffusivity for both phases (D (s,g),i = 0).The last two assumptions imply that the problem is governed by advection and by the thermal conductivity of the solid only.In most industrial applications; however, the thermal conductivity of the solid can also be neglected (convective dominated problem), and in Section 3, the conductivity in the gas phase is neglected to simplify the analytical solution.
The assumption of constant specific heats causes an overestimation of the temperatures of both phases, and it is retained for simplicity, even if not strictly needed by the numerical method used to solve the equations.
Further assumptions are used to simplify the expression of the chemical reaction rates between the phases.

6.
Heterogeneous oxidation and pyrolysis are the only two processes considered here.7.
Pyrolysis is a single reaction in which the fuel is converted into gas and char.The rate is chemically controlled and is written using the classical Arrhenius rate: in which Y s,F is the mass fraction of the fuel in the solid.8.
The pyrolysis energy q p is negative, and the associated source term is present in the gas-energy equation only.This assumption is justified when the enthalpy of the gas alone determines the pyrolysis process.The gas phase has a higher temperature than the solid, and heat is transferred from the gas into the particle.A thin film on the particle surface is heated, and chemical bonds will be promptly disrupted before the heat can be transferred into the particle's internal region.Only the gas takes part energetically in the pyrolysis process.

9.
In the modeling of the oxidation rate, both kinetics and diffusion of oxygen in the bulk are considered.The reaction order for the oxygen is assumed equal to one, and the mathematical expression of the reaction rate takes the classical form: where R and D are the kinetic and diffusion reaction velocity per unit area, respectively.The term S c represents the area of contact between the phases per unit volume, and Y s,C is the mass fraction of the char in the solid phase.The kinetics velocity R is expressed in an Arrhenius form, while the diffusion velocity is defined as a power of the temperature T. 10.The heat of combustion q c is positive, and the associated source term is split between the gas and the solid using a factor f , expressing the amount of heat absorbed directly by the solid phase.

The Equations in Non-Dimensional Form
The system of Equations ( 3)-( 7) can be rewritten in non-dimensional form using the following procedure:

•
The length L of the reactor is used to introduce the variable • The non-dimensional temperature is defined as: where T 0 is the inlet temperature of the gas (usually room temperature) and T * a normalizing temperature taken T * =1000 • C. The ignition temperature is not used since its value changes depending on the reactor input conditions, and it is not a property of the fuel alone.As it has been chosen here, a fixed temperature allows for a quick comparison among the temperature profiles for different configurations.

•
The length of the reactor and the initial velocity of the solid are used to normalize the time: • The pressure is normalized using the atmospheric pressure: • The densities of both phases (ερ g and (1 − ε)ρ s ) are normalized to ρ g and ρ s using the initial density of the respective phase, namely at z = 0 for the gas and at z = L for the solid.

•
The fluxes Φ are also normalized using the initial values of the respective streams as: • The mass sources are normalized using the mass flow rate of the solid: • The heat of reactions (q p and q c ) are both normalized using the temperature difference ∆T * = T * − T 0 and the specific heat of the solid phase: Since the properties of the solid are used to normalize the sources and the heat of reactions, (see Equation (24)) both also present in the gas phase equations, new parameters must be introduced which describe the ratio of several properties in the two phases: express the ratio of the specific heats of the phases, express the ratio of the mass flows, express the ratio of the heat capacities.The parameter g introduced in Equation ( 27), can be related to the air stoichiometric ratio, commonly used in industrial applications, which is defined as the ratio between the air flow rate and the minimum (stoichiometric) air flow rate: Considering that Φ 0 g,O2,min = χ O2 Φ1 s,C,min , using the stoichiometry for the oxygen χ O2 after few steps, it is derived that: Using all the previously defined quantities, the non-dimensional equations can be expressed in the following form where τ is the non-dimensional time and x is the nondimensional form of the original variable z: where h r,a (T s ) = c p,r,a • (T s − T 0 ) is the enthalpy of the reaction products, c r,a = c p,r,a c p,s is the ratio of the products and the solid phase specific heats, and χ i,a are stoichiometric coefficients.
The following important variables naturally arise from the previous procedure: The variable ζ in Equation ( 37) appears in all gas phase equations since the normalization of the time has been performed using the velocity U s of the solid (see Equation ( 21)).Further, the variable g (defined in Equation ( 27)) appears in all the gas phase equations.The presence of this variable in the equations for the gaseous phase can be explained by considering that the source terms have also been normalized using the initial mass flow rate of the solid phase (see Equation ( 24)).The thermal coupling of the two systems is described by the parameter α.An increase in this value causes a decrease in temperature differences between the solid and the gas.This parameter determines the amount of heat recirculated inside the reactor, and, as a consequence, higher values of α cause an increase in the peak temperature.From Equation (38), it is evident that α scales linearly with the length of the reactors; therefore, industrial reactors such as the one used in the EL process tend to have a more tight thermal coupling (if other parameters remain similar).
The parameter δ, the inverse of the thermal Peclet number, describes the importance of heat transfer by diffusion and radiation compared to heat transferred by convection.If the parameter's value is smaller than one, then the reactor is advection dominated.An increase in the value of δ causes more heat to be transferred away from the reaction region and, therefore, a decrease in the solid system's peak temperature.Equation (39) shows that the parameter δ scales inversely with the length of the reactor.It is then expected that conduction and radiation play a lesser role in large industrial reactors than in small experimental reactors such as the one used in the MA experiments.
The equation for the pressure can be written finally in the following way: where The adiabatic flame temperature, used in combustion, is the temperature reached by a reacting mixture under adiabatic conditions.In terms of the non-dimensional quantities introduced so far, θ ad can be expressed by: where Y 1 s,C is the initial mass fraction of the carbon in the fuel.It is worth noticing that the adiabatic flame temperature is not the highest temperature that a counter-current system can achieve since heat loss from the reaction zone is recirculated back into it, as shown already in [54].
In Table 1, values for the non-dimensional parameters are reported in the case of the Ecoloop process (labeled as EL) and in the case of the experiments carried out by Mandl (labeled as MA).The values of the parameters that describe the pyrolysis and the oxidation kinetics are derived from the analysis performed in [30].In the case of the EL process, thermogravimetry data presented in [41] are used.In Table 1, other significant parameters are presented that are discussed later in this paper (see Section 3).
Table 1.Values of the non-dimensional variables used for the simulation of the MA [30] and EL process [41].From the data presented in the table, it is essential to note that in the case of the EL process, the parameter δ is smaller than unity, and, as already pointed out in the report [41], the convective processes largely determine the heat transfer in this reactor.Convection plays a significant role also in the MA process, but since the parameter δ is larger by a factor of 100, effects due to conduction and radiation in the solid phase are observable.A second significant difference between the two processes under consideration is presented by the values of the parameters ṁ, for both oxidation and pyrolysis.Since the EL reactor is fed with fuels and inert material, ṁ has smaller values than in the case of the MA experiments.

Oxidation Pyrolysis
The previous equations can be rearranged to express the global mass balance and the global energy balance.After summing up Equations ( 31) and (32), and considering steady state, the following equation is obtained: In the same way, after summing up Equations ( 33) and (34), in steady state: After integration of the previous equations over the entire x-domain, and taking care of the boundary conditions, the following relation is obtained: where: The relation (44) and (45) are used later for deriving an analytical solution.Equation (45) represents the well-known linear relation between the solid's exit temperature and the gas's exit temperature.When the gas leaves the reactor with a low temperature, the solid possesses a high temperature.Correspondingly, the solid leaves at a low temperature when the gas leaves the reactor at a high temperature.The exact linear behavior is just a consequence of the assumptions of constant thermal specific heats.At the same time, the inverse behavior of the exit temperatures also remains in a more general description of the reactor (see [30]).

The Mathematical Solution
The mathematical problem under investigation requires the solution of a system of partial differential equations with particular characteristics that make the problem difficult to solve by standard numerical methods (recognized in [30,51]).Specifically, the system is non-linear, stiff, and advection dominated.The last characteristic should be considered with particular caution because the strong coupling between the phases moving in opposite directions does not allow a straightforward application of the upwind scheme for the advection term.Perturbations that the upwind scheme handles on a single phase propagate backward on the second phase due to strong coupling and, therefore, can give rise to numerical instabilities.
Several discretization schemes have been tested, but the results are not detailed in this article.The method used for the calculations, which proved to be the most accurate, the most stable, and the fastest, is based on the method of lines [55,56].The spatial derivatives are discretized using finite differences based on spectral differentiation matrix [57,58].According to the theory, the domain x = [0, 1] is discretized using a Chebyshev grid.The resulting set of coupled ordinary differential equations (in the time variable) are solved using a stiff differential equations solver with a variable-order method available in the MATLAB suite codes [59,60].
If the steady state solution is of interest, additional difficulties arise, specifically, the fact that boundary conditions are specified at different endpoints depending on the phase to be calculated.Several approaches have been proposed.

•
The shooting method in which the differential equations are written as a system of first-order ordinary differential equations.Boundary conditions at one interval endpoint must be guessed or solved to meet the conditions at the other endpoint.The advantage of this method is that the solution of the ODE system can be easily calculated using highly accurate methods (such as Runge-Kutta methods), where the accuracy can be controlled by grid refining.

•
The method proposed by Schwabauer in [30] in which the gas phase and the solid phase are solved separately and iteratively until convergence is achieved.

•
The solution of a time-dependent problem (or pseudo-time dependent) for a long enough time interval until steady-state is reached.
Three kinds of difficulties arise when using the shooting method in a parametric study such as this one.The first disadvantage lies in the exponential form of the solution.Minor deviations from the correct unknown boundary conditions quickly produce nonphysical values of the solution, and without particular constrains the problem cannot be solved.The second disadvantage lies in the number of boundaries that must be sought.In simple applications such as the one used in Section 4.1, the number of unknowns at the boundary x = 0 is small, (only two unknowns boundary conditions are used in Section 4.1), the solid exit temperature and the solid exit mass flow.In more complex applications, the number of unknowns becomes larger and larger depending on the detail of the description of the solid system.Finally, when the reactor is operated in sub-stoichiometric conditions, the amount of volatile and fixed carbon (mass fraction w, for example) is equal to zero at the solid outlet.The unknown values cannot be adjusted from that boundary, so a procedure to choose the position and also an additional procedure for choosing the point to shoot at, must be introduced.Because of these drawbacks, the shooting method has been used only for simple modeling as in Section 3.

Solution Using the Flame Sheet Approximation
The steady-state solutions of the system of Equations ( 31)- (36) are derived using the flame sheet approximation in which the total reaction rate ṁ is expressed using two Dirac's deltas of the ignition temperatures: The parameters m o and m p represent the stoichiometric masses of the solid transferred from the solid phase to the gaseous phase during combustion and pyrolysis, respectively.According to the previous equation, the oxidation process is activated when the char reaches the ignition temperature T ig while the pyrolysis process is started when the gas reaches the pyrolysis temperature T p .The choice of this particular temperatures is in accordance with the discussion made in the description of the pyrolysis process in Section 2.2.
To further clarify this point, it is crucial here to mention that the ignition temperate in a counter-current reactor is different than the ignition temperature used in flame sheet models for reactors in co-flow configuration.The reason for the difference lies in the product Y s,C • p O2 present in the reaction rate.Since the distribution of oxygen and solid fuels inside the oxidation zone differ in co-current and counter-current flow, the ignition temperature used here is higher than the standard ignition temperature.More details are given in Section 5.4.
In the flame sheet approximation, the reaction front and the pyrolysis zone (also a sheet in this approximation) are positioned at x = x * and x = x p , respectively.Discontinuities for all the thermodynamic variables, specifically for the char and the oxygen mass fractions and the temperature profiles, are present at those two positions.The position of both reaction sheets must be determined as a function of the ignition and pyrolysis temperatures.
The temperature and the mass flux in the region x < x * (below the oxidation region) are denoted as θ − and φ − , respectively, the temperature and the flux in the region x * < x < x p (above the oxidation region) are indicated as θ + and φ + , respectively, and the temperature and the fluxes in the region x > x p (above the pyrolysis region) are denoted as θ ++ and φ ++ , respectively.The discontinuity in the thermodynamic variables due to oxidation is positioned where the solid temperature equals the ignition temperature: while the discontinuity due to pyrolysis is determined when: For all x different than both reaction zone positions, the system of Equations ( 31)-(36) simplify drastically.The equations for the fluxes and those for the components are (ordinary derivatives are used instead of the partial derivatives since the system is supposed in steady state): These equations imply that each component's fluxes and mass fraction remain constant throughout the reactor, with jumps on the sheet positions only.
The equations for the temperatures are slightly more complex: From the previous equations, the following properties can be inferred: 1.
The solutions of the equations for the total mass, Equations ( 50), are straightforward.

2.
Since the details of the flame front is not resolved in this approximation, the equations for the temperatures are decoupled from the Equations ( 51) for the fluxes and the concentrations.

3.
The equations for the temperatures can be analytically solved (see later).

4.
The position x * and x p are determined by the temperature profile as a function of the parameters of the problem.
The equation system must be closed with jump conditions at x = x * and x = x p which, for the fluxes at the oxidation front, are: and similar relations on the pyrolysis sheet: If δ = 0 and c r = 1 the jump conditions for the solid temperature are (here the equations for the oxidation front only): where: If δ = 0 and also c r = 1, the temperatures are continuous, and jumps are present in the first derivatives only: The jump conditions for the gas temperature are derived from the global energy balance (43): where q a is the non-dimensional enthalpy of reaction for the oxidation and for pyrolysis, respectively.The equations are too cumbersome to be reported and they do not add clarity to the text.

The Temperature Profile
Before analyzing the relations for the position of both reaction zones x * and x p , a short discussion of typical temperature profiles is of importance for further understanding.Even if, as already mentioned before, the profiles can be analytically calculated, the general formula is not given here because the analytically derived relations are rather involved and do not add clarity to the results presented in this paper.Therefore, to simplify the mathematical description and the final formulas, the analysis was first conducted for the case without pyrolysis and without conduction in the solid phase.
Equations ( 52) and ( 53) describe the temperature profiles when no reactions are present.These equations are identical to the ones valid in a counter-current heat exchanger, and the results are similar to those of the theory of those devices.Excluding the particular case in which w = φ g φ s , the temperature profiles are described by exponential functions.If the difference in temperature is defined as ∆ = θ s − θ g , the Equations ( 52) and ( 53) can be rearranged into the following two independent equations: where the parameter k in Equation ( 66) is defined as: The oxidation and the pyrolysis positions divide the reactor into three different regions with three different k i because the fluxes φ i are different in each region.When pyrolysis is not present, as considered in this section, two regions are only present and two k i can be defined, denoted as k 1 , and k 2 in the range x < x * and x > x * , respectively.
Using Equation ( 66) and the boundary values θ 0 g ≡ θ g (x = 0) and θ 1 s ≡ θ g (x = 1), the temperature profiles can be analytically calculated using the function F(s; k) = e k•s : where, using: and the values of w 1 and w 2 in Equations ( 68) and (69) are given by: Inserting them into Equation (67), it is possible to write the values of k 1 and k 2 in the following form: The temperature profiles for the whole x-range are calculated with the help of the two jump conditions given by Equations ( 58) and (59) that determine a relation between θ 0 g and θ 1  s as a function of the θ ig and x * .The concavity of the temperature profiles in each region of the reactor is determined by the factor k and the temperature difference ∆.In particular (see also Equation ( 66)), if k < 0, the temperature difference between the phases decreases exponentially, while if k > 0, the temperature difference increases exponentially as x increases.The concavity of the temperature profiles are determined according to the following cases: Four typical profiles can be identified, and typical temperature profiles are shown in Figure 1.The profiles correspond to configurations in which pyrolysis is inactive, and the flame sheet is positioned at x = 0.4.According to the previous explanation, in both profiles on the right of Figure 1 (the profile TC2 and the profile TC4), the concavity of the region x > x * is positive, while on the profiles on the left, the concavity in the same x-region is negative.Consequently, in the profiles TC2 and TC4, the temperature of the solid before the runaway (the solid phase move from right to left) is primarily low, growing exponentially near the ignition front only.This effect is more pronounced when the absolute value of parameter k 2 is large.
Both profiles on the left of Figure 1, present an opposite behavior.In particular, for x > x * the concavity of the temperature profile is negative, and the solid temperature quickly increases around x = 1.If the absolute value of k 2 is large (compared to one), an asymptotic temperature is quickly reached.In this case, the solid phase before the runaway point experiences a relatively high constant temperature plateau.The next section shows that the runaway can also be achieved in the upper part of the reactor.In such a case, the temperature profiles are similar but spatially mirrored from the ones shown in Figure 1.The solid phase can experience a high-temperature plateau for w < w 1 (because k > 0) after the runaway while its temperature is quickly decreed w > w 2 (because k < 0).

Conditions for the Existence of Each Thermal Condition
After determining the characteristics of the temperature profile in the previous section, the existing conditions of each profile are determined in this section.The inequality k 1 > 0 is satisfied when: while the inequality k 2 > 0 is satisfied when: In terms of w * 1 and w * 2 , the thermal state TC1 is characteristic for w < w * 1 and w < w * 2 at the same time, and that can be achieved with a small value of w.The thermal state TC2 requires w > w * 1 and w > w * 2 at the same time, and that can be achieved with a large value for w.The case of TC3 requires w * 2 < w < w * 1 , with the consequence that w * 2 < w * 1 .Considering the definition of w * 1 and w * 2 given by the Equations ( 73), these inequalities hold only if: This condition, together with the definition of the thermal state TC3 (w < w * 2 ), requires the specific heat ratio c < 1.Similarly, the thermal state TC4 can be obtained only if c > 1.The profiles for the thermal states TC1 and TC2 can be achieved with every value of c.
It is worth noticing, that the thermal coupling α does not play any role in the determination under which of the four thermal states the reactor operates.
The conditions derived for each thermal state can be recast in terms of stoichiometry number λ and composition of the fuel Y 1 s,C using the relation ( 79) and (80).The relation and the relation k 2 > 0 implies: The results of Equations ( 79) and (80) are presented in Figure 2, whose curves are derived for complete oxidation of pure carbon to carbon monoxide with standard air.The black and the thin magenta lines correspond to line of constant k 1 and k 2 , respectively.with the values in the following ranges k 1 ∈ [±0.5, ±1.0, ±2.0, ±4.0, 8.0] and k 2 ∈ [−0.95, −0.9, −0.8, −0.6, −0.4,0.5, 1.0, 1.05, 1.1, 1.2, 1.5, 2, 4, 8].The thick lines present the borders for the various thermal conditions.The left plot is representative of a c < 1, and the right is representative of c > 1.It is evident from the plots that the two thermal conditions TC3 and TC4 can be achieved in a narrow range of parameters, while in the most common practice, thermal condition TC1 is achieved.Condition TC2 is most likely representative when the fuel's fixed char is highly diluted with inert material.The decrease in the air excess ratio λ or the specific heat capacity ratio c both favor the achievement of the thermal condition TC2.

Relaxation of Various Assumptions
Before concluding the discussion of the temperature profiles, a few comments must be made, relaxing some assumptions used previously.Firstly, the calculations in this section have been performed using the flame sheet approximation without conduction in the solid phase (the parameter δ has been taken equal zero).Already Figure 1 presents, together with the profiles with δ = 0, profiles when a thin oxidation zone is present and with thermal conduction accounted in the solid phase.The calculations have been performed using δ = 10 −4 and parameters for the oxidation source term chosen in such a way that (see also Section 4.1) the reaction front thickness equals ∆x * = 10 −2 .As a consequence of these improvements, the temperature profiles are continuous functions.The energy released by oxidation is partially removed from the reaction region, and the gas and the solid temperatures outside the oxidation region increase, which is also observed in [54] for stationary adiabatic premixed flames in porous inert media.However, if δ is increased to larger values, as is shown in Section 4.2, energy is definitely removed from the reaction zone, and the temperatures decrease.The similarity of the profiles calculated with and without approximations is evident, demonstrating the correctness of the conclusions obtained using the reaction sheet approximation.
The second aspect worth mentioning concerns the possibility of composition dependent specific heats.If the specific heats c p are constant, the conditions w > 1 and w < 1 are determined by the initial conditions only and are not altered when crossing the reaction zones.In some reactors, the specific heat changes drastically crossing the reaction zones, mixing the conditions mentioned above.For example, this physical change is observed in Limekilns where the solid phase consists of a mixture of limestone and char.The char possesses a high specific heat, while the lime at the end of the shaft is small.It is common to have w > 1 in the upper part of the reactor and w < 1 after the oxidation zone.
The last aspect being commented on concerns the effect of the thermal coupling α.In Figure 3, the temperature profile previously presented as a TC1 example is shown for different values of the thermal coupling, namely α = [2.0,4.0, 6.0, 8.0] for a fixed position of the ignition (x * = 0, 2 in the figure).The position of the flame remains unchanged only if the ignition temperature T ig = T + s (x = x * ) increases as the thermal coupling increase.The more relevant case in which the ignition temperature is fixed is analyzed later.According to Equation (73), a stronger thermal coupling increases the values of both k i , shown in the figure with the increased curvature of the temperature profiles for both the gas and the solid phase.The temperatures in the reactor before the runaway (x > x * ) increase, but, since the inlet solid temperature is fixed, the influence of α is more pronounced near the oxidation zone and less pronounced at the gas exit (x = 1).The gas temperatures in the region x < x * are almost unaffected by α because the solid temperature is reduced as the thermal coupling increases.This last property is not as general as the first one, depending on which thermal condition is achieved in the reactor and the values of the parameters k 1 .More commonly, the temperatures increase also in this region as α increases, as shown in Figure 3

The Position of the Runaway
In the previous sections, the positions of the oxidation fronts have been chosen by appropriate ignition temperatures (around T ig = 0.6 in the case of Figure 1).In real systems, the position of the ignition x * is not a free parameter, but it is linked, using the jump conditions to the ignition temperature T ig .In this section, the position of the oxidation front is presented as a function of the ignition temperature.The shooting method is used to solve the equations and to capture all the solutions to the problem.The stability of those solutions will be investigated in Section 5.1.
The diagrams in Figure 4 show the runaway position x * as a function of the ignition temperature T ig for the parameters listed in Table 2. Together with the parameter α ∈ [10.0, 6.0, 3.0, 1.0], the parameter c, has also been varied with steps of ∆c = 0.1 for TC1, TC2 and TC4, and ∆c = 0.05 for TC3.Several general characteristics are evidential in the results: • Only materials with a finite range of ignition temperatures can reach a stable position in the reactor.

•
As a consequence of the assumption that the ignition temperature is equal to a specific solid temperature, an oxidation front is present even for small T ig , and it is located at a position near x = 1, the upper part of the reactor.• If a stable oxidation front is stabilized in the upper part of the reactor, an increase of T ig causes a decrease of x * , and the reaction front moves deeper into the reactor.

•
The ignition temperatures that cause the reaction front to be at the bottom part of the reactor x * ∼ x = 0, are higher than the ones which stabilize the front in the upper part.

•
With an increase of α, the maximum T ig clearly increases, and this means that, for materials with a fixed ignition temperature, the oxidation front shifts toward the two ends of the reactor as thermal couplings is increased.

•
In most of the conditions, two possible positions exist for a given T ig .Only for small thermal couplings α one position is present in the upper part of the reactor.
It is important to notice the differences in ignition temperatures allowed by the system and how the runaway is a function of the two chosen parameters α and c.While for each thermal condition (TC), an increase in the thermal coupling causes an increase in ignition temperature, an increase in the specific heat of the solid phase c causes a different trend for each thermal condition.In fact, an increase in the parameter c causes the stabilization of materials with a higher ignition temperature.In Thermal Condition 2 only, an increase of the c value cause the stabilization of fuels with a lower ignition temperature.
The form of the curves is also different for each thermal condition.In Thermal Conditions 1 and 3, the curves are smooth, and stabilization in the central region of the reactor is achieved at somewhat similar T ig .Conversely, in the conditions TC2 and TC4, the runaway in the reactor center are achieved by materials with ignition temperatures in a larger range.This different behavior can be of practical importance, considering that, as pointed out in Section 3.1, a runaway in the reactor's central region ensures a better process performance.
The qualitative features presented above can be completed by introducing other characteristic parameters.To introduce them, the same data as in Figure 4 are plotted in Figure 5 where the non-dimensional heating value q is presented as a function of the runaway position.The diagram (a) on the left presents the data for several chosen values of α = [0.2,1.0, 2, 0, 5.0, 10.0] and at a constant T ig = 0.7.In contrast, the diagram on the right presents the same information but for different ignition temperatures (T ig = [0.7, 1.0, 1.2, 1.5]) at a fixed value of α = 2.The following conclusions can be drawn from the figure : 1.
A critical parameter q cr (T ig ) exists for which there is no stabilization when q < q cr .This condition describes the situation in which the exothermic reactions and the thermal coupling cannot heat up the solid to the ignition temperature.2.
If q > q cr the number of solutions depend on a particular value α * : (a) If α < α * , only one runaway position is possible for the allowed range of ignition temperatures.For values of q slightly higher than q cr , the oxidation front stabilizes in the vicinity of x = 0. (b) If α > α * , the system respond to a variation of q in a more complex way.It exists a value q * for which, when q cr < q < q * , two positions of the front are allowed, while if q > q * , only one position is allowed at the reactor top.When q is slightly higher than q cr , the front stabilizes in a position near the center of the reactor.The curves (b) on the right in Figure 5, show that an increase in the parameter q cause a shift of the curves towards higher ignition temperatures without changing the qualitative behavior of the T ig − x relation.This observation means that oxidation stabilizes at the same position materials with higher q if the ignition temperature is higher, while if the T ig remains the same, the stabilization front moves towards the border of the reactor.
The thin lines in the same figure (b) in the right are obtained for c r = 1.5, a value characteristic of materials for which the specific heat of the products is higher than the specific heat of the material itself (see the definition c r = c p,r c p,s ).This condition can be achieved if the solid materials are based on a mixture of char with high c p -values and inert with low c p -values.This is the case for the EL configuration [41].The figure clearly shows that materials with higher c r have lower ignition temperatures at x = 0.For these materials, according to the flame sheet model, higher q-values cause stabilization on the top of the reactor only.Partial oxidation due to a finite reaction rate analyzed in Section 4.1, stabilization can also be achieved on the bottom part of the reactor.

Improvements of the Flame Sheet Solution
The flame sheet approximation is based on the assumption that reactions are infinitely quick.For real fuels, oxidation and pyrolysis reaction rates are always temperature and concentration-dependent, and the flame sheet approximation is typically not valid.Indeed, rates for heterogeneous reactions are intrinsically slower than homogeneous reactions, and gasification reactors are usually operated at low temperatures.

The Role of the Finite Reaction Rate
Several processes determine reaction rates for heterogeneous fuel conversion, namely, bulk diffusion, pore diffusion, and kinetic (see [61]).A pretty general relation containing the three processes mentioned above, can be obtained in the case of unity reactant reaction order: with η being the effectiveness, limiting the kinetic reaction due to pore diffusion, and D and K being the rate per unit area of particle surface due to bulk diffusion and kinetic, respectively.The kinetic reaction rate K is usually written in Arrhenius form with an exponential dependency on the temperature, while the bulk diffusion reaction rate D has a power low dependency.Both reaction rates linearly depend on the reactant concentration, because of the unity reaction order for the oxidizer: The effectiveness η can be expressed as a function of the Thiele modulus T h , defined as the ratio between the kinetic reaction rate and the diffusion rate [61]: The effectiveness is much smaller than one for particles whose pore radios are small enough (and long enough) to create a resistance to diffusion of the oxidizer into the particle's interior.From the phenomenological point of view, this limitation can be seen as a reduction in available area and a decrease in the activation energy as well as the oxygen reaction order [61].In the limit of small effectiveness (large Thiele numbers), the phenomenological activation energy E is half of the intrinsic activation energy.This limit is approached for larger particles, with diameters larger than several hundred micrometers.The diffusion rate limits the total rate at higher temperatures and for large particles with a diameter larger than several millimeters.
A finite rate reaction alters in a few aspects the thermodynamic behavior of the reactor seen until now.Firstly, the oxidation and pyrolysis zones have a finite width.The corresponding temperature profiles have been already discussed above.Secondly, the fluxes φ i and the composition of both the solid and the gaseous phase become a continuous function, quickly varying inside the reaction regions.The solution of the reaction zone structure requires the solution of the mass balance equations coupled with the energy balance equations.
Figure 6 shows temperature profiles for TC1 in which the oxidation rate gives rise to a finite reaction zone width used to label the profiles.The parameters of the reactions have been adjusted so that the center of the oxidation zone is positioned at x * = 0.1.The broadening of the oxidation zone is evident, but the most important ones to note are the changes in the temperature level.All the temperatures are generally lower, both in the oxidation region and in the region above it, with a decrease in the solid temperature (the red line in the figure) before the onset of the oxidation, identified by a temperature still denoted as the ignition temperature T ig .Similar behavior can also be observed in the other thermal cases, but the profiles are not reported here.The possibility of partial oxidation has an higher impact on the thermal structure since partial oxidation implies that less energy is released into the reactor.In fact, as shown before, the amount of energy injected into the system has an impact on the position of the oxidation zone.It has been shown in Section 3.2 that the presence of the oxidation region in the bottom part of the reactor is possible only in a limited range of the non-dimensional calorific value q.As this parameter assumes a larger (and more realistic) value, the only stable position appears to be the upper part of the reactor (see Figure 5).A finite rate for oxidation opens up the possibility of partial oxidation with some fuel leakage from the bottom part of the reactor.The energy released in the reactor is then reduced, and an oxidation front can also be stabilized near the reactor bottom.
Figure 7 (left) shows the position of the oxidation zone near the bottom of the reactor (x = 0) as a function of the non-dimensional calorific value q for several widths of the oxidation region.For reference, Figure 5 contains the same information derived using the flame sheet approximation.The lowest green curve in Figure 7 presents the results using the flame sheet approximation already given in Figure 5.As seen using the flame sheet approximation, at high non-dimensional calorific values q, the temperature is high enough that the oxidation region is pushed out of the reactor.In those conditions, no oxidation can be stabilized at the bottom of the reactor.When a finite rate reaction is applied, in the same conditions, partial oxidation can cause a reduction in the energy released inside the reactor with a decrease in the peak temperature, and stabilization can be achieved for larger q.The temperature profiles of the solid phase in the vicinity of the bottom of the reactor (x = 0) are shown in Figure 7 (right).The green curve presents the temperature profile obtained by the flame sheet approximation with a sheet positioned exactly at x = 0.The decrease in temperature due to the partially oxidized material is visible.
Finally, a less important consequence of partial oxidation is the possible change in ratio between gas and solid fluxes at the bottom part of the reactor, so that the thermal condition (TC) can be different than the one expected from the flame sheet approximation.The relation between the factors k and the fluxes has already been presented in Equation ( 67) and it is here reported again for clarity: At x = 0 and x = 1, using a new parameter C b , defined as the oxidized carbon in the process, and considering that φ s (0) = 1 − C b and φ g (1) = 1 + g • x, can be written: In the limit of no combustion, nothing happens inside the reactor, mathematically expressed by: The analysis of both, Equations ( 85) and (86) demonstrate the following:

The Role of the Effective Solid Conduction
As shown in Table 1, the applications considered in this article are characterized by small values of the parameter δ with respect to unity.The mathematical problem is, therefore, a advection-dominated problem.For completeness, however, in this section, we briefly discuss what happens if the conduction in the solid phase becomes significant.Not all the details will be given, but only a few summary points.
In Figure 8, on the left diagram, the temperature profiles, and on the right, the concentration C of the reactive fraction of the solid phase are presented.The profiles are calculated for a fixed stoichiometric number λ = 1.5 with δ as a variable parameter.The diagrams correspond to results for the TC1, but similar behavior is also found in other thermal conditions.Shortly, the following observation can be made when δ increases: 1.
The heat is removed from the highest temperature region and stored in the central part of the reactor; 2.
For sufficiently large δ, the reactor freezes, and the oxidation zone is quenched; 3.
The ignition temperature increases with the critical parameters α and q; 4.
The oxidation region widens, and if partial oxidation occurs (as in the case shown in the figure), the burnout increases; 5.
The steady state position of the oxidation front moves towards the upper part of the reactor.
In pyrolysis mode (under stoichiometric conditions) a similar behavior can be observed with a few minor differences.

The Role of Pyrolysis
Pyrolysis affects reactor stabilization usually to a lesser extent than oxidation.The pyrolysis reaction rate is expressed in Arrhenius form: where T g is the temperature of the gas phase (see the discussion in Section 2.2).The parameters A and E are fuel-specific, and their non-dimensional values are presented in Table 1 for both materials, which are used in the EL process and in the MA process.Pyrolysis reactions primarily cause changes in the fuel material's structure, composition, and reactivity (see, for example, [62-64]).However, the pyrolysis reaction velocity and its dependency on temperature have an impact on the thermal structure of the reactor.All the details leading to complications will be skipped in this article, mentioning only some aspects relevant to the reactor stability.Following the main line of this paper, pyrolysis causes changes in the following balances: 1.
Both mass fluxes φ i are changed.The gas-phase flux increases after the region where pyrolysis is active, while the solid-phase flux decreases.The reactor is divided into three regions instead of the two regions as considered in the previous sections.

2.
The total physical energy available in the system, previously released by oxidation only, will decrease due to endothermic pyrolysis reactions.Consequently, the temperatures decrease.
The ignition temperature is plotted against the position of the oxidation region in Figure 9.The oxidation is treated according to the flame sheet approximation, therefore the figure corresponds to Figure 4 where the different colors plot the results for the four thermal couplings α ∈ [1,3,6,10].Pyrolysis occurs at a temperature θ p = 0.4, and the other parameters are listed in Table 1.The thick lines show the results without pyrolysis, and the curve has already been presented in Figure 4.In contrast, the thinner lines are calculated for several fractions of solid phase converted, namely m p ∈ [0.02, 0.04, 0.06].
It is clear from the figures that in the case of TC1 (top left), no significant changes in the position of the oxidation region can be observed, while in the other thermal conditions, more substantial differences are observed, especially for larger thermal coupling α.Excluding the thermal case TC1, the position of the oxidation front moves downstream of the reactor while keeping the ignition temperature T ig constant.Materials with larger T ig are allowed to stabilize at the lower part of the reactor.In contrast, due to decrease in temperature at the upper part, no thermal stabilization is possible for larger x.

Discussion of the Results
Before discussing the application of the results to the two study cases (see Section 5.5), some important aspects of the model need to be discussed.First, the problem of the stability of the obtained solutions is addressed.Secondly, three particular aspects of the modeling are discussed: 1.
The effects of mathematical description of the pyrolysis reactions; 3.
The difference between ignition temperature in a co-flow and counter-current flow configuration.

The Stability of the Steady Solutions
The analysis presented in the previous sections demonstrates that the oxidation front is steady, for most thermal conditions, in two positions, one at the top and one at the bottom of the reactor.The investigation of the stability of these two thermal states is performed by solving the full unsteady Equations ( 31)- (36).As discussed in Section 2.4, the solution of the equations mentioned above has been achieved using a finite difference schema with spatial discretization performed with the spectral method.
The first and foremost significant result is that a dependency on the initial air stoichiometry ratio λ appears.Fuel-rich reactor (gasification mode) configurations are stabilized at the bottom of the reactor while fuel-lean reactors (combustion mode) at the top of the reactor.When one of the reactants is less available for reactions, an oxidation front is stabilized near the source of that reactant.This behavior is mathematically explained by observing the speed of the oxidation front, which is a function of the ratio λ, and it is negative for λ < 1 and positive for λ > 1.Thus, in combustion mode, the front moves upwards, while in gasification mode it moves downwards.
The proof of this statement is straightforward when the reacting zone is thin compared with the length of the reactor.Unlike the flame speed in the conventional combustion cases (see for an application [54]), in which the reactants are perfectly mixed before the onset of reactions, in opposed flow reactors, the movement of the oxidation front is determined based on the mass balance of char and oxygen only.
From the balance equations for the char (Y s,C is the mass fraction of the char in the solid phase) and the oxygen (Y g,O2 is the mass fraction of the oxygen in the gas phase), Equation ( 35) and ( 36) can be written in their conservative form: Subtracting the first equation from the second one (after multiplication by g): If the profile of each quantity involved (here shortly denoted as Ψ a ) does not change form but translates with a constant velocity U, Ψ a , it can be expressed using only a single variable ς = x − U • τ as Ψ a (x, τ) = Ψ a (x − U • τ).The time derivative and the spatial derivative are expressed using ς in the following way: If the reaction zone is thin compared with the reactor length, at x = 0 no char is present, and at x = 1, no oxygen is present, regardless of the stoichiometric coefficient λ.By integrating Equation (91) between 0 and 1, and imposing the boundary conditions at both extremes, the following relation can be obtained: which describes the behavior of the front as a function of the air ratio mentioned above.Despite the drastic assumptions of the previous calculation, the relationship (93) can also be successfully applied in the two study cases under consideration.The contour of the temperature in the plane x − t is presented in Figure 10 for both λ = 1.5 (a) on the left, and for λ = 0.5 (b) on the right.The MA configuration [51] is used since the oxidation region is particularly thin compared with the reactor size and the approximations better holds.In both diagrams, after ignition at t = 0 a transient period is present, which lasts differently for the two configurations.Specifically, when λ = 1.5, the ignition region is carried down by the solid material, and only after the front has reached the lower part of the reactor, does the true oxidation front move to the stable position located at the top of the reactor.In the case of λ = 0.5, the oxidation front is immediately transported to the lower part of the reactor where it stabilizes.The comparison of the results of this calculations and the values given by Equation (93) shows an excellent agreement.
The oxidation front is stabilized when the thin reaction zone approximation is no longer applicable and the the presence of the boundary conditions becomes important.In particular, temperature becomes an essential variable and stabilization is achieved when the temperature decreases enough to reduce the rate of the oxidation reaction.When λ is equal to one, every stationary position identified earlier is stable, and the position at which the reactor stabilizes depends on the mode of ignition.
Before closing the section, it is worth emphasizing that the air ratio λ is calculated here (see Equation ( 29)) using the char and the oxygen available at the respective boundaries.In a real gasifier, this number is not fixed for two reasons.Firstly, because of low temperatures, the pyrolysis reaction rates responsible for char production can be low, and the expected amount of char cannot be formed (see later for more details).Secondly, and most importantly, gasification reactions (C reacting with CO 2 and H 2 O) can be activated, reducing char concentration in the oxidation region.In both cases, the air ratio λ will increase, possibly changing the stoichiometry of the oxidation reaction.If one of both mechanisms is efficiently active, the oxidation front will switch from the lower to the upper part, drastically changing the reactor's temperature profile.

Oscillations before Flame Quenching
In the previous sections, it has been shown that for the reactor stabilization, it is required a minimum specific heat q min and its value depends on the thermal coupling α as shown in Figures 5 and 7.If the non-dimensional heating value is slightly smaller than q min (α), oscillations can take place, as shown in Figure 11.The treatment of the oscillations, presented in this section, does not want to be an exhaustive investigation, but it serves only as a reference to complete the description of the thermal behavior.
Figure 11 presents the time variations of the solid temperature at the outlet (θ s (x = 0) in blue) and the maximal solid temperature (in red) for two different thermal couplings.The left diagram represents the results for reactors with smaller thermal coupling (in the figure α = 5), while the right plot represents the results for a higher thermal coupling (in the figure α = 30).After a transient period at short non-dimensional times, the proper oscillations appear.The data in the figure arise in a fuel-rich reactor (pyrolysis mode), but similar oscillations are also noticed in a fuel-lean reactor (combustion mode).In all cases, the reactor is dynamically stabilized by the re-ignition of the fuel in a position further away from the position x = 0 when the solid temperature is not high enough to maintain a stable oxidation.During the quenching of the oxidation region (the red curve and the blue curve in the left plot are superimposed), less gaseous products are formed.Since the gas is colder than the solid, the cooling effect for the solid phase decreases.The temperature of the solid phase tends to slightly increase (the red curve bends up while the blue curve still descends), and, eventually, the ignition temperature is reached away upstream (the red curve is higher than the blue one).
With an increase in the thermal coupling (see Figure 11 on the right), the period of the oscillations increases, and the amplitude decreases.For a tighter thermal coupling, no oscillations are allowed anymore.As an example, with the parameters used in this section, a maximum value of α = 50 has been found for the maximum thermal coupling.
If q further decreases below a second critical value, the oxidation reactions quench entirely, and the reactor freezes.
Di Blasi has also found oscillatory behavior in [48] where the duration of the cycle is in the range 25 to 100 ms.In Figure 11, the cycle has period of ∆τ = 1.4 (left) and ∆τ = 0.25 (right), both on the same order of magnitude of τ = 1.The variable τ is a non-dimensional time normalized using Equation ( 21) therefore the period of oscillation in real time depends on the mass flow of the solid material and the reactor height.In the case of EL operating condition τ = 1 correspond to t = 15 h while in the MA experiment correspond t = 1 h.
The cause of these "dynamic patterns of the combustion/gasification zone" (as cited in the reference) has been identified in "the variations in the rates of solid/gas heat transfer".In the present work, the parameter identified by Di Blasi is the coefficient h(ε) in Equations (3) and ( 4).With the mathematical apparatus developed in this work, the identified parameter is part of the non-dimensional α number.Di Blasi discussed the sequence of events leading to an oscillating behavior in terms of particle diameters, basically caused by a reduction in the particle size at the minimum value.This explanation is partial since there is no variation of diameters in the present model.In a more general explanation, it has been proven that the cause of the oscillations lies in the relationship q < q min (α); therefore, it lies in a combination of thermal characteristics of the feedstock (q) and thermal coupling between the phases.

On the Modeling of Pyrolysis
Until this section, the fuel has been modeled to be composed of an inert fraction, a fraction that takes no part in the pyrolysis (the char), and a fraction that, thermally treated, is transformed into a gas (the volatile fraction).The fractions of inert, char, and volatile are fixed at the beginning of the calculations, independently from the thermal history of the solid material.In particular, a fixed amount of reactive char is always available in the oxidation region.In opposition to this simple model, there is experimental evidence (see [65,66] for recent results) that the amount of volatile depends on the sample thermal history.A theoretical explanation of this evidence, also supported by arguments for chemical structure changes, is based on the fact that char, the fraction of the solid phase available for oxidation, is formed during the pyrolysis process itself.As a consequence, an alternative, more accurate but still simple modeling approach, assumes that solid material (F) available in the raw material, during pyrolysis, undergoes the following splitting: where char C and gaseous volatile V are formed.As before, the reaction rate is still described by a rate in Arrhenius form.Variations of this model allow the dependency of the volatile yield on the heating rate [67].However, the effects of those models on the thermal structure of the reactor are not deeply investigated here.
Using this more complex pyrolysis model, it is not guaranteed that all the fuel F is converted into reactive char C and volatile V when the temperatures are too low.If the pyrolysis is modeled this way, reactors running with thermal parameters near the critical values found in the previous sections are unstable.In particular, if the pyrolysis zones and the combustion zone are not clearly separated (as in the EL setup) the oscillations observed in Section 5.2 disappear.

The Ignition Temperature T ig
The last point to be discussed concerns the ignition temperature T ig .As already mentioned, this temperature cannot be identified with the ignition temperature found in a co-current flow configuration or a well-mixed reactor.In fact, in a counter-current reactor, the char fraction and the oxygen fraction profile, which appear in the reaction rate, differ from the same profile in a parallel flow reactor.
Figure 12 presents a typical temperature variation of a solid material slowly heated up (T Ramp in the figure).The solid-phase temperature is shown for a co-current configuration (the blue lines) and a counter-current configuration (the magenta lines).In both cases, the horizontal lines identify the ignition temperature for the runaway.The figure shows that the ignition temperature is higher in the counter-current flow configuration than in a parallel flow configuration.In the same figure, the product of the fraction of char and partial pressure of oxygen is also plotted, showing the difference that causes an increase in the ignition temperature.

Application to the Test Cases
The simplified model described until now can be applied to two processes, the EL process [38] and the MA process [51] (see Table 1 for the parameters of both processes), both used as test cases.Figure 13 presents the temperature profiles for the MA process (Left in the figure) and for the EL process (right in the figure).A model based on the flame sheet approximation (the thick lines in the figure) is compared with the one based on finite rate chemistry with effective conduction in the solid phase (the thin lines in the figures).In case of the MA process, a measured temperature profile is also available [51], and temperature data are reported in the figure mentioned above as dots, together with the predictions from the models.The good agreement between the two models, already noticed previously in the discussion of Figure 1, validates the utilization of the flame sheet approximation.This validation is reinforced by the good agreement with the measured temperature in the case of the MA process.It is important to point out that both models used for this discussion remain simplified models, with the consequence that no perfect agreement can be reached.Specifically the temperature dependency of the specific heats of both the gas and the solid material has been disregarded.For a more accurate modeling approach with a tighter agreement with measurements the reader is referred to the original paper [51] and to the successive analysis [30].In both processes, the oxidation front is stabilized at the bottom of the reactor, which is a consequence that both reactors are gasifiers and they are operated in under-stoichiometric conditions.
A notable difference between the two temperate profiles is that, in case of MA reactor, the pyrolysis region is positioned at the top of the reactor, while in the case of the EL process, the pyrolysis is just above the oxidation front.This observation is a consequence of the different thermal conditions TC for the reactors and the difference in the nondimensional heating value q.In fact, the oxidation zone of the MA experiments is in thermal condition TC1, while the reactor proposed in EL is in thermal condition TC2 (see Table 1).Both characteristics can be traced back to the presence of a considerable quantity of inert material in the EL process.
Recalling the analysis of thermal profiles in Section 3.1, a negative value of the parameter k 2 imposes an exponentially increasing solid temperature profile (decreasing x, i.e., following the movement of the solid).The exponential increase forces large temperature gradients with the consequence that conditions, where the pyrolysis is active, are found near the oxidation front, as seen in the EL case.On the contrary, a positive value of the parameter k 2 imposes a milder increase in the solid temperature, and significant temperature changes are observable only in the upper region of the reactor.Conditions for pyrolysis are found near the upper part of the reactor, as seen in the MA configuration.This observation is general, and the two proposed reactor configurations are just a realization of the underlying different thermal structures.
To make this statement clearer, Figure 14 presents the position of the pyrolysis zone as a function of the oxidation zone for calculations performed with parameters similar to MA and EL, while the pyrolysis temperature is θ p = 0.4 for both processes.The curves are calculated with four different thermal couplings (α = [1,3,6,10]) and c = [0.6,0.7, 0.8, 0.9] for the EL configuration and c = [1.4,1.5, 1.6, 1.7] for the MA configuration.In case of a process similar to the MA process (on the left), a slight influence of the thermal coupling α is observed, the pyrolysis zone is positioned for all the coupling investigated at the top of the reactor (see also Figure 9 and the corresponding discussion).The plot on the right in Figure 14, representative of a thermal condition similar to the EL process, the thermal coupling α does not influence much the relative positions of both reaction zones.Moreover, the pyrolysis zone is tightly correlated with the position of the oxidation, and both reaction zones are always positioned near each other.
In industrial applications, the proximity of both regions can be a disadvantage.In fact, in the case of low char reactivity, oxygen may diffuse from the oxidation region and reacts with the pyrolysis gas, reducing the efficiency of the gasifier.The conditions and the consequences of overlapping between the oxidation and the pyrolysis regions in the EL process have also been proven by simulations with a more detailed model (see [30,41]) and the discussion is not here elaborated further.

Conclusions
The paper presents a simplified model to investigate the configurations of stationary conditions in a counter-current reactor.The stability of the stationary configurations are also analyzed.The flame sheet approximation is used at first because of its numerical simplicity and the model is used primarily to identify non-dimensional parameters responsible for the thermal structure of the reactor.In a second step, a more complex model with temperature dependent conductivity and reaction rates in the solid material are taken into account.The differential equations are solved using the finite difference method based on spectral approach.
The following conclusions can be listed after the analysis of the results: 1.
There are four possible thermal structures in a counter-current reactor; 2.
The constraints on the thermal structure allow, for most of the parameters, only two steady-state positions of the oxidation front, one in the upper part of the reactor and one in the lower part of the reactor; 3.
An increase in the heating value of the solid material or an increase in thermal coupling results in forcing both steady state positions to drift towards the extremities of the reactor; 4.
The position of the oxidation zone in the upper part of the reactor is stable in fuel-lean conditions (combustion mode), while the position in the bottom part is stable in fuel-reach conditions (gasification mode); 5.
A critical minimum heat value of the feedstock exists under which the reactor cannot be stabilized by an oxidation front; this value is primarily a function of the thermal coupling between the phases; 6.
Thermal oscillations can exist for heating values smaller than the critical heating value if enough char is available in the oxidation zone and if the thermal coupling is not too strong.
These unique results offer a theoretical basis for the interpretation of reactors operation, which is also useful for the conception of future reactors but the model lacks the accuracy needed for detailed predictions.Not all the simplifications have been removed as, for example, some thermodynamic quantities are still held constant.However, a model capable of being detailed, must implement more advanced correlations following the investigation in [30,48] for example.In particular, the gasification reactions and the water-gas shift reaction must be implemented if the product composition is of interest, as shown in the previously cited work.Such a model has the ability to predict the thermodynamic structure for a given real material not only qualitatively, but also quantitatively.
From the scientific point of view, the following phenomena needs to be further investigated: 1.
Interaction between the combustion zone and the pyrolysis zone.In the previous sections, it has been demonstrated that there is an interaction between the two regions.As a unique result, it has been shown that the thermal structure of the reactor is primarily responsible for this interaction, but a more detailed analysis taking into account more precise reaction rates needs to be performed.

2.
The existence of possible oscillations was discovered in [48] and in this work a more accurate theoretical explanation has been provided.Although probably not relevant for industrial applications, the study of oscillations in an analysis of reactor dynamic phenomena is of interest.In particular, the amplitude and frequency of the oscillations depend on the operating conditions of the reactor, the property of the feedstock, and the detail of the pyrolysis reactions.

Figure 1 .
Figure 1.Typical temperature profiles present in a counter-current reactor in which the ignition is on the lower side of the reactor x = 0.4.The red lines indicate the temperature of the solid phase and the blue lines indicate the temperature of the gas phase.• Thermal state TC1, in which k 1 < 0 and k 2 < 0, at the top left part of Figure 1.• Thermal state TC2, in which k 1 > 0 and k 2 > 0, at the top right part of Figure 1.• Thermal state TC3, in which k 1 > 0 and k 2 < 0, at bottom left part of Figure 1.• Thermal state TC4, in which k 1 < 0 and k 2 > 0, at the bottom right part of Figure 1.

Figure 2 .
Figure 2. Existence conditions for the four thermal cases as a function of the char concentration in the fuel Y 1 s,C and the air excess ratio λ.The black and the thin magenta lines correspond to k 1 and k 2 , respectively, while the think lines present the borders for the various thermal conditions.Plot (a) is representative of a c < 1, and Plot (b) is representative of c > 1.

Figure 3 .
Figure 3.Typical temperature profiles expected in the counter current gasifier for α = [2.0,4.0, 6.0, 8.0].Pyrolysis is not present, and the reactor works under Thermal Condition 1 (a) and Thermal Condition 3 (b).The red lines indicate the temperature of the solid phase and the blue lines indicate the temperature of the gas phase.

Figure 4 .
Figure 4. Ignition temperature as a function of the position for the four thermal conditions.The curves are presented for four different values of the parameter α [10.0, 6.0, 3.0, 1.0] and for different values of the parameter c.

Figure 5 .
Figure 5. Runaway position as a function of the non-dimensional fuel heating value q.(a) The plot on the left shows the runaway position for several α parameters and constant T ig = 0.7, while (b) the plot on the right shows the runaway position for several ignition temperatures at a fixed thermal coupling α = 2.The thin lines are calculated with c r = 1.5.

Figure 6 .
Figure 6.Gas and solid phase temperature profiles for several reaction zones widths.The calculations have been performed in TC1.The blue lines show the gas temperature and the red lines the solid temperature.

Figure 7 .
Figure 7.The same as Figure 5 with a finite rate oxidation.
the reactor is in TC1 for all the values of C b ; • If w > 1 and c > 1 we have the sequence TC2 → TC3 → TC1 for a decreasing C b ; • If w > 1 and c < 1 we have the sequence TC2 → TC4 → TC1 for a decreasing C b .

Figure 8 .
Figure 8.(a) Solid temperature profiles and (b) char fraction profiles for different choices of the parameter δ.The calculation has been performed in TC1 and α = 10.

Figure 9 .
Figure 9. Ignition temperature as a function of the position for the four cases when pyrolysis is present in the process.The other conditions are the same as for the calculations presented in Figure 4.

Figure 11 .
Figure 11.Time variations of the solid temperature at the outlet (θ s (x = 0) in blue) and the maximal solid temperature (in red).(a) results for thermal coupling α = 5; (b) results for thermal coupling α = 30.

Figure 12 .
Figure 12.Comparison between the ignition temperature in counter-current conditions (the magenta lines) and in co-flow conditions (the blue lines).The product Y s,C • p O2 is also shown.

Figure 13 .
Figure 13.Profiles of solid-phase Temperature (in red) and gas-pahse temperaturte (in blue) (a) of the MA process and (b) of the EL process.The thick lines are calculated with the flame front approximation while the thin lines are calculated with the finte rate model.

Figure 14 .
Figure 14.Position of the Pyrolysis region as a function of the position of the combustion region for (a) the MA process and (b) for the EL process.The curves are calculated with four different thermal couplings (α = [10, 6, 3, 1], the black, blue, rad and magenta lines respectively) and c = [0.6,0.7, 0.8, 0.9] for the EL configuration and c = [1.4,1.5, 1.6, 1.7] for the MA configuration.

Table 2 .
Parameters used in the calculation whose results are presented in Figure4.