1 Introduction

Integration of biomass gasification and a combined cycle power plant (Bio-IGCC) is an attractive option for balancing the future power grid, which is expected to have a large fraction of intermittent wind and/or solar power (Varone and Ferrari 2015). The advantage with Bio-IGCC is that it is carbon neutral, that it can be built in large scale (> 100 MW), and that it is freely dispatchable. Thus, it can be used to compensate for large amplitude dips in the overall power generation when the wind stops blowing or when the sun goes down as well as provide a reliable base power generation. Moreover, Bio-IGCC has a very high power efficiency (Saini and Pollock 2012), which is important both for economic reasons and to leverage the available biomass resource, which is limited by sustainability constraints.

The main difference between Bio-IGCC and conventional IGCC is the composition of the fuel that is burned in the gas turbine. The fuels that are commonly used in conventional IGCC are characterised by high Higher Heating Values (HHV) and moderate flame speeds while the syngas that is used in Bio-IGCC has significantly lower HHV and higher flame speed. In addition, there is a significant variability of the syngas composition in Bio-IGCC depending on the feedstock and the gasifier type. One way to enable the use of existing gas turbines in Bio-IGCC would be to convert the syngas to Substitute Natural Gas (Thunman et al. 2018). However, that would lead to an unacceptable energy loss compared to burning the raw gas directly. Hence, fuel flexible gas turbines that are capable of combusting all types of biomass syngas are an enabling technology for Bio-IGCC.

The data in Table 1 shows the gas composition for three candidate gasifiers, two entrained flow (Weiland 2015; Carlsson 2011) and one dual fluidised bed (Thunman et al. 2018) gasifiers, with three different feedstocks (black liquor, wood and bark). As can be seen from the table, both the raw gas major species H2, CO and CO2 and the minor species, including CH4, C2 hydrocarbons and tar species, differ significantly between the three different types of gasifiers. Other types of gasifiers and feedstocks can also be used in a Bio-IGCC power plant but the range of syngas in Table 1 is a good approximation of the possible composition space.

Table 1 Syngas properties from three biomass gasification processes (Thunman et al. 2018; Weiland 2015; Carlsson 2011)
Full size table

At present, the most interesting option for Bio-IGCC, in countries with a large forest product industry, is Black Liquor Gasification (BLG) since it can be operated in large scale (100–500 MW thermal) at a single pulp mill and the global potential is very large (Jafri et al. 2019). BLG is based on oxygen blown entrained flow gasification of small (~ 100 μm) black liquor droplets at approximately 1050 °C and 30 bar. The syngas is quench cooled with water to avoid clogging. The syngas from oxygen blown and pressurised black liquor gasification (see Table 1) is characterised by a very low concentration of tar and soot and an unusually high concentration of H2S. The major species concentrations are almost equal with about 1/3 each for H2 and CO2 while the molar concentration of CO is slightly lower at 28%. The only significant minor species are CH4 and H2S at about 1–2%. Ideally, the H2S should be removed before combustion, to simplify flue gas cleaning, while the CH4 should be left in the gas to maximise its heating value. This can be achieved with existing amine wash technology, which at the same time can remove all or part of the CO2, to increase the heating value of the syngas. The resulting heating value of cleaned syngas where all of the CO2 has been removed will increase by about 33% while the corresponding Wobbe Index will increase by more than 80% (see Table 1). However, the BLG syngas Wobbe index is still only about 1/3 of that in methane, even after removal of all the CO2 in the syngas.

Solid biomass gasification is also an important option, since this makes it possible to expand the feedstock base and to implement large scale Bio-IGCC in countries without pulping industry. The first example of solid biomass gasification in Table 1 is the oxygen blown and pressurised Entrained Flow Gasification (EFG) of finely ground (~ 200 μm) dry fuel particles at approx. 1400 °C and 10 bar. The syngas produced is quenched cooled with water to avoid clogging. A small amount of nitrogen is added in certain sensitive positions for safety reasons. The gasifier operates at slagging conditions with molten ash and the resulting high temperature typically results in a syngas with a CO concentration above 40% while the H2 and CO2 concentrations are around half of that. This syngas is similar to Black Liquor (BL) syngas in that it contains very little tar. The only significant minor species is CH4 and there are only trace amounts of sulfuric species in the syngas due to the low sulphur content in the feedstock. The only cleaning that would be necessary before combustion in a gas turbine would be removal of large particles and aerosols, which can be done efficiently with existing technology. The resulting Wobbe index for the clean syngas is about 1/5 of that in methane.

For the Dual Fluidised Bed gasifier (DFB), hot sand from an air-fluidised bed combustor is fed to a steam-fluidised bed reactor together with dry wood chips or pellets. The fuel is partially gasified in the first bed and the residual char is fed back to the combustor together with cooled sand to generate the heat for the process. It operates at non-slagging conditions with dry ash and the necessary low temperature (app. 800–900 °C, depending on the ash melting temperature) results in a hydrogen rich gas with up to 10% methane. It has the advantage that although it uses air as oxidant the syngas is virtually nitrogen free, which reduces the CAPEX of the gasifier by eliminating the need for an air separation unit and increases the heating value compared to traditional air blown gasification (Higman and Burgt 2011). However, the syngas would need both tar reforming and cleaning before it can be used in a gas turbine, which can be done with existing hot gas filtering and catalytic tar reforming (Pfeifer and Hofbauer 2008). The resulting values for clean syngas in Table 1 are calculated assuming that all the tar species are converted to syngas and all minor species are eliminated in the same operation. The corresponding Wobbe Index for the clean syngas is 23.62 MJ/Nm3, which is about 1/2 of that in methane.

The Wobbe index differences make it clear that the general behaviour of biomass syngas, when it is used as a fuel in a gas turbine, will be significantly different to that of typical fossil fuels, e.g. natural gas. Moreover, the differences in laminar flame speed of the syngas compared to fossil fuels will have a direct impact on the combustion behaviour (Taamallah et al. 2015). As an example, the flame speed for BL syngas, with and without CO2 removal, is compared to methane for different equivalence ratios in Fig. 1. The flame speeds and equivalence ratio were calculated with the software Cantera for all cases and considers oxidation of C, H and S (which is not included in the present paper), all other species are considered inert. Looking at stoichiometric conditions, the syngas derived from BLG, has a flame speed of up to 118 cm/s, whereas the flame speed of methane in the same conditions is 38 cm/s.

Fig. 1
figure 1

Laminar flame speed versus equivalence ratio for syngas from black liquor gasification (see Table 1), both with and without removal of CO2. The corresponding flame speed for pure methane is also included for comparison. The flame speeds were calculated with Cantera (Goodwin et al. 2020)

Full size image

The implication of the data in Table 1 and Fig. 1 is that gas turbines intended for Bio-IGCC must be optimised with respect to the expected syngas properties. Among the tasks of this optimisation the adaptation of the combustor hardware to maintain a reliable and low emission combustion is a very important one. This optimisation can be done with experimental methods but experiments in full scale are very costly and it is therefore beneficial to do most of the optimisation with computer simulations and to do a small set of advanced experiments to finetune the design and to provide data for optimisation of the turbine section. The current trend in computer simulations of combustion is to use Large Eddy Simulation (LES). In the context of combustion in a premixed gas turbine environment, it is common to use high swirl to stabilise the flame. The resulting flow field is very complex and challenging for turbulence modeling since it involves phenomena, e.g. vortex breakdown and precessing vortex core (Syred 2006) that are difficult to model. Hence, simulations of combustion of a complex fuel under gas turbine conditions are very demanding and time consuming, even on a super-computer. An attractive alternative that requires significantly less computational power is to solve the Reynolds Averaged Navier–Stokes (RANS) equations combined with a suitable combustion model (Benim and Pfeiffelmann 2019). This method does not resolve all the fine details that a LES simulation does, neither is it capable of predicting all of the phenomena that a LES simulation can capture, e.g. flashback and thermoacoustic instabilities. However, the significantly faster turnaround times compared to LES make it possible to try many more design alternatives and to find an approximate global optimum design that can be used as the starting point for LES simulations or experiments.

Regardless of whether the approach is based on LES or RANS it is necessary to model the chemical reactions. Mansourian and Kamali investigated syngas turbulent combustion using OpenFoam (Mansourian and Kamali 2017). They employed Renormalization Group – Large Eddy Simulation (RNG-LES) turbulence model, paired with kinetic mechanisms like Skeletal and GRI-MECH 3.0 (Smith et al. 2018). For the turbulence–chemistry interaction they used the extended Eddy Dissipation Concept extinction from Aminian et al. (2016). Their findings suggest that GRI-MECH 3.0 produces more accurate results than skeletal kinetic mechanisms by providing fewer computational errors in the simulations if the fuel contains species like H2, CO, CH4, CO2 and H2O.

Detailed chemical schemes are computationally expensive but can be economically implemented through Flamelet combustion models. These models assume that the chemical states of turbulent flames are similar to those of laminar flames, and pretabulate the chemistry using 0D or 1D geometry (Oijen et al. 2016). In addition, the state of the flame can be expressed based on a few control parameters (mixture fraction, reaction progress variable and enthalpy). Hence, there is no need to solve transport equations for all the species in the detailed mechanism, which will significantly reduce the computational time and memory requirements.

The aim of this study is to compare the combustion behaviour of a wide range of syngas fuels that can be expected from biomass gasification, with a special focus on how the laminar flame speed influences the flame shape and position. Before the comparison can be made the model must be validated against experiments. The TECFLAM (Schmittel et al. 2000) swirl stabilised burner was chosen as suitable test case, in the premixed, atmospheric configuration (Meier et al. 2000; Schneider et al. 2005; Ketelheun et al. 2013) (see Fig. 2). The flame in this setup is stabilised over a water-cooled bluff body, with a premixed fuel–air mixture issuing from an annulus concentric to the bluff body. Moreover, there is an annular co-flow of air with an axial velocity of 0.5 m/s, to shield the flame from external disturbances. There are ample experimental data available in three operating conditions (30, 90 and 150 kW). In the present paper, the case of 30 kW, with a Reynolds number of 10.000, an equivalence ratio of 0.833 and a swirl number of 0.75, is selected for validation (Schneider et al. 2005). A comparison of different eddy viscosity models (k–ω SST, Realizable k–ε 2 layer), Reynolds stress models (RSM), and combustion models, Flamemelet Generated Manifold (FGM) and Eddy Break Up (EBU) is made in Sect. 3.1. After the selection of the optimum model in Sect. 3.1, the resulting model is used for a series of simulations with different fuels to characterise the flame response to different flame speeds. Finally, the results are discussed with respect to how the flame shape is affected by the type of fuel and how the laminar flame speed can be used as an indicator of the flame behaviour.

Fig. 2
figure 2

Schematic of the premixed TECFLAM configuration adapted from Ketelheun et al. (2013)

Full size image

2 Numerical Setup

All simulations were carried out with the commercial software STAR CCM + (Simcenter STAR-CCM+ User Guide 2020), which in the present setup uses a finite volume discretisation method and a segregated flow model. This flow model solves the equations for velocity in three directions as well as pressure sequentially and then iterates over the equations to achieve the necessary coupling between them. A collocated variable arrangement is used combined with the SIMPLE algorithm (Simcenter STAR-CCM+ User Guide 2020; Ferziger et al. 2020). A tight coupling between pressure and velocity in the same cells is achieved by using Rhie-Chow (pressure weighted) interpolation in the convective terms (Ferziger et al. 2020). The code can use both structured and unstructured grids and in the present case, an unstructured polyhedral mesh was used. The discretisation of the equations was done with second order accurate approximations. Experimental velocity data and the turbulent kinetic energy that were used as boundary conditions at the premixed inlet were extracted from measurements taken 1 mm above the burner inlet (Schneider et al. 2005). The turbulent dissipation (ε) and the specific dissipation rate (ω) that were also used as boundary conditions, for respective models, for the premixed inlet were calculated from the data of Schneider et al. (2005). The experimental data at the inlet were interpolated to the numerical grid and used, together with an assumption of perfect premixing between air and fuel and constant wall temperatures, as boundary conditions for the simulations. Hence, the interior volume of the burner was excluded from the simulations. The simulation domain can be viewed in Fig. 3. Assuming rotational symmetry, a 30° slice of the full 360° domain was used together with periodic boundary conditions. This approach was selected to reduce the computational load but still resolving the 3-dimensional velocity field accurately. The inlet inner and outer radius are 15 mm and 30 mm respectively, the burner mouth extends to 100 mm radially. The height of the domain is three times the co-flow width (3 × 600 mm) to avoid recirculation through the outflow boundary, which would be incompatible with the Neumann boundary condition at the outlet. Room temperature at 300 K is set for both the premixed inlet and co-flow. The utilised mesh comprises of more than 2 million cells and is refined from the premixed inlet up to 150 mm upstream with a 160 mm radius.

Fig. 3
figure 3

The 30° slice of the cylindrical geometry that was used for the simulations is shown on the right. On the left, a magnified view of the domain near the burner inlet is shown. The inlets are shown in red

Full size image

The boundary conditions for the velocity at the inlets were interpolated from experimental data (Schneider et al. 2005). Field functions for the premixed inlet were written based on the experimental data and radial distance to simulate the experiment as close as possible. On the solid walls blended wall functions (Simcenter STAR-CCM+ User Guide 2020) were applied. On the water-cooled bluff body a constant temperature was used in the wall function, while the outer boundary of the simulation was treated as an adiabatic solid wall with a no slip condition. On the outlet a Neuman zero gradient condition was applied for all variables except the pressure which had a constant value Dirichlet condition. The ambient pressure was atmospheric and all gases were assumed ideal and incompressible. All combustion simulations were started from a converged cold flow solution to reduce the total computational time.

2.1 Combustion Models

For methane combustion two models were investigated, the standard EBU model (Poinsot and Veynante 2012; Magnussen and Hjertager 1977) and the FGM (Oijen et al. 2016). The EBU combustion model was selected as a baseline model for this study as it is a widely used industry standard. EBU is a reactive species transport model which assumes that the reaction rate is governed by the turbulent mixing rate. For standard EBU the reaction rate \(r_{F}\) expression is as follows (Magnussen and Hjertager 1977):

$$ r_{F} = - \frac{\rho }{{W_{F} }}\left( {\frac{1}{{{\rm T}_{{turb}} }}} \right)A_{{ebu}} \min \left[ {Y_{F} ,\frac{{Y_{o} }}{{s_{o} }},B_{{ebu}} \left( {\frac{{Y_{{p1}} }}{{s_{{p1}} }} + \frac{{Y_{{p2}} }}{{s_{{p2}} }} \cdots \frac{{Y_{{pn}} }}{{s_{{pn}} }}} \right)} \right], $$
(1)

where \(\rho\) is the density, \(W_{F}\) is the molecular weight of fuel, \({\rm T}_{{turb}} \) is the turbulent time scale (equal to k/ε when using the k–ε turbulence model), \(A_{{ebu}} \) and \(B_{{ebu}} \) are EBU model coefficients with values of 4 and 0.5 respectively, \(Y\) is the mass fraction, \(s_{o}\) and \(s_{{pi}}\) are:

$$ s_{o} = {\raise0.7ex\hbox{${\upsilon _{o} m_{o} }$} \!\mathord{\left/ {\vphantom {{\upsilon _{o} m_{o} } {\upsilon _{F} W_{F} }}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${\upsilon _{F} W_{F} }$}}, $$
(2)
$$ s_{{pi}} = {\raise0.7ex\hbox{${\left| {\left. {\upsilon _{{pi}} } \right|} \right.m_{{pi}} }$} \!\mathord{\left/ {\vphantom {{\left| {\left. {\upsilon _{{pi}} } \right|} \right.m_{{pi}} } {\upsilon _{F} W_{F} }}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${\upsilon _{F} W_{F} }$}}, $$
(3)

where \(\upsilon\) is a stochiometric coefficient, \(m\) is molar mass and \(W\) is the molecular weight. The subscripts o, F and p refer to oxidiser, fuel and products respectively. Furthermore, it uses limited complexity reaction mechanisms, hence, it tends to over- or under-predict certain species. The chemical scheme chosen for EBU simulations was the Jones–Lindstedt (JL), 6 species, 4-step reaction mechanism (Jones and Lindstedt 1988) shown in Table 2.

Table 2 JL reaction scheme (Jones and Lindstedt 1988)
Full size table

As mentioned in the introduction, flamelet models parameterise species based on progress variable, mixture fraction and enthalpy. The mixture fraction \(Z\) is defined as (Van Oijen and De Goey 2000; Lehtiniemi et al. 2006):

$$ Z = \frac{{m_{F} }}{{m_{F} + m_{o} }}, $$
(4)

where \(m_{F}\) and \(m_{o}\) are the mass of fuel and mass of oxidiser respectively. In a non-premixed configuration, a mixture fraction value of 1 would refer to the fuel stream and a value of 0 to the oxidiser stream. In this study, since there is an air co-flow, two streams were used, one for pure air and one for the premixed fuel corresponding to a mixture fraction value of 0 and 1 respectively.

For the definition of the progress variable, the chemical enthalpy formulation was used, which is based on the chemical formation enthalpies of each species (Lehtiniemi et al. 2006). The values of these enthalpies were obtained through the thermodynamic properties file imported alongside the reaction mechanism (Van Oijen and De Goey 2000; Lehtiniemi et al. 2006), in this case, GRI-MECH 3.0. Similar to the mixture fraction range, the FGM progress variable has values from 0 to 1 representing an unburnt or burnt condition respectively. For the progress variable (denoted as c) definition based on enthalpy, the equations:

$$ \rho \frac{{\partial \varUpsilon _{\iota } }}{{\partial t}} + \rho V_{x} \frac{{\partial \varUpsilon _{\iota } }}{{\partial x_{x} }} - \rho D\frac{{\partial ^{2} \varUpsilon _{\iota } }}{{\partial x_{x} ^{2} }} = \omega _{i} , $$
(5)

for species and

$$ \rho \frac{{\partial h}}{{\partial {\text{t}}}} + \rho V_{x} \frac{{\partial {\text{h}}}}{{\partial x_{x} }} - \rho V_{x} \frac{{\partial ^{2} \varUpsilon _{\iota } }}{{\partial x_{x} ^{2} }} = \frac{{\partial {\text{p}}}}{{\partial {\text{t}}}} + q_{R} , $$
(6)

for enthalpy, and finally the progress variable:

$$ c = \frac{{\mathop \sum \nolimits_{{i = 1}}^{N} h\left( {Y_{i} \left( Z \right),\tau } \right) - \mathop \sum \nolimits_{{i = 1}}^{N} h_{u} \left( {Y_{i} \left( Z \right),0} \right)}}{{\mathop \sum \nolimits_{{i = 1}}^{N} h_{b} \left( {Y_{i} \left( Z \right),\tau _{\infty } } \right) - \mathop \sum \nolimits_{{i = 1}}^{N} h_{u} \left( {Y_{i} \left( Z \right),0} \right)}} $$
(7)

are used. Where ρ is the density, p is the pressure, \(Y_{i}\) is the mass fraction of the ith species, D is the diffusion coefficient, \(V\) is the velocity, \(\omega _{i} \) the reaction rate of the ith species,\(q_{R}\) is the volumetric heat loss term, h is the enthalpy, \(\tau\) is the flamelet time and \(Z\) is the mixture fraction. Where the subscripts u and b refer to the initial unburnt and burnt states respectively. Equations (5)–(7) come from reference (Lehtiniemi et al. 2006), for more information the reader is directed there.

The progress variable source term, \(\dot{\omega }_{y}\), was calculated according to the Coherent Flame Model (CFM) (Meneveau and Poinsot 1991). The assumption in this model is that combustion occurs in the flamelet region and the reaction rate is defined as the product of the flame surface density, the laminar flame speed and the unburnt density.

$$ \dot{\omega }_{y} = - \rho _{u} S_{u} \Sigma , $$
(8)

where \(\rho _{u}\) is the unburnt density, \(S_{u}\) is the laminar flame speed and \(\Sigma\) is the flame area density. For all methane simulations the laminar flame speed was calculated through Gülders equation (Gülder 1991):

$$ {\text{S}}_{u} = {\text{ZW}}{\upphi }^{\eta } \exp \left[ { - {\upxi }\left( {{\upphi } - 1.075} \right)^{2} } \right]\left( {{\raise0.7ex\hbox{${{\text{T}}_{u} }$} \!\mathord{\left/ {\vphantom {{{\text{T}}_{u} } {{\text{T}}_{0} }}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${{\text{T}}_{0} }$}}} \right){\upalpha }\left( {{\raise0.7ex\hbox{${\text{P}}$} \!\mathord{\left/ {\vphantom {{\text{P}} {{\text{P}}_{0} }}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${{\text{P}}_{0} }$}}} \right)^{{\upbeta }} , $$
(9)

where Z, W, η, ξ, α and β are fuel dependent constants with values 1, 0.422, 0.15, 5.18, 2 and 0.5 respectively. For the syngas case the laminar flame speed was calculated from 1D simulations with Cantera (Goodwin et al. 2020) using the GRI-MECH 3.0 as the chemical mechanism.

Simulations were carried out both with and without radiative heat transfer to assess the importance of radiation in the flame. The radiation was modelled with the Discrete Ordinates Method utilising an improved grey gas approximation (Simcenter STAR-CCM+ User Guide 2020). The model used for radiation was the correlated k-distribution method, which is a refinement of the weighted sum of gray gases method. In the k-distribution model, a spectral reordering was used to define the absorption coefficients of H2O and CO2. Unfortunately, the inclusion of radiation in the model adversely impacts the computational time (24 h for convergence versus 8 h).

2.2 Selection of Simulated Cases

The syngas compositions used for the simulations were representative of BL (Raw and Clean), entrained flow of bark powders and Dual Fluidised Bed (DFB) gasification technologies, all according to Table 1. These gas mixtures cover the range of major and minor species concentrations that can be expected from biomass gasifiers that are suitable for use in Bio-IGCC. In addition, in order to test the hypothesis that the dominant parameter that controls the flame shape is the laminar flame speed, a hydrogen enriched methane mixture of 40% H2 enrichment was also included. This value of hydrogen enrichment was selected to yield the same laminar flame speed as for the entrained flow bark syngas, but with a significantly different HHV. In Table 3, the exact composition for each of the cases is presented as well as the adiabatic flame temperature, the HHV and the laminar flame speed at an equivalence ratio of 0.833. Besides the practical relevance of the fuels, the selection also results in a wide range of laminar flame speeds, ranging from 29 cm/s for pure methane to 88 cm/s for BL Clean.

Table 3 Composition and properties at an equivalence ratio of 0.833 for the gas mixtures that were selected for the simulations
Full size table

For the reactive simulations of syngas, and as is required by the CFM, the laminar flame speed had to be defined as an input. This was done in tabular form of laminar flame speed versus equivalence ratio. Figure 4 shows a comparison of the syngas cases and methane based on laminar flame speed from 1D simulations (Goodwin et al. 2020). Notice the similarity between 40% hydrogen enrichment of methane and entrained bark powder syngas for equivalence ratios around 0.83.

Fig. 4
figure 4

Laminar flame speed versus equivalence ratio for all the cases considered

Full size image

3 Results and Discussion

3.1 Validation

3.1.1 Cold Flow

A set of cold flow simulations was run to select the best alternative of three candidate turbulence models, based on the agreement with cold flow experimental data in the TECFLAM burner. The axial velocity from experimental cold flow data (Schneider et al. 2005) was used for the assessment of the models. Experimental data is shown in filled circles in Fig. 6, whereas red, green and blue solid lines represent the k–ω SST, Realizable k–ε 2 layer, and Reynolds Stress Turbulence Model, respectively (Wilcox 1994; Pope 2000). In Fig. 5 it is evident that all of the selected models agree well with the experiments in the vicinity of the mouth of the burner at 0 < z < 20 mm. However, from 20 mm upwards significant deviations from the experimental results become obvious for the RSM and Realizable k–ε 2-layer model. On the other hand, the k–ω SST model agrees quite well with the experiments, albeit with a slight under-estimation of the maximum axial velocity. The peak velocity at 60 mm above the premixed inlet is 20% lower than that of the experiments. For all the results presented, r is the radial position, Ro is the radius of the inlet and H is the axial position.

Fig. 5
figure 5

Cold flow comparison of different turbulent models at different axial locations. Notice the reversed flow close to the axis at all axial positions, which is a result of the swirling flow

Full size image

In order to interpret the reactive flow results it is instructive to have an understanding of the general behaviour of the non-reacting flow. Based on the chosen flow model, the flow pattern can be described as an annular swirling jet that issues from the burner exit. A short distance downstream from the burner exit the flow (see Fig. 11) forms a central recirculation zone enclosed by the annular swirling jet. On the outside of the annular swirling jet, near the premixed inlet, there is a second toroidal recirculation zone driven by the co-flow that separates at the outer edge of the horizontal surface between the premixed and co-flow inlets of the burner (see Fig. 11). The horizontal width of the recirculation zone is most likely proportional to the width of the horizontal surface.

The swirl number for the experiments and the simulations was also compared to facilitate the choice of a turbulence model. Table 4 below shows the comparison of the swirl number for different axial positions. The swirl number was calculated based on (Weber and Dugué 1992; Zhao and Weber 1992), by omitting the pressure term. From the comparison of Table 4, it can be seen that the k–ω SST model agrees satisfactorily well with experiments for small axial positions whereas the differences become larger for 30 and 60 mm above the burner inlet. On the other hand, the RSM agrees better with experiments than the k–ω SST for higher axial positions. The Realizable k–ε 2 layer model provides the poorest predictions between the three alternatives. Based on this and the axial velocity comparison the k–ω SST turbulence model was chosen as the best option for the combustion simulations discussed in the next section.

Table 4 Swirl number comparison at different heights
Full size table

3.1.2 Reacting Flow

In a first step, the best combustion model was found by comparison against the TECFLAM experiments. The first assessment of the performance of the candidate models was to compare, the axial velocity (Schneider et al. 2005), temperature and species concentration data (Gregor et al. 2009) with experiments. Due to thermal expansion, the gas velocity and spreading angle of the annular jet are significantly different from the non-reacting case. Figure 6 illustrates the differences between EBU and FGM combustion models. It is clear from the figure that FGM agrees much better with the experimental measurements. However, the magnitude of the reverse flow velocity at the axis is significantly higher with FGM than in the experiments.

Fig. 6
figure 6

FGM, EBU and experimental data axial velocity comparison

Full size image

In Fig. 7 the temperature profile for the EBU model, in contrast to experiments, exhibits two peaks, while FGM agrees much better with the experimental results except close to the axis where FGM predicts an almost constant temperature while the experiments exhibit a slight dip. The explanation for this qualitative as well as quantitative difference between the two models, is the M-shaped flame front in the EBU simulations whereas the FGM model predicts a V-shaped flame front (see Fig. 8). As a consequence of the M-shaped flame front, a line plot of temperature along a horizontal line will cut through the flame front at two points and exhibit two peaks. The fact that the experimental flame front is V-shaped leads to the conclusion that FGM is the more realistic of the two models.

Fig. 7
figure 7

FGM, EBU and experimental data temperature comparison

Full size image
Fig. 8
figure 8

EBU (left) and FGM (right) temperature fields side by side

Full size image

The peak temperature with the FGM model matched excellently with the experiments, when radiation is included (see Fig. 8). On the other hand, if radiation was excluded the simulated peak temperature became around 200 K higher than in the experiment. However, even when radiation was included there was a small mismatch with experiments at 10 mm above the burner where the experiments experience a sharp increase in temperature close to the axis. It is clear from Fig. 7 that the experimental temperature profiles diffuse in the radial direction much faster than the simulations, which is surprising since eddy viscosity models often tend to exaggerate turbulent diffusion compared to experiments (Pope 2000).

In an independent work, Pampaloni et al. (Pampaloni et al. 2017) investigated methane combustion in the TECFLAM geometry, with the inclusion of the swirler geometry. The RANS result of Pampaloni et al. differed to the present results by always having the maximum slope of the temperature profile outside of the experimental maximum, in contrast to the present case that always had its maximum inside the experimental maximum slope. With respect to the steepness of the slope both the results of Pampaloni et al. and the present results are much steeper than the experimental profile. The explanation for the difference in location of the maximum gradient is most likely due to the different inlet boundary conditions but it is obvious that both cases fail to capture the strong diffusive spread of the temperature profile that occurs in the experiment. One significant feature of the results in Pampaloni et al. is that the maximum value of the radial temperature profile was 200–300 K higher than the experiments. This stems from the absence of radiation in the simulations.

A possible explanation for the differences in the slope between the experimental temperature profile and the present case could be that the turbulent Schmidt and Prandtl numbers that are part of the model are improperly set. However, exploratory simulations with a range of values for these dimensionless numbers, and also the Lewis number, did not have a significant impact on the results.

An investigation of the importance of the co-flow boundary conditions was also carried out with a non-uniform co-flow but with the same average velocities (0.5 m/s). This was done by including the upstream geometry for the co-flow according to the description of the experiment in Schneider et al. (2005). However, no significant improvements were observed. The best improvement that could be obtained for diffusion was if the co-flow velocity was quadrupled (2 m/s), which resulted in temperature results that were almost identical to the experimental values but since this is an unrealistic change it will not be discussed further. The most likely explanation to the under-diffusion in the numerical results for the post flame region is that the eddy viscosity RANS formulation is unable to capture the underlying complex phenomena. Hence, modelling with a more advanced LES or possibly an URANS model could potentially provide better agreement due to a better mixing description (see Kuenne et al. 2017). Revisiting Fig. 8, it is evident that the EBU model, bearing in mind all of its shortcomings, is more diffusive than FGM. This trait of FGM has also been reported by Samiran et al. (2019), where inefficient diffusion would lead to underpredicted NOX and in Donini et al. (2015) who conducted a study about including differential diffusion to FGM but concluded the effect was minor.

Looking at the results from Palanti et al. (2018) and Pampaloni et al. (2017) for smaller radial distances than about 20 mm the core temperature is almost constant. The results in the present paper are similar for small axial distances, but for a distance of 60 mm they start to show a similar behaviour to the experiments, with a trough in the profile at the axis, followed by a small increase in temperature with radial position and finally a strong gradient towards lower temperatures at the outer edge of the annular jet. However, the steepness of the gradient in the simulation is still much larger than the experiments.

The depth of the trough in the temperature profile at the axis is most prominent at a height of 10 mm above the burner inlet (Fig. 7). Most likely, the cooling effect from the water-cooled bluff body is the source of this behaviour and the RANS simulations are underestimating this effect.

To assess the numerical accuracy of the simulations, Richardson extrapolation was used to estimate errors (Roache 1994). The variable used for the extrapolation was a line average from the axis to a radial distance of 90 mm for the temperature at 120 mm above the burner inlet. The choice of this value was made based on its sensitivity to mesh changes. Tabulated in Table 5, are the parameters for the error estimation. The original grid used for all the simulations is Case 1 in Table 4 and \(h\) is the typical edge length of a mesh element estimated from the number of grid points as follows:

$$ H = \left[ {\frac{1}{N}\mathop \sum \limits_{{i = 1}}^{N} \varDelta V_{i} } \right]^{{1/3}} , $$
(10)

where \(\varDelta V_{i}\) is the volume of the ith cell and \(N\) is the total number of cells (Celik et al. 2008).The estimated relative error for the finest grid based on this variable is 0.38% and the apparent order of the scheme is approximately 1.87, which agrees well with the formal second order accuracy of the basic discretisation scheme. Judging from the standard deviation of the measured temperatures, which on average is 320 K at 30 mm above the burner inlet, it appears that the simulation results are within an acceptable range for the Case 1 grid.

Table 5 Parameters for error estimation
Full size table

Figure 9 shows a comparison of the CH4 mole fraction with experiments (Gregor et al. 2009). The FGM results agree well with experiments close to the burner while the molar concentration of the fuel is overpredicted further downstream, indicating a slower rate of reaction than in the experiments. It was assumed that this stemmed from the chemical kinetic mechanism that was used; nevertheless, no other mechanisms have been tested in this work. Notice that the resolution in the experiments seems to be relatively low, indicated by the stepwise variation of the molar fraction, and that more accurate experimental data are needed for further development of the model.

Fig. 9
figure 9

Experimental methane mole fraction measurements comparison with FGM simulations

Full size image

3.2 Effect of Laminar Flame Speed

The goal of the paper was to investigate how the combustion behaviour in a swirl stabilised flame is altered when methane is replaced by the syngas from a biomass gasifier. Since FGM with k–ω SST turbulence model provided satisfactory results for methane, it was assumed that the same model would be able to predict the combustion of syngas with the same burner and outer geometry. The assumption behind the model for syngas combustion is that a validated model for methane combustion should be able to predict the general behaviour of a more complex mixture of gases as long as the methane model is based on a chemical scheme in which these species are present.

Due to the large variability of biomass syngas composition (see Table 1) several syngas compositions were investigated to see how much the detailed composition influences the results. The exact gas mixtures that were used in the simulations are listed in Table 3. One hypothesis that was investigated was that the flame shape ought to be quite sensitive to changes in the flame speed but relatively insensitive to the HHV. Hydrogen enriched methane was added in order to have a case with similar laminar flame speed to one of the syngas alternatives but with a significantly higher HHV.

The axial velocities of the different fuels are compared with methane in Fig. 10. Close to the burner (10 and 20 mm) the peak axial velocities were found to be almost identical for all fuels except for BL clean syngas which had a slightly higher velocity at the outer edge of the annular jet. Moreover, for larger radial positions, where reactions are slow due to the low temperature (cf. Fig. 7), the axial velocity was almost identical for all fuels. At higher axial positions (30 and 60 mm) methane experiences a significantly lower axial velocity peak than the rest of the fuels, probably as a result of turbulence chemistry interactions which are tied to the slower rate of heat release with methane as a fuel. The higher axial velocities associated with the syngas fuels (and enriched methane), increase axial momentum and decrease the swirl number, effectively weakening the IRZ. In addition, it was found that the axial velocities of the entrained bark powder syngas and hydrogen enriched methane, which have the same laminar flame speed at the relevant equivalence ratio, are almost identical.

Fig. 10
figure 10

Axial velocity comparison of all the simulated cases

Full size image

In Fig. 11 the streamlines of the flow are shown, together with a fixed isosurface for the OH mole fraction to indicate the flame position, for the different cases. It can be seen that the flame shape, indicated by the OH isosurface is clearly affected by the gas composition and that a higher flame speed seems to correlate well with the flame being bent outwards and the root of the flame moving closer to the burner with increasing flame speed.

Fig. 11
figure 11

Streamlines indicating the position of the recirculation zones. OH isosurfaces indicate the position of flame fronts. The limiting values for the OH isosurfaces (40–80% of the maximum value) in this figure are selected to outline the flame brush

Full size image

The OH mole fraction range that was used to define the flame front was 40–80% of the maximum mole fraction of OH (case specific) to generate a well-defined region. A point to remember here is that the reaction rate that the CFM uses (see Eq. 7) is directly related to the laminar flame speed. Methane with a flame speed of approximately 29 cm/s (at Φ = 0.833) produces a flame front that can be described as a V-shape, with a relatively small angle to the axial direction and with a higher axial position than the syngas flames. The syngas mixtures from Entrained bark powders, BL Raw and 40% hydrogen enriched methane have much higher laminar flame speeds than methane and will as a result be able to propagate significantly longer upstream. These three syngas mixtures have similar flame fronts, which most likely stems from their almost identical laminar flame speeds of 40.62, 45 and 41 cm/s for Entrained, BL Raw and hydrogen enriched methane respectively (see Table 3). These flame fronts were stabilised closer to the inlet, had a larger angle to the axial direction (wider flame) and were also shorter compared to the methane flame. The differences between the syngas cases start to become more visible for the DFB and the BL Clean cases that have higher laminar flame speeds of 51 and 88 cm/s. These flames were bent even further away from the axial direction and had a stronger curvature than the other fuels. The overall effect of an increasing laminar flame speed appears to be a gradual transition from a V-shape towards an M-shape with increasing flame speed. In order to enhance the differences between all cases and to illustrate the effect of increasing flame speed on the flame shape, the flame fronts of methane, hydrogen enriched methane, DFB and BL Clean are plotted together in Fig. 12. The BL Raw flame front lies in between DFB and hydrogen enriched methane, and the Entrained flow syngas flame front lies exactly on the hydrogen enriched methane, as would be expected if the laminar flame speed is the controlling parameter of the flame shape. These two cases are not included in Fig. 12 to avoid cluttering.

Fig. 12
figure 12

Superimposed flame fronts of 4 different cases

Full size image

Figure 13 attempts to characterise further the bending of the flame fronts by showing the OH mole fraction in three axial positions (30, 40 and 50 mm above the burner inlet) for methane, hydrogen enriched methane, DFB and BL Clean. For the methane case, the 30 mm plot passes close to the root of the flame and thus, in comparison with the 40 and 50 mm line plots a peak with a smaller amplitude appears. For hydrogen enriched methane, the radial plots cross the flame front at intermediate heights (neither top nor root, see Fig. 12), hence, all plots have similar peaks. For the DFB case the 50 mm plot passes close to the top boundary of the flame front, showing a lower peak than the 30 and 40 mm plots. The most interesting one is BL Clean for which the 50 mm plot line is outside the flame front as defined in Fig. 13. However, there is a small perturbation of OH at the outer limit of the plot from the wake of the flame. This means that BL Clean produces a shorter flame than the other cases presented in the figure. In addition, the 40 mm plot experiences a less sharp peak than the rest of the cases in that height. This is due to the higher end of the flame front being almost in line with the 40 mm plot. These observations lead to the hypothesis that a higher laminar flame speed, not only should shorten the flame front but also result in bending of the flame brush towards the inlet. Moreover, a higher flame speed seems to yield a more curved flame brush. Fuels with a higher laminar flame speed than BL Clean, would likely continue this trend, with an even shorter flame front that has a larger angle to the axial direction and an even stronger curvature that would make the flame resemble the letter M.

Fig. 13
figure 13

OH mole fraction at 3 different axial positions for the 4 cases of Fig. 13

Full size image

In Fig. 14 the temperature and OH mole fraction fields have been plotted for all the cases. BL Clean can be seen to have a larger high temperature area due to its higher reactivity and flame speed that leads to earlier heat release. While a prediction of the peak temperature in the flow could be based on 1D calculations of the adiabatic flame temperature, the 1D simulations would not give any information about the position and size of the hot zone, which can differ significantly for different fuels, as is shown in Fig. 14. Furthermore, there are two distinct modes of convective cooling of the flame that are impossible to predict with 1D simulations. One stems from the interaction with the surrounding cold co-flow and the other one stems from the interaction with the inner recirculation zone, which contains combustion products that are colder than the ones located immediately after the flame front. Nonetheless, the recirculated products have a high enough temperature to ignite the fresh premix coming from the inlet and maintain the flame.

Fig. 14
figure 14

Temperature and OH mole fraction fields for the different cases. The temperature is plotted as a continuous colour scale while the OH mole fraction is plotted as isolines

Full size image

An interesting observation is the similarity between the temperature fields for 40% hydrogen enriched methane and syngas from Entrained bark powders (see Fig. 15). The results for the temperature at two axial positions are almost identical for both fuels in spite of the large difference in heating value and different reaction pathways. The differences in the combustion behaviour of these two fuels are reflected by the amount of OH present in the flame front (see Fig. 14), which for the hydrogen enriched methane is higher due to the higher amount of H2 present in its composition (see Table 3).

Fig. 15
figure 15

Temperature comparison of hydrogen enriched methane and entrained bark powders at different axial positions

Full size image

4 Conclusions

Premixed combustion of methane in the TECFLAM experiment has been simulated with several turbulence and combustion combinations. The results were compared with experimental data and it was found that the k–ω SST turbulence model and the FGM combustion model combined with the GRI-MECH 3.0 chemical scheme yielded a good agreement with the experiments. It was also concluded that it is necessary to include the effects of radiative heat transfer to be able to predict the peak temperature of the flame. The largest discrepancy between the simulations and the experiments is that the experimental temperature profiles are spreading faster in the radial direction than the simulation, resulting in more smeared out profiles.

The validated model was used to explore the behaviour of biomass syngas in the same burner. The first step in the exploratory study was to compare the flow field from simulations with 4 typical syngas mixtures, 40% hydrogen enriched methane and pure methane. A qualitative difference between the methane and the rest of the gas mixtures was found to be that the magnitude of the reverse flow was significantly higher for the syngas cases than for methane. The syngas cases were found to have increased axial momentum and reducing swirl intensity, which weakens the IRZ. The outer recirculation zone was unaffected by the changes in fuel composition.

The simulation results based on the CFM seemed to support the hypothesis that the laminar flame speed of the premix is the dominant parameter that controls the flame shape. This hypothesis was tested further by comparing the flame shape from the two distinctly different gas mixtures that had the same laminar flame speed but a large difference in heating value and composition, and the outcome gave further support of the hypothesis. Overall, based on the current model, it can be concluded that a higher laminar flame speed will result in shorter flames, a bended flame front away from the axial direction and an increased curvature of the flame front.