A lattice Boltzmann model for reactive mixtures.

A new lattice Boltzmann model for reactive ideal gas mixtures is presented. The model is an extension to reactive flows of the recently proposed multi-component lattice Boltzmann model for compressible ideal gas mixtures with Stefan-Maxwell diffusion for species interaction. First, the kinetic model for the Stefan-Maxwell diffusion is enhanced to accommodate a source term accounting for the change in the mixture composition due to chemical reaction. Second, by including the heat of formation in the energy equation, the thermodynamic consistency of the underlying compressible lattice Boltzmann model for momentum and energy allows a realization of the energy and temperature change due to chemical reactions. This obviates the need for ad-hoc modelling with source terms for temperature or heat. Both parts remain consistently coupled through mixture composition, momentum, pressure, energy and enthalpy. The proposed model uses the standard three-dimensional lattices and is validated with a set of benchmarks including laminar burning speed in the hydrogen-air mixture and circular expanding premixed flame. This article is part of the theme issue 'Progress in mesoscale methods for fluid dynamics simulation'.

Introduction

The lattice Boltzmann method (LBM) is a recast of fluid dynamics into a fully discrete kinetic system for the populations of designer particles, which are associated with the discrete velocities fitting into a regular space-filling lattice. As a result, the kinetic equations for the populations follow a simple algorithm of ‘stream along links and collide at the nodes in discrete time t’. LBM has been successfully applied to a range of problems in fluid dynamics including but not limited to transitional flows, flows in complex moving geometries compressible flows, multiphase flows and rarefied gas, to name a few [1,2].

Nevertheless, in spite of extensive development, the multicomponent reactive mixtures so far resisted a significant advancement in the LBM context. Arguably, one of the main reasons was the absence of a thermodynamically consistent LBM for mixtures. Early approaches such as [3,4] suffer many limitations such as incompressible flow restriction, constant transport properties and rudimentary diffusion modelling.

As a remedy, a number of recent works [5-10] abandoned the construction of a kinetic model or LBM for multicomponent mixtures in favour of a so-called hybrid LBM where only the flow of the mixture is represented by an (augmented) LBM equation while the species and the temperature dynamics are modelled by conventional macroscopic equations. While the hybrid LBM approach can be potentially useful, in particular for combustion applications, our goal here is to retain a fully kinetic model and LBM for multicomponent reactive mixtures.

Recently, we proposed a novel lattice Boltzmann framework for compressible multi-component mixtures with a realistic equation of state and thermodynamic consistency [11]. The strongly coupled formulation consists of kinetic equations for momentum, energy and species dynamics and was validated for a variety of test cases involving uphill diffusion, opposed jets and Kelvin–Helmholtz instability. This extends the LBM to realistic mixtures and opens the door for reactive flow applications with a fully kinetic approach, which is the subject of this paper. We propose a fully kinetic, strongly coupled lattice Boltzmann model for compressible reactive flows as an extension of [11]. To that end, a generic M-component ideal gas mixture is represented by two sets of kinetic equations. A set of M kinetic equations is used to model species undergoing Stefan–Maxwell diffusion. This set is now extended to include the reaction source term. Furthermore, the mixture is described by a set of two kinetic equations, where one accounts for the total mass and momentum of the mixture and another for the total energy of the mixture. The kinetic equation for the mixture energy is extended to also include the internal energy of formation in addition to the sensible internal energy. Thus, the approach presented here can accurately model a reactive M-component compressible mixture with M+2 kinetic equations. The system is fully coupled through mixture composition, momentum, pressure and enthalpy. The thermodynamic consistency of the model allows us to automatically account for the energy changes due to chemical reactions. The Stefan–Maxwell diffusion is retained and thus complicated phenomena such as reverse diffusion, osmotic diffusion or diffusion barrier can be captured, as was already demonstrated in the non-reactive case in [11].

The outline of the paper is as follows. In §2, we extend the lattice Boltzmann model of [11] to the reactive multicomponent mixtures. This is achieved by supplying a reaction source term to the kinetic equations for the species in such a way that the Stefan–Maxwell diffusion mechanism already implemented by the model stays intact. In §3, we extend the two-population lattice Boltzmann model for the mixture flow and energy to include the enthalpy of formation of chemically reacting species. Thanks to the thermodynamic consistency featured by the original model [11], this final step completes the construction of the lattice Boltzmann model for the reactive mixtures. The derivation follows the path presented in detail in [11], and we indicate the differences brought about by the thermodynamics of the chemical reaction. In §4, we outline the coupling of the lattice Boltzmann solver with the open source chemical kinetics package Cantera. Validation of the model is presented in §5 with the simulation of detailed hydrogen/air combustion mechanism and the discussion is provided in §6.

Lattice Boltzmann model for the species

The composition of a reactive mixture of M ideal gases is described by the species densities , , while the mixture density is . The rate of change of due to chemical reaction satisfies mass conservation, Introducing the mass fraction , the molar mass of the mixture is given by where is the molar mass of the component a. The equation of state of the mixture provides a relation between the pressure P, the temperature T and the composition, where is the specific gas constant of the mixture and is the universal gas constant. The pressure of an individual component is related to the pressure of the mixture P through Dalton’s Law of partial pressures, , where the mole fraction of a component is related to its mass fraction as . Combined with the equation of state (2.2), the partial pressure takes the form , where is the specific gas constant of the component. A kinetic model for the Stefan–Maxwell diffusion in the non-reactive mixture was introduced in [11]. Here, we extend the formulation [11] to include the reaction. To that end, we write the kinetic equation for the populations , , of the component a, corresponding to the discrete velocities , , Here are the binary diffusivity coefficients. The species’ densities and partial momenta are, respectively, The momenta of the components sum up to the mixture momentum, At variance with the non-reactive mixture [11], kinetic equation (2.3) includes a source term which implements the rate of change of due to the reaction and satisfies the following conditions, The kinetic model (2.3) is realized on the standard three-dimensional lattice with the discrete velocities . As in [11], the equilibrium and the quasi-equilibrium in (2.3) are constructed using the product-form [12]: we define a triplet of functions in two variables, and , The equilibrium and the quasi-equilibrium populations are evaluated as the products of the functions (2.6), with and , respectively, and with in both cases, and The reaction source term in (2.3) is also represented by the product-form similar to (2.7), The analysis of the hydrodynamic limit of the kinetic model (2.3) follows the lines already presented in [11]. Note that the constraint on the momentum of the source term (2.5) is required. The balance equations for the densities of the species in the presence of the source term are found as follows: where the diffusion velocities, , satisfy the Stefan–Maxwell constitutive relation, Summarizing, kinetic model (2.3) recovers both the Stefan–Maxwell law of diffusion and the contribution of the species mass change due to chemical reaction, as presented in equation (2.10).

Derivation of the lattice Boltzmann equation from the kinetic model (2.3) proceeds along the lines of the non-reactive case [11]. Upon integration of (2.3) along the characteristics and application of the trapezoidal rule, we arrive at a fully discrete lattice Boltzmann equation The discrete equation (2.12) is used in the actual numerical implementation. Unlike equation (2.3), the discrete equation does not contain species mass in the denominator. Therefore, no special treatment is required for species masses going to zero. The short-hand notation for the inter-species interaction term and the relaxation parameters are, where the characteristic times and the relaxation times are related to the binary diffusivities, Furthermore, the quasi-equilibrium populations in the expression (2.13) depend on the diffusion velocity . The latter are found by solving the linear algebraic system for each spatial component The linear algebraic system was already derived in [11] for the non-reactive mixtures and is not altered by the presence of the reaction source term. The equilibrium population and the reaction source term in (2.12) and (2.13) are evaluated at the mixture velocity . Summarizing, the lattice Boltzmann system (2.12) delivers the extension of the species dynamics subject to the Stefan–Maxwell diffusion to the reactive mixtures. We proceed with the extension of the flow and energy dynamics of the mixture.

Lattice Boltzmann model of mixture momentum and energy

The mass-based specific internal energy and enthalpy of a species a are,

where and are, respectively, the energy and the enthalpy of formation at the reference temperature , while and are specific heats at constant volume and at constant pressure, respectively. The internal energy and the enthalpy of a mixture are While the sensible heat was considered in the non-reactive case [11], by taking into account the heat of formation we immediately extend the model to reactive mixtures. As in [11], we follow a two-population approach. One set of populations (f-populations) is used to represent the density and the momentum of the mixture Another set (g-populations) represents the total energy, A coupling between the mixture and the species kinetic equations is established through energy since the mixture internal energy (3.2) depends on the composition. Furthermore, the temperature is evaluated by solving the integral equation, cf. (3.1) and (3.2), The temperature is used as the input for the equation of state (2.2) and hence in the equilibrium, the quasi-equilibrium and the reaction source term of the species lattice Boltzmann system which leads to a two-way coupling between the species and the mixture kinetic systems. As in [11], the lattice Boltzmann equations for the f- and g-populations are realized on the discrete velocity set and where relaxation parameters and are related to the viscosity and thermal conductivity. The equilibrium f-populations in (3.6) are evaluated using the product-form, with and in (2.6), The last term in (3.6) is a correction needed to compensate for the insufficient isotropy of the lattice in the compressible flow setting [11,13-15]: is the vector with the components while the components of vectors are defined as The equilibrium and the quasi-equilibrium g-populations, and in (3.7), are defined with the help of Grad’s approximation [16], and Here, the weights are the products of the one-dimensional weights , , and is the lattice reference temperature. The equilibrium mixture energy flux and the second-order moment tensor in (3.11) and (3.12) are and where is the specific mixture enthalpy (3.2). The quasi-equilibrium energy flux in (3.12) has the following form: The first two terms in (3.15) include the energy flux and the pressure tensor , Their contribution maintains a variable Prandtl number and is patterned from the single-component case [13]. The remaining two terms in the quasi-equilibrium energy flux (3.15), and pertain to the multicomponent case. The interdiffusion energy flux is where the diffusion velocities are defined according to equation (2.15). The flux (3.17) contributes the enthalpy transport due to diffusion and hence it vanishes in the single-component case but is significant in reactive flows. Finally, the correction flux , which also vanishes in the single-component case, is required in the two-population approach to the mixtures in order to recover the Fourier law of thermal conduction, see [11] for details, Prefactors featured in (3.17) and (3.18) were found in [11] based on the analysis of the hydrodynamic limit of the lattice Boltzmann system (3.6) and (3.7) and are not affected by the present reactive mixture case. Second-order accurate isotropic lattice operators proposed in [17] were used for the evaluation of spatial derivatives in the correction flux (3.18) as well as in the isotropy correction (3.9). Following [11], the continuity, the momentum and the energy equations for a reactive multicomponent mixture [18] are obtained as follows:

The pressure tensor in the momentum equation (3.20) reads where the dynamic viscosity and the bulk viscosity are related to the relaxation parameter , where is the mixture specific heat at constant volume. The heat flux in the energy equation (3.21) reads The first term is the Fourier law of thermal conduction, with the thermal conductivity related to the relaxation parameter , where is the mixture specific heat at constant pressure. The second term in (3.24) is the interdiffusion energy flux. The dynamic viscosity and the thermal conductivity of the mixture are evaluated as a function of the local composition, temperature and pressure using the chemical kinetics solver Cantera [19], wherein a combination of methods involving interaction potential energy functions [20], hard sphere approximations and the methods described in [21,22] are employed to calculate the mixture transport coefficients. Finally, in accordance with a principle of strong coupling [11], the excess conservation laws arising due to a separated construction of the species diffusion model in §2 and the two-population mixture model are eliminated by removing one set of species populations (here, the component M), Thus, the component M is not an independent field any more but is slaved to the remaining species and the mixture f-populations. Summarizing, the thermodynamically consistent framework of [11] allows for a straightforward extension to reactive mixtures provided the sensible energy and enthalpy are extended to include the energy and the enthalpy of formation.

Coupling between lattice Boltzmann and chemical kinetics

In this work, the lattice Boltzmann code is coupled to the open source code chemical kinetics solver Cantera [19]. The Cantera solver is supplied with the publicly accessible GRI-Mech 3.0 mechanism [23] as an input data file. The communication between the lattice Boltzmann solver and the Cantera chemical kinetics solver can be summarized as follows:

An input from the lattice Boltzmann solver to Cantera is provided during the collision step in terms of internal energy, specific volume and mass fractions.

Cantera internally solves numerically the integral equation (3.5) and thus the temperature at that state is obtained.

Cantera calculates the production rates of species and the transport coefficients including dynamic viscosity, thermal conductivity and the Stefan–Maxwell diffusivities as a function of the current state.

The temperature obtained from Cantera is used to evaluate the equilibrium and quasi-equilibrium moments and populations. The transport coefficients are used to calculate the corresponding relaxation times and thus the collision step is complete.

Other thermodynamic parameters necessary for the simulations such as the specific heats and molecular masses are also obtained through Cantera. The reference standard state temperature is and the reference standard state pressure is . The data required by the lattice Boltzmann solver during runtime is obtained by querying Cantera through its C++ API using the ‘IdealGasMix’ and ‘Transport’ classes. A flowchart of the coupling between the lattice Boltzmann solver and the chemical kinetics solver Cantera is sketched in figure 1. The diffusivities obtained from Cantera are converted to lattice units through the non-dimensional Prandtl number and the Lewis number . Here, is the thermal diffusivity of the mixture.

Figure 1

Flowchart for the coupled lattice Boltzmann–Cantera solver. (Online version in colour.)

Results

To test the coupling between the lattice Boltzmann solver and the chemical kinetics solver, we begin by verifying the perfectly stirred reactor set-up with homogeneous hydrogen/air mixture. Next, as a second validation, probing the basic validity of our model, we compute the flame speed in a premixed hydrogen/air mixture with the reactive Stefan–Maxwell formulation in a wide range of equivalence ratios . Subsequently, in order to test the isotropy of the model, the problem of outward expanding circular flame [24,25] is solved for the premixed hydrogen/air mixture. For both test cases, we use the detailed chemical kinetics mechanism [26] involving the following nine species: , , , , , , , , . It is worthwhile to mention that the model is not restricted to the detailed mechanisms. Reduced mechanisms available in the literature such as the five-species propane mechanism have also been tested with this model. In this paper, we will restrict ourselves to the more interesting detailed hydrogen/air mechanism which forms sharper and faster propagating flames. While this benchmark not only probes the model’s behaviour in two dimensions, it is also a stringent isotropy test where it is crucial that the circular shape of the flame is preserved and not contaminated or distorted by the errors of the discrete numerics on the underlying Cartesian grid. Finally, the model’s ability to capture nonlinear instabilities is probed by simulations of wrinkled flames, which form as a result of monochromatic perturbations.

Perfectly stirred reactor

To verify the coupling between the lattice Boltzmann solver and the chemical kinetics solver, we simulate a constant volume reactor. The set-up consists of a three-dimensional domain of lattice points with periodic boundary conditions in all directions. Stagnant hydrogen/air mixture at an equivalence ratio is initialized with the pressure and the temperature . The results are compared to the ‘IdealGasReactor’ class of Cantera using the ‘advance’ method to march the state of the reactor forward in time. Figure 2a shows the history of the temperature and the mass fraction of hydroxide OH. Accurate match with the results obtained from Cantera shows that the coupling has been done correctly. For this problem, since all the boundaries are periodic, the total energy of the system has to remain constant while no kinetic energy should develop over time. Figure 2b verifies that in the absence of kinetic energy, the total energy of the solver equals the kinetic energy. Also, the total energy does not deviate from its initial value in this set-up, as expected.

Figure 2

Simulation of hydrogen/air constant volume perfectly stirred reactor. (a) Time evolution of temperature and OH mass fraction, (b) time evolution of the kinetic, internal and total energy. All quantities are scaled by the initial total energy . (Online version in colour.)

Laminar flame speed

In order to validate our model, we calculate the burning velocity of a hydrogen/air mixture in a one-dimensional set-up. As illustrated in figure 3, the set-up consists of a one-dimensional tube initialized with unburnt mixture at throughout from the left end up to 80% of the domain towards the right. The remaining 20% of the domain is initialized with the adiabatic flame temperature and with the equilibrium burnt composition at the respective equivalence ratio. The pressure is initialized uniformly at . Zero gradient boundary conditions are used at both ends for all variables using equilibrium populations. At the left end, the velocity is imposed to be zero so that the flame propagates from right to left into the stationary unburnt mixture. The set-up is used to calculate the burning velocity of the premixed system. Nitrogen is considered as an inert gas and thus does not split or form any radicals like nitrous oxides. However, the heat capacity of the inert gas has a strong influence on the flame temperature and consequently on the burning velocity. This is naturally accounted for in the formulation. The burning velocity is measured for various equivalence ratios ranging from to . We use the laminar flame thickness at for defining the reference length, where . In order to accurately calculate the burning velocity, we use a long domain of , which corresponds to lattice points. In order to avoid the effect of the boundaries and transients due to initial acceleration, the flame speed is measured when the flame front approaches the middle of the domain. The results are compared to the data provided by [27] from multiple experimental and computational sources in figure 4. It can be seen that flame speed computed by our model agrees well with the data available in the literature. Although there is considerable dispersion in the literature for the flame speed values for fuel-rich mixtures , the location of the peak burning velocity between and has been correctly captured. As is evident in figure 5, the profiles of temperature and mass fractions for show a good match compared to the corresponding solution obtained from the ‘FreeFlame’ class of Cantera. This test case indicates that the present model is a promising candidate for simulating reactive flows with the lattice Boltzmann method.

Figure 3

Set-up for the burning velocity simulation. (Online version in colour.)

Figure 4

Burning velocity versus equivalence ratio for the nine-species hydrogen/air mixture detailed chemistry [27]. (Online version in colour.)

Figure 5

Profiles of temperature and mass fractions for one-dimensional planar flame at . (Online version in colour.)

Circular expanding premixed flame

After confirming the one-dimensional behaviour of the model, we compute the two-dimensional circular expanding flame in a premixed hydrogen/air mixture with detailed chemistry. Similarly to the study [24,25], due to symmetry, only a quarter of the flame is solved. Symmetry boundary conditions are used on the left and bottom edges of the square domain while the characteristic based outlet boundary conditions [28,29] are imposed at the right and top edges of the domain. The bottom left corner is initialized with a burnt quarter sector at the adiabatic flame temperature corresponding to the equivalence ratio . The rest of the domain is initialized with an unburnt mixture at the temperature . The composition in the burnt section is set to the equilibrium composition and the pressure in the entire domain is initialized to a uniform pressure . For this premixed initial condition, the burning velocity is obtained as from solving a one-dimensional flame propagation set-up in Cantera. The flame thickness at these initial conditions is obtained as . A square domain with the side was considered, which corresponds to lattice points. The radius of the region initialized with the burnt equilibrium conditions is .

The characteristic flame transit time is defined as [25]. Contours of temperature, velocity and mole fractions of oxygen and the hydroxide radical are shown at in figure 6b. As can be verified from figure 6b, the solution is not contaminated by numerical noise or anisotropies and the contours do not contain any other spurious features. The thin interface of the hydroxide radical at the flame front is captured correctly and the curvature of the flame is maintained. This is in contrast to e.g. [25], where the errors of the underlying numerical discretization leading to a spurious behaviour were reported when using Cartesian grids.

Figure 6

Premixed hydrogen/air circular outward expanding flame. (a) Contour plot of the mole fraction of at obtained by reflecting about the left edge and the bottom edge of the domain, (b) Contours of temperature, mole fractions of , OH and velocity at . (Online version in colour.)

Next, we study the response of this set-up to a deterministic perturbation to validate the model with the Direct Numerical Simulation (DNS) of [25]. The initial circular profile of the flame is perturbed with a sinusoidal profile, where corresponds to the number of modes of the perturbation per sector of the flame and is the amplitude of the perturbation. The evolution of the perturbation is shown in figure 7. The heat release rate, , is a measure of the reactivity of the mixture. As it is evident in figure 7a, during the initial stages of the evolution, the perturbed modes are continuous and the heat release rate is uniform along the circumference of the flame. As explained in [25], the reactivity and therefore the heat release rate reduces at the crest due to diffusion and more consumption of the deficient reactant. This, along with the hydrodynamic instability due to the density ratio and the thermal-diffusive instability due to the heat and mass imbalance of the deficient reactant, leads to splitting of the peak of the crests into smaller cells, as is visible in figure 7b. A snapshot of the temperature contours over time shown in figure 7c verifies that the splitting of the flame indeed occurs from crests. Therefore, the splitting stems from the deterministic perturbation as expected, and not because of numerical noise. The mean radius of the flame is calculated by integrating along the flame front circumference, Here A is the circumferential length and R is the distance of the mean temperature isoline from the centre. On fitting , the growth exponent was found to be , in agreement with the results from DNS in the literature wherein the value of the exponent was found to be between almost linear [25] and 1.25 [30]. The local displacement speed [24,25] is calculated as With the local flame normal , the absolute propagation speed is calculated as . The density weighted displacement speed is defined as , where is the density of the unburnt mixture. The flame speeds are calculated as a mean over the flame interface isoline of in a way similar to equation (5.2). After the initial transients, the absolute propagation speed was found to reach a value of , whereas the density weighted displacement speed was found to fluctuate about . The corresponding values from the DNS results [25] are and , respectively. The difference could be attributed to a number of factors including the type of grid, resolution, type of diffusion model, etc. Overall, the results agree well with the DNS [24,25].

Figure 7

Contours of temperature and heat release rate. (a) , (b) , (c) line contours of form to . The domain has been reflected about the left and the bottom edge for plotting. (Online version in colour.)

Conclusion

In this paper, we proposed a lattice Boltzmann framework to simulate reactive mixtures. The novelty of the model lies in the fact that temperature and energy changes due to chemical reaction are handled naturally without the need of additional ad-hoc modelling of the heat of reaction. This was possible because of the thermodynamic consistency of the underlying multi-component model [11], which was extended to compressible reactive mixtures. The species interaction is modelled through the Stefan–Maxwell diffusion mechanism which has been extended in this work to accommodate the creation and destruction rates of the species due to chemical reaction. Computational efficiency has been achieved through reduced description of energy which makes it possible to describe the physical system by only M+2 kinetic equations instead of 2M kinetic equations while retaining necessary physics such as the inter-diffusion energy flux. The model has been realized on the standard D3Q27 lattice, which not only reduces the computational costs compared to multispeed approaches but also possesses a wide temperature range, which is crucial for combustion applications.

The proposed model was validated in one and two dimensions with the 9-species 21 steps detailed hydrogen–air reaction mechanism. The accuracy of the model was assessed by calculating the burning velocity of a premixed hydrogen-air mixture in one dimension. The calculated flame speed agrees well with the results in the literature. The ability of the model to capture complex physics was tested by simulating a two-dimensional expanding circular flame. The circular flame simulation exhibited good isotropy and low numerical noise. The set-up was then subjected to monochromatic perturbations in order to study the evolution of the perturbed flame. Good agreement with DNS simulations demonstrates viability of the proposed LBM for complex reactive flows.