Consistent lattice Boltzmann model for reactive mixtures

A new lattice Boltzmann model (LBM) for chemically reactive mixtures is presented. The approach capitalizes on the recently introduced thermodynamically consistent LBM for multicomponent mixtures of ideal gases. Similar to the non-reactive case, the present LBM features Stefan--Maxwell diffusion of chemical species and a fully on-lattice mean-field realization of the momentum and energy of the flow. Besides introducing the reaction mechanism into the kinetic equations for the species, the proposed LBM also features a new realization of the compressible flow by using a concept of extended equilibrium on a standard lattice in three dimensions. The full thermodynamic consistency of the original non-reactive multicomponent LBM enables to extend the temperature dynamics to the reactive mixtures by merely including the enthalpy of formation in addition to the previously considered sensible energy. Furthermore, we describe in detail the boundary conditions to be used for reactive flows of practical interest. The model is validated against a direct numerical simulation of various burning regimes of a hydrogen/air mixture in a microchannel, in two and three dimensions. Excellent comparison in these demanding benchmarks indicates that the proposed LBM can be a valuable and universal model for complex reactive flows.


Introduction
The lattice Boltzmann method (LBM) models fluid flow using a fully discrete kinetic system of designer particles with discrete velocities c i , fitting into a regular space-filling lattice. In the LBM, the kinetic evolution equation for the populations f i (x, t) follows a simple algorithm of "stream along links c i and collide at the nodes x in discrete time t". Since its inception (Higuera & Jiménez 1989;Higuera & Succi 1989), LBM has evolved into a versatile tool for the simulation of complex flows including but not limited to turbulent flows (Dorschner et al. 2016, compressible flows (Frapolli et al. 2016;Dorschner et al. 2018;Xu & Sagaut 2013;Yang et al. 2018;Lin & Luo 2018), multiphase flows (Mazloomi et al. 2015(Mazloomi et al. , 2017Wöhrwag et al. 2018), rarefied gas (Shan et al. 2006) and nanoflow (Montessori et al. 2016;Montemore et al. 2017). While the majority of the LBM development concerns single-component fluids, the case of mixtures, and especially of reactive mixtures, remains an active area of research (Yan et al. 2013;Lin & Luo 2018;Hosseini et al. 2018Hosseini et al. , 2019Hosseini et al. , 2020Feng et al. 2018;Tayyab et al. 2020Tayyab et al. , 2021. Recently, in Sawant et al. (2021a), we revisited the LBM construction for a compressible multicomponent mixture, focusing on a thermodynamically consistent coupling between diffusion and momentum and energy transfer. The species kinetic equations recovered the Stefan-Maxwell diffusion with barodiffusion in the hydrodynamic limit. In addition, we also validated and derived from our kinetic model approximate diffusion models such as Curtiss-Hirschfelder and generalized Fick (Kee et al. 2003;Poinsot & Veynante 2005;Giovangigli 2012). A mean-field was introduced for the lattice Boltzmann formulation of the mixture momentum and energy using a two-population lattice Boltzmann equation for the mixture. The mean-field approach consists of two lattice Boltzmann equations, one for the mixture density and momentum and another one for the energy with the help of a modification of the non-equilibrium fluxes. The two-population mixture LBM and the lattice Boltzmann scheme for the species kinetic equations were realized on the standard three-dimensional lattice. The resulting LBM provides a reduced description of the M -component mixture with M + 2 tightly coupled lattice Boltzmann equations.
In this paper, we extend the two-population mixture LBM to reactive flows and show viability and accuracy for practical applications. To that end, we propose novel boundary conditions for walls as well as inlets and outlets. Furthermore, unlike previous realizations Sawant et al. (2021a,b), we use the extended lattice Boltzmann method (Saadat et al. 2021) for the mean field model. For validation, we start with simulations of a perfectly stirred reactor and a one dimensional laminar flame. Subsequently, combustion of a lean hydrogen/air mixture in microtubes is simulated in two and three dimensions. A variety of different flame dynamics such as the periodic ignition-extinction, stable symmetric V-shaped flames and asymmetric flames are captured by the reactive LBM and our simulations are in quantitative agreement with the direct numerical simulations (DNS) of Pizza et al. (2008bPizza et al. ( , 2010. This demonstrates that the proposed is a viable alternative for the simulation of reactive flows.
The paper is structured as follows: We begin with a recap of the nomenclature and the kinetic system for the species in section 2. This section presents the discrete lattice Boltzmann equations for the reactive species and their implementation on the standard lattice. The section closes with a short discussion on time integration of the reaction mass source term. Next, we turn our attention to describing the mean field approach for modeling the momentum and energy of the reactive mixture in section 3. Here, we discuss the adoption of the extended lattice Boltzmann method for reactive flows and present the realization on standard lattice with the two-population approach. The section closes with a brief outline of the resultant macroscopic Navier-Stokes equations in the continuum limit followed by validation with the perfectly stirred reactor and one dimensional laminar flame. Having completely described the dynamics in bulk of the multicomponent fluid, we proceed to formulate the boundary conditions for the combined model in section 4. In section 4.2, we discuss the equivalent of no-slip adiabatic boundary condition using the popular bounce-back method. This is followed by a technique to implement an isothermal wall based on the Tamm-Mott-Smith boundary condition in the lattice Boltzmann framework. Next, a realization for applying the inlet flux boundary condition is discussed in section 4.3. In the section 4.4, a convective boundary condition, which can be used in conjunction with the characteristics based boundary condition, is provided to approximate the species mass fractions at the outlet. With the resultant model, we compute the combustion of a premixed hydrogen/air mixture flowing through hot microtubes in section 5. Validation is performed in different regimes by changing the inlet velocity. The 2D simulation exhibits rich flame dynamics by undergoing repetitive extinction-ignition, forming stable V-shaped flame and stable asymmetric flames. Finally, a 3D open flame is computed in a microtube.

Kinetic equations for the species
The nomenclature follows Sawant et al. (2021a). The composition of a reactive mixture of M components is described by the species densities ρ a , a = 1, . . . , M , while the mixture density is, (2.1) The rate of change of species densities due to reaction,ρ r a , satisfies mass conservation, Introducing the mass fraction, Y a = ρ a /ρ, the molar mass of the mixture m is given by where m a is the molar mass of the component a. The ideal gas equation of state (EoS) provides a relation between the pressure P , the temperature T and the composition, where R = R U /m is the specific gas constant of the mixture and R U is the universal gas constant. The pressure of an individual component is related to the pressure of the mixture through Dalton's law of partial pressures, P a = X a P , where the mole fraction is X a = mY a /m a . Combined with the equation of state (2.3), the partial pressure takes the form P a = ρ a R a T , where R a = R U /m a is the specific gas constant of the component. In the kinetic representation, each component is described by a set of populations f ai corresponding to the discrete velocities c i , i = 0, . . . , Q − 1. The species densities ρ a and the partial momenta ρ a u a are defined accordingly, while partial momenta sum up to the mixture momentum, (2.6) Following Sawant et al. (2021a,b), the kinetic equations for the species can be written as, where D ab are Stefan-Maxwell binary diffusion coefficients, while the reaction source term satisfies the following conditions, consistent with (2.4): (2.9) We now proceed with specifying the equilibrium f eq ai , the quasi-equilibrium f * ai and the reaction source termḟ r ai .

Standard lattice and product-form
Kinetic model (2.7) is realized on the standard discrete velocity set D3Q27, where D = 3 stands for three dimensions and Q = 27 is the number of discrete velocities, (2.10) In order to specify the equilibrium f eq ai , the quasi-equilibrium f * ai and the reaction source termḟ r ai in (2.7), we first define a triplet of functions in two variables, ξ α and P αα , Ψ 0 (ξ α , P αα ) = 1 − P αα , (2.11) and consider a product-form associated with the discrete velocities c i (2.10), (2.14) All pertinent populations to be encountered in this paper shall be determined by specifying the parameters ξ α and P αα in the product-form (2.14). To that end, the equilibrium and the quasi-equilibrium populations are found by setting, (2.16) in the former, and (2.18) in the latter cases: (2.20) Reaction terms are specified with the product-form (2.14) using the equilibrium parameters (2.16), (2.21) Analysis of the hydrodynamic limit of the kinetic model (2.7) follows the lines already presented in Sawant et al. (2021a). The balance equations for the densities of the species in the presence of the source term are found as follows, where the diffusion velocities, δu a = u a − u, satisfy the Stefan-Maxwell constitutive relation, (2.23) Summarizing, kinetic model (2.7) recovers both the Stefan-Maxwell law of diffusion and the composition change due to chemical reaction, as presented in equation (2.22).

Lattice Boltzmann equation for the species
Derivation of the lattice Boltzmann equation from the kinetic model (2.7) proceeds along the lines of the non-reactive case already presented in detail by Sawant et al. (2021a). Upon integration of (2.7) along the characteristics and application of the trapezoidal rule to all relaxation terms on the right hand side except for the reaction term, we arrive at a fully discrete lattice Boltzmann equation for the species, (2.24) Here δt is the lattice time step, the equilibrium populations are provided by Eq. (2.19), while the relaxation parameters β a ∈ [0, 1] are, β a = δt 2τ a + δt . (2.25) Their relation to the Stefan-Maxwell binary diffusion coefficients is found as follows: Introducing characteristic times, the relaxation times τ a in (2.25) are defined through mixture-averaging, (2.27) Furthermore in (2.24), the quasi-equilibrium relaxation term F ai is spelled out as follows, (2.28) Here the quasi-equilibrium populations f * bi are defined by the product-form (2.20), subject to the following parameterization, (2.29) where the second-order accurate diffusion velocity V b is the result of the lattice Boltzmann discretization of the kinetic equation and is found by solving the M × M linear algebraic system for each spatial component, (2.31) The system (2.31) has been derived in Sawant et al. (2021a) and is not altered by the presence of the reaction. In our realization, we solve (2.31) with the Householder QR decomposition method from the Eigen library (Guennebaud et al. 2010). All the elements of the lattice Boltzmann equation (2.24) described so far are identical to those already present in the non-reactive case of Sawant et al. (2021a). Finally, the reaction term in (2.24) is represented by an integral over the characteristics, R r ai = δt 1 0ḟ r ai (x + c i sδt, t + sδt)ds. (2.32) Taking into account the structure of the reaction term (2.21), we use a simple explicit approximation for the implicit term (2.32), Reaction ratesρ r a are obtained from the open source chemical kinetics package Cantera (Goodwin et al. 2018) as a function of mixture internal energy U and composition, ρ r a =ρ r a (U, ρ 1 , . . . , ρ M ). In order to mitigate the stiffness of the reaction rates for detailed reaction mechanisms, we introduce a time step δt r = δt/l, where l = 1, 2, ... and evaluate (2.33) by forward Euler in l sub-steps, Note that, during sub-iterations, the energy remains fixed although the temperature changes, in general. In other words, at each grid point, sub-iterations implement a zerodimensional perfectly stirred reactor. Execution time for sub-steps increases by about 6% for l = 2 and by 15% for l = 4. In this paper, we use l = 2, which is small enough that the integration error does not influence the flow solution but still reduces the computational complexity by roughly half due to the larger time step of the fluid solver, δt = 2δt r .
Summarizing, the lattice Boltzmann system (2.24) delivers the extension of the species dynamics subject to the Stefan-Maxwell diffusion to the reactive mixtures. We now proceed with setting up the lattice Boltzmann equations for the mixture momentum and energy.
3. Lattice Boltzmann model of mixture momentum and energy 3.1. Double-population lattice Boltzmann equation The mass-based specific internal energy U a and enthalpy H a of the species are, where U 0 a and H 0 a are the energy and the enthalpy of formation at the reference temperature T 0 , respectively, while C a,v and C a,p are specific heats at constant volume and at constant pressure, satisfying the Mayer relation, C a,p − C a,v = R a . Consequently, the internal energy ρU and enthalpy ρH of the mixture are defined as, While the sensible heat was considered in the non-reactive case (Sawant et al. 2021a), by taking into account the heat of formation we immediately extend the model to reactive mixtures. Same as in Sawant et al. (2021a), we follow a two-population approach. One set of populations (f -populations) is used to represent the density and the momentum of the mixture. Below, we refer to the f -populations as the momentum lattice. The locally conserved fields are the density and the momentum of the mixture, Another set of populations (g-populations), or the energy lattice, is used to represent the local conservation of the total energy of the mixture, (3.8) Since the mixture internal energy (3.3) depends on the composition, the species kinetic equations become coupled with the kinetic equations for the mixture to be introduced shortly. Conversely, the temperature is evaluated by solving the integral equation, cf.
(3.1) and (3. (3.9) The temperature evaluated by solving (3.9) is used as the input in the equation of state (2.3) elsewhere in the species lattice Boltzmann system. This furnishes a two-way coupling input between the species and the mixture kinetic systems. Similar to Sawant et al. (2021a), the lattice Boltzmann equations for the momentum and for the energy lattice are patterned from the single-component developments and are realized on the D3Q27 discrete velocity set. While the prototype single-component LBM used in Sawant et al. (2021a) was that of Saadat et al. (2019), here we take advantage of a more recent proposal by Saadat et al. (2021). It is noted that, while both these singlecomponent models are essentially equivalent, the recent formulation is more compact in its formulation and simpler in terms of implementation. Following the more recent proposal, the mixture lattice Boltzmann equations are written, where relaxation parameters ω and ω 1 are related to the mixture viscosity and thermal conductivity, and we proceed with specifying the pertinent populations in (3.10) and (3.11).

Extended equilibrium for the momentum lattice
The extended equilibrium populations f ex i in (3.10) are specified by the product-form (2.14), with the parameters identified as ξ α = u α and P αα = P ex αα , where the extended parameter P ex αα reads, corresponds to the conventional product-form equilibrium, The effect of extension, featured by the second term in (3.13), is to correct for the incomplete Galilean invariance of the standard D3Q27 velocity set (2.10). With the original formulation of the mixture momentum lattice in Sawant et al. (2021a), a similar correction was achieved by augmenting Eq. (3.10) with an additional forcing term which required evaluation of second-order derivatives in space. In the present formulation, the correction of Galilean invariance is achieved by the extended equilibrium which requires evaluation of only a first-order derivative, cf. Eq. (3.13), a more local operation.

Equilibrium and quasi-equilibrium of the energy lattice
Turning our attention to the energy lattice, the corresponding equilibrium and quasiequilibrium populations in (3.11) are evaluated along the lines of Saadat et al. (2021): Let us introduce linear operators O α , acting on any smooth function A(u, T ) according to a rule, The equilibrium populations g eq i are specified with an operator version of the productform (2.14). To that end, we consider parameters ξ α and P αα as operator symbols, With the operators (3.17) and (3.18) substituted into the product form (2.14), the equilibrium populations g eq i are compactly written using the energy E as the generating function, It is straightforward to verify by a direct computation that the equilibrium (3.19) satisfies the necessary conditions to recover the mixture energy equation as in Sawant et al. (2021a), namely, the equilibrium energy flux q eq and the flux thereof R eq , where H is the specific mixture enthalpy (3.4). Finally, the quasi-equilibrium populations g * i differs from the equilibrium g eq i by the energy flux only Sawant et al. 2021a;Saadat et al. 2021), ( 3.22) were q * is a specified quasi-equilibrium energy flux, All contributions on the right hand side of (3.23), except for the vector q ex , were already introduced in Sawant et al. (2021a) and do not alter under the present modifications: The two first terms in (3.23) maintain a variable Prandtl number and include the energy flux q and the pressure tensor P , The interdiffusion energy flux q diff , where the diffusion velocities V a are defined by Eq. (2.31), contributes the enthalpy transport due to diffusion, cf. (Sawant et al. 2021a). Moreover, the correction flux q corr is required in the two-population approach to the mixtures in order to recover the Fourier law of thermal conduction (Sawant et al. 2021a), Finally, the term q ex in the quasi-equilibrium flux (3.23) is required for consistency with the extended equilibrium (3.12), and is similar to its single-component counterpart (Saadat et al. 2021). Components of the vector q ex follow the structure of (3.13), Spatial derivatives in the correction flux (3.27) and in the isotropy correction (3.13) and (3.28) were implemented using isotropic lattice operators (Thampi et al. 2013).

Mixture mass, momentum and energy equations
With the equilibrium and quasi-equilibrium populations specified, the hydrodynamic limit of the two-population lattice Boltzmann system (3.10) and (3.11) is found by expanding the propagation to second order in the time step δt and evaluating the moments of the resulting expansion. Analysis is standard, details can be found in Sawant et al. (2021a) and Saadat et al. (2021), here we present the final result. The continuity, the momentum and the energy equations for a reactive multicomponent mixture (Williams 1985;Bird et al. 2007) are, respectively, Here, the pressure tensor π in the momentum equation reads, where the dynamic viscosity µ and the bulk viscosity ς are related to the relaxation parameter ω, ,v is the mixture specific heat at constant volume. The heat flux q in the energy equation (3.31) reads, (3.35) The first term in (3.35) is the Fourier law of thermal conduction, with thermal conductivity λ related to the relaxation parameter ω 1 , where C p = C v + R is the mixture specific heat at constant pressure. The second term in (3.35) is the interdiffusion energy flux. With the thermal diffusivity α = λ/ρC p and the kinematic viscosity ν = µ/ρ, the Prandtl number becomes, Pr = ν/α. For this reactive formulation, the local dynamic viscosity µ(x, t) and the thermal conductivity λ(x, t) of the mixture is evaluated as a function of the local chemical state by using the chemical kinetics solver Cantera (Goodwin et al. 2018). Cantera employs a combination of methods such as interaction potential energy functions (Kee et al. 2003), hard sphere approximations, the methods described in Wilke (1950) and Mathur et al. (1967) to calculate the mixture transport coefficients. In summary, by virtue of thermodynamic consistency of the lattice Boltzmann model for mixture momentum and energy (Sawant et al. 2021a), the extension to the reactive case requires merely an upgrade of the sensible heat by the heat of formation. The proposed realization also takes into account the revised formulation of the two-population LBM for compressible flow (Saadat et al. 2021). We proceed to finalizing the model development by specifying the coupling between the lattice Boltzmann models for the species and the mixture momentum and energy, as well as the coupling to the external chemical kinetics solver.

Coupling between the species and the mixture subsystems
With the two subsystems, the species and the mixture, first constructed independently from each other and after that being coupled weakly in the way described in Sawant et al. (2021a), we are left with two independent definitions of the mixture density and the mixture momentum: On the one hand, the mixture density ρ (3.5) and the mixture momentum ρu (3.6) are defined as the moments of the f -populations on the momentum lattice. On the other hand, the same quantities are defined with the species populations as the sum of partial densities and partial momenta. The number of the conservation laws for the species subsystem is M + D, while for the mixture subsystem it is D + 2. The total number of the conservation laws in the weakly coupled combined system is M + 2D + 2. Thus, the weakly coupled system is in excess of D + 1 conservation laws. This redundancy is eliminated by removing one set of species populations (here, the M th ) and writing, (3.37) As a consequence, the M th component is not an independent field anymore but is slaved to the remaining species and mixture populations. The number of independent conservation laws in the resulting strongly coupled system is M +D+1, which corresponds to the locally conserved fields, ρ 1 , . . . , ρ M −1 (2.4), ρ (3.5), ρu (3.6) and ρE (3.7) . While the assignment of the slaved component M is not unique, it is advisable to select the component which carries the majority of mass in the mixture. The coupling (3.37) reduces the number of lattices from M + 2 to M + 1.

Coupling between lattice Boltzmann and chemical kinetics
The lattice Boltzmann code is coupled to the open source code chemical kinetics solver Cantera (Goodwin et al. 2018). The Cantera solver is supplied with the publicly accessible GRI-Mech 3.0 mechanism (Smith et al. 1999) as an input. The communication between the lattice Boltzmann solver and the Cantera is summarized as follows: (i) During the collision step, the lattice Boltzmann solver provides internal energy, specific volume and mass fractions to set the chemical state in Cantera. (ii) Cantera numerically solves the integral equation (3.9) to find the temperature at that state. (iii) The production rates of speciesρ r a , transport coefficients including dynamic viscosity, thermal conductivity and the Stefan-Maxwell diffusivities are obtained from Cantera as a function of the current state. (iv) In the lattice Boltzmann solver, the temperature is used to evaluate the equilibrium and quasi-equilibrium moments and populations. The transport coefficients are used to calculate the corresponding relaxation times.
Other thermodynamic parameters necessary for the simulation such as the specific heats and molecular masses are also obtained through Cantera. The reference standard state temperature is T 0 = 298.15 K and the reference standard state pressure is P 0 = 1 atm. The data required by the lattice Boltzmann solver during runtime is obtained by querying Cantera through its C++ API. In all cases considered in this paper, we use the detailed mechanism of hydrogen/air combustion (Li et al. 2004) Sawant et al. (2021a), acoustic scaling is used for conversion of length and time between the physical and the lattice units. The speed of sound at a specified reference composition and specified temperature (typically, at the unburnt mixture state) is used to make the velocity non-dimensional. The characteristic length in the respective setup is used to rescale the length.
We shall now proceed with a validation of the coupled reactive flow lattice Boltzmann model in two test cases. The perfectly stirred reaction (PSR) simulation is selected to validate the multistep approach to the evaluation of the reaction term (2.34) while the laminar flame speed simulation is to probe the coupling of the new formulation of the mixture momentum and energy LBM of sec. 3.2 and 3.3.

Perfectly stirred reactor
A constant volume PSR is simulated using LBM with a three-dimensional domain consisting of 4 × 4 × 4 nodes. Periodic boundary conditions are used in all directions. The computational domain is initialized with a stagnant and homogeneous hydrogen/air mixture at an equivalence ratio φ = 1, pressure P in = 1 atm and temperature T in = 1400 K. Fig. 1 shows the evolution of the temperature and of the hydroxide mass fraction in the reactor over time. The results from the lattice Boltzmann model are compared to the solution produced by the ideal gas constant volume reactor from Cantera. The time integration in Cantera is performed through its built-in 'advance' function. Accurate match with the results obtained from Cantera verifies that the coupling and the multistep time integration of the reaction term is correct. Since all the boundaries are periodic in this setup, the total energy of the system must remain constant. Also, due to the completely homogeneous initial condition, no kinetic energy should develop over time. Fig. 1 verifies that in the absence of flow, the total energy not only equals the internal energy but it also remains constant in time, as expected.

Laminar flame speed
For a further validation, we calculate the burning velocity of a hydrogen/air mixture. The setup consists of a one-dimensional tube initialized with unburnt mixture at T u = 300 K throughout from the left end up to 80% of the domain towards the right. The remaining 20% of the domain are initialized with the adiabatic flame temperature T af = 1642.49 K and with the equilibrium burnt composition at the respective equivalence ratio. The pressure is initialized uniformly at P in = 1 atm. The inlet and the outlet boundary condition used in this case will be explained below in sec. 4. At the left end, the inlet velocity is set to u in = 10 cm/s so that the flame propagates from right to left against the unburnt mixture.
We use the laminar flame thickness δ f = 0.043 cm at φ = 0.5 for defining the reference length, where δ f = (T af − T u )/ max(|dT /dx|). The domain size is N ≈ 23δ f with a resolution of 34 points per flame thickness. As evident in Fig. 2, the profiles of the temperature and the mass fractions for φ = 0.5 compare well with the solution obtained from the 'FreeFlame' solver of Cantera. The burning velocity is found to be S L = 59 cm/s which is in good agreement with the reference result of Pizza et al. (2008a), i.e. 58 cm/s. To summarize, the basic validation of the proposed LBM for reactive mixtures is considered successful. We now proceed with specifying various boundary conditions for the multicomponent LBM, needed for most of practical applications.

Nomenclature
Boundary conditions for multi-component LBM are scarce in the literature. In order to facilitate the explanation, we use the cartoon in Fig. 3, which represents a rectangular grid and empty circles represent grid points (nodes) which are part of the computational domain. The boundaries are marked by coloured dotted lines, where the colour reflects either the wall, the inlet or the outlet. The boundaries do not belong to the computational domain and therefore do not participate in the collision and the advection operations. During the advection step, a node at location x performs the following operation for each of the populations f i , (4.1) Equation (4.1) is a mathematical expression for the free streaming of a population f i by jumping a distance c i δt to a new node. Since we do not need to discuss the collision step in this section, the times t and t − δt simply indicate the post-and the pre-advection states, respectively. In Fig. 3, each population f i is represented by its corresponding discrete velocity vector c i (link) by an arrow pointing in the direction of its propagation. Solid arrows represent the post-advection populations that arrived from a node belonging to the computational domain. Dotted arrows represent the post-advection populations that have arrived from one of the boundaries and carry with them the information about the fluid properties at the boundary. These populations will be referred to as incoming populations since they enter the domain from the boundaries. In the lattice Boltzmann method, the boundary conditions are applied by specifying the incoming populations.
The nodes which are adjacent to the boundaries and therefore require such description for incoming populations will be referred to as the interface nodes. We denote D the set of the incoming velocities at the interface node. Finally, the rest of the velocities at the interface node c i / ∈ D are the outgoing velocities. Below, the equilibrium form shall be used to evaluate a variety of incoming populations. In order to keep the discussion concise, we shall display the dependence of pertinent equilibria on the respective control parameters as follows: Here Y = {Y 1 , . . . , Y M } stands for the totality of mass fractions. Dependence on the mixture composition Y in the energy lattice equilibrium (4.4) is manifest in the operational definition (3.19) through the mixture-averaged gas constant R in the operators (3.16) as well as in the mixture energy (3.8). The composition dependence Y enters the momentum lattice equilibrium (4.3) through the gas constant R, cf. Eqs. (3.14) and (3.15). We now proceed to derive the wall, the inlet and the outlet boundary conditions for the multicomponent LBM.

Bounce-back boundary condition
Bounce-back (BB) is a widely used wall boundary condition in the lattice Boltzmann method (Ladd 1994). For the incoming populations at interface node x, the bounce-back rule reads, (4.5) Here D is the set of incoming velocities shown by grey dotted arrows in Fig. 3. When applied on the momentum lattice, the bounce-back rule (4.5) results in the no-slip boundary condition at a half-way distance between the wall and the interface nodes (Ziegler 1993;Boyd et al. 2004). On the energy lattice, the bounceback boundary condition conserves the total energy and leads to zero heat flux, thereby representing an adiabatic wall (He et al. 1998). While simple and efficient, the bounceback boundary condition (4.5) is limited as it does not allow to impose a target value for the velocity at a prescribed wall location nor to implement a target wall temperature. Since these are the cases typical of many applications, including the ones considered below, a so-called Tamm-Mott-Smith (TMS) boundary condition of  shall be adapted to the multicomponent mixture.

Tamm-Mott-Smith wall boundary condition
Let u tgt , u tgt a and T tgt be the target values of the flow velocity, species velocity and the temperature, respectively, to be imposed at the interface node x. Moreover, the outgoing populations f i , g i and f ai , where i / ∈ D, are obtained in the propagation step (4.1) and assumed known. The TMS construction of the incoming populations f TMS i , g TMS i and f TMS ai , where i ∈ D, executes the following steps: (i) Perform bounce-back on the momentum lattice to find the densities ρ bb , Perform bounce-back on the species lattices to find the mass fractions Y bb , Note that the bounce-back operation is used solely for computing the density ρ bb (4.6) and the mass fractions Y bb (4.7), in order to satisfy mass conservation at the boundary. However, the incoming populations are not set to the bounce-back values f bb i , rather, they are defined with the subsequent steps of the TMS algorithm.
(ii) The bounce-back density ρ bb (4.6) and mass fractions Y bb (4.7), together with the target velocities u tgt , u tgt a and temperature T tgt , uniquely specify the equilibrium states f tgt i , g tgt i and f tgt ai at the interface node, f tgt ai = f eq ai (ρ bb Y bb a , u tgt a , T tgt ).
(4.10) (iii) With the incoming populations set to the target equilibrium, we find the local density ρ loc , flow velocity u loc , mass fractions Y loc , species velocities u loc a and temperature T loc at the interface node, i∈D f tgt (4.14) i∈D C a,v (T )dT + 1 2 ρ loc u loc 2 . (4.15) We remind that Eq. (4.15) is an integral equation to be resolved for the temperature T loc , cf. sec. 3.1, Eq. (3.9). With the local parameters, the following equilibrium populations are uniquely specified on the momentum, energy and species lattices, f loc ai = f eq ai (ρ loc Y loc a , u loc a , T loc ).
(4.18) (iv) Finally, we update all populations of the momentum, the energy and the species lattices at the wall interface node as: Comments are in order. The TMS boundary condition in step (iv) sets the flow variables at the interface nodes to ρ bb , Y bb , u tgt a , u tgt and T tgt . While the same is achieved by the target equilibrium at step (ii), the corresponding equilibrium boundary condition is insufficient as it is prone to generating spurious shocks, cf. . For this reason, a non-equilibrium part of the incoming populations is taken into consideration and modelled with the local state in step (iii). Note that, while the latter also uses the equilibrium form, it is evaluated at different (local) values of flow variables and thus describes a non-equilibrium state relative to the target equilibrium. The presence of two different equilibrium states in the resulting populations motivated  to naming the algorithm in analogy to the bimodal Tamm-Mott-Smith shock wave approximation for the Boltzmann equation (Mott-Smith 1951).
With the exception of walls aligned with the Cartesian LBM grid, target parameters u tgt , u tgt a and T tgt at the interface nodes are obtained by interpolation between the values of the corresponding fields at the wall, u wall , u wall a and T wall , and the data at the surrounding fluid nodes. Interpolation is performed following the procedure described in Dorschner et al. (2015), examples shall be demonstrated below in sec. 5.2 for stationary no-slip walls, u wall = 0, subject to a temperature profile. Finally, the impermeable wall boundary condition is imposed on the species populations by setting the species velocity at the wall as u wall a = u wall . Zero flux of species at the wall is implied by the absence of diffusion velocity in the species equilibrium velocity u wall .

Inlet
The flux boundary condition is widely used to model the inlet in multicomponent flows (Kee et al. 2003;Pizza et al. 2010Pizza et al. , 2008aGoodwin et al. 2018). The rationale of this boundary condition is that it prescribes only the incoming mass fluxes of species ρ in a u in . Because only the incoming mass flux is prescribed and not the mass itself, the composition at the inlet interface node is not fixed to the incoming composition Y in . This degree of freedom is necessary as light species such as hydrogen have the capability to diffuse fast enough and thus are able to propagate upstream into the inlet. Therefore, the composition at the inlet interface node is not a fixed set of parameters but is rather a result of a balance between the mass flux inside the domain and the inlet mass flux. Below, we establish the flux boundary condition for the multicomponent lattice Boltzmann setting.
In Fig. 3, the inlet boundary is represented by a dotted vertical blue line. The inlet boundary condition is applied on the interface nodes where incoming discrete velocities c i , i ∈ D, are represented by dotted blue arrows. With the inlet data for mass flux ρ in u in , composition Y in and temperature T in , populations at the interface node are derived in the following steps: (i) The inlet density ρ in and composition Y in , together with the inlet velocity u in and temperature T in , uniquely specify the inlet equilibrium populations f in i , g in i and f in ai at the inlet interface node, (4.24) (ii) With the incoming populations set to the inlet equilibrium and the outgoing populations known, we find the local density ρ loc , flow velocity u loc , composition Y loc and temperature T loc at the interface node, (4.28) With the local parameters, the following equilibrium populations are uniquely specified on the momentum, energy and species lattices, f loc ai = f eq ai (ρ loc Y loc a , u loc , T loc ).
(4.31) (iii) Replacing the local flow velocity u loc and temperature T loc with the target values u in and T in , the following target equilibrium populations are identified, (4.34) (iv) Finally, all populations at the inlet interface nodes are updated as follows: It is straightforward to verify that the populations (4.35), (4.36) and (4.37) at step (iv) imply the target values u in and T in for the velocity and temperature at the interface node, respectively. At the same time, the composition and density at the interface node are identified as Y loc and ρ loc , respectively. The latter are derived in step (ii) by taking into account the outgoing populations and are different, in general, from the inlet values Y in and ρ in . Thus, the outgoing populations contribute to the balance between the incoming and outgoing mass fluxes as required by the flux boundary condition. It is instructive to compare with the TMS wall boundary condition of sec. 4.2.2 where the local composition at the interface node was determined by the bounce-back step, Eqs. (4.6) and (4.7). In the present case, at step (ii), the local composition is computed using the equilibrium at the inflow properties for the incoming populations. Note that, while the inlet composition Y loc is already determined at step (ii), the purpose of the remaining steps is to enforce the inlet velocity u in and temperature T in . Hence, whenever the inlet temperature and velocity need not be strictly imposed, it is sufficient to terminate the algorithm at step (ii) and to apply the local equilibria (4.29), (4.30) and (4.31). With this simplification, the velocity and temperature acquire local values u loc and T loc , respectively, rather than the target values u in and T in . The latter simplification was validated in Sawant et al. (2021a) with the simulation of diffusion in opposed jets. The mass fractions at the inlets of both jets matched the reference solution by Cantera, which employs a macroscopic realization of the flux boundary condition. While the simplified inlet realization (4.29), (4.30) and (4.31) can be regarded as a good approximation to the macroscopic flux boundary condition, in this work we rather use the inlet populations (4.35), (4.36) and (4.37) to ensure that the inlet velocity and temperature are imposed exactly.

Outlet
Unlike the inlet and the wall, the values of the macroscopic state variables are usually unknown at the outlet. To that end, we apply the Local One Dimensional Inviscid (LODI) approximation by Poinsot & Lele (1992). LODI is based on the characteristics of compressible Euler equations, i.e. Eqs. (3.29), (3.30) and (3.31) without dissipation terms. LODI boundary condition allows both the pressure fluctuations travelling as sound waves as well as the convection disturbances travelling as entropy waves to exit the computational domain with minimum reflection (Poinsot & Lele 1992). The LODI approximation is derived for a single-component fluid and therefore predicts the outlet density ρ out , velocity u out and the temperature T out which can be directly used in the present mean field formulation of the mixture. In addition, we need also to specify the composition Y out at the outlet. Consistent with the LODI approximation, we use the advection part of the species equation (2.22) which is discretized at the outlet interface node with forward Euler scheme to give, where mass fraction Y out a (x, t − δt) is known from the previous time step, while u out is the LODI outlet velocity. The gradient ∇Y out a is evaluated by backward finite difference. Armed with the outlet data, we proceed to specify the populations at the outlet interface node, following essentially the steps already familiar from the wall and inlet construction: (i) Outlet data Y out a , ρ out , u out and T out uniquely specifies the equilibrium populations f out i , g out i and f out ai at the outlet interface node, f out i = f eq i (ρ out , Y out , u out , T out ), (4.39) g out i = g eq i (ρ out , Y out , u out , T out ), (4.40) f out ai = f eq ai (ρ out Y out a , u out , T out ). (4.41) (ii) With the incoming populations set to the outlet equilibrium, we find the local density ρ loc , flow velocity u loc , mass fractions Y loc , and temperature T loc at the outlet interface node, Based on these local parameters, the local equilibrium populations are uniquely specified on the momentum, energy and species lattices, f loc ai = f eq ai (ρ loc Y loc a , u loc , T loc ).
(4.48) (iii) Finally, all populations of the momentum, energy and species lattices at the outlet interface node are updated as: With populations (4.49), (4.50) and (4.51) at step (iii), the macroscopic fields at the outlet are set to the target values ρ out , u out , T out and Y out , as prescribed by LODI approximation and (4.38). Although the same is achieved by the equilibrium populations at step (i), the non-equilibrium part of the incoming populations is taken into consideration and modelled with the local state in step (ii). Thus, the present construction of the outlet is similar to the TMS wall boundary condition of sec. 4.2.2.

Wall-bounded reactive flow
In order to test the proposed boundary conditions, we perform the computation of combustion in microtubes. The results are validated with the direct numerical simulation of Pizza et al. (2008b) for 2D microchannels and with Pizza et al. (2010) for a 3D microtube. The setup involves combustion of a premixed hydrogen/air mixture in a tube over a range of inlet velocities u in of the unburnt mixture. The fuel-lean unburnt mixture of equivalence ratio φ = 0.5 at a temperature T u = 300 K and pressure 1 atm enters a microchannel with l/d = 10. Here, l is the length of the tube and d is its diameter. The mixture gets ignited due to hot isothermal walls which are maintained at a temperature of T w = 960 K. The wall temperature is increased from 300 K at the inlet to 960 K using a hyperbolic tangent profile at a distance of about l/20 from the inlet. For this premixed initial condition, the burning velocity is obtained as S L = 59 cm/s and the flame thickness is obtained as δ f = 0.043 cm from solving a 1D flame propagation setup with the lattice Boltzmann method in section (3.8).

Premixed hydrogen/air flames in a microchannel
For the 2D simulations, we choose a channel diameter d = 1 mm, which corresponds to a width of 2.3256 δ f in terms of flame thickness. The spatial resolution corresponds to approximately 15 nodes per flame thickness. As studied in Pizza et al. (2008b) for the same channel width, the flame exhibits different dynamics depending on the inlet velocity. At low inlet velocity of about u in = 10 cm/s, periodic ignition and extinction of the flame is observed. The inlet velocity is then progressively increased until the oscillatory behaviour ceases and a stable flame can be sustained near the inlet of the channel. A further increase of the inlet velocity to u in = 75 cm/s results in a symmetric "V-shaped flame" in the channel, the flame being concave towards the unburnt mixture. Finally, at inflow higher than u in = 165 cm/s, stable asymmetric flames are formed which shift downstream with increasing inlet velocity. The Reynolds number corresponding to the inlet velocity varies between Re = 5.36 to Re = 88.52, the reference length being the channel width and the reference viscosity corresponding to the viscosity at inlet composition and temperature.

Periodic ignition and extinction
The fluid in the bulk of the domain is initialized with the inlet unburnt composition and the inlet velocity is set to u in = 10 cm/s. The initial temperature of the fluid in the bulk follows the wall temperature profile. As the fresh mixture passes between the heated walls, the reactants break into radicals which build up in the channel over time. This build up of radicals is associated with a long period of inactivity after which the mixture achieves a radical runaway and eventually a thermal runaway, leading to ignition. The mixture ignites at some distance downstream as seen in Fig. 4. The Fig. 4 shows the hydroxide mass fraction which we will use as a marker to represent the "flame" itself. The flame first forms a concentrated nearly circular structure which then propagates in both upstream as well as the downstream direction as visible in Fig. 4. This flame splitting occurs as the flame consumes the relatively fresh mixture in both possible directions. The frames in Figs. 5 show the mass fraction of hydrogen which is the deficient reactant. The flame then propagates and splits, consuming the deficient reactant in its path. The part of the flame travelling upstream is extinguished at the cold inlet whereas the part travelling downstream exits the channel through the outlet. Subsequently, the channel is again filled with the fresh mixture from the inlet and the process repeats periodically. In this regime, the maxima of all the species mass fractions as well as that of the temperature is located on the centreline of the channel. The behaviour is consistent with the DNS of Pizza et al. (2008b) and the subsequent simulations of Alipoor & Mazaheri (2016). This phenomenon which is also referred to as a flame with repetitive extinction-ignition (FREI) has also been observed in methane/air combustion experiments of Maruta et al. (2005) and numerical simulations of Norton & Vlachos (2003). The periodicity of the ignition- extinction behaviour has been presented through the variation of the integrated heat release rate with time in Fig. 6. In the figure, the heat release rate has been normalized with respect to the heat release rate of the unburnt state. The ignition events are seen produce a rise in the integrated heat release rate by 8 orders of magnitude. The peaks are localized in time with an average frequency of approximately 111 Hz. This is in good agreement with the frequency of 106.9 Hz reported in Pizza et al. (2008b). Table 1 shows the convergence of the ignition-extinction frequency with resolution. The frequency changed by only 2.7% with an increase in the spatial resolution by 50%. Therefore, the computations have been considered to be converged with respect to the resolution. The maximum velocity of the upstream propagation of the flame is found to be 16.8 m/s and that of the downstream propagation is found to be 20 m/s from the LBM simulations. For comparison, the maximum upstream propagation speed is reported to be 15 m/s in Pizza et al. (2008b). Overall, the LBM results quantitatively agree well with the DNS results.

V-shaped stable flames
Using the solution from the ignition-extinction regime as an initial condition, the inlet velocity is progressively increased to u in = 75 cm/s. In this regime, there is a sufficient flow of fresh mixture to sustain combustion and therefore a stable flame is formed in the channel. As evident in Fig. 7, the flame assumes a "V-shaped" structure which is concave towards the unburnt mixture. At this inlet velocity, the structure of all the species is symmetric about the centreline. The maxima of the mass fraction of all the species is located on the centreline except for the hydrogen radical. The hydrogen radial has a high molecular diffusivity, causing it to shift away from the channel centreline. This is evident from the line contours in Fig. 8. The heat release rate contours in Fig. 9 show a localized  heat release at the upstream interface of the flame. Also, the heat release rate contour follows a concave curvature that is similar to mass fraction contours of the hydroxide and the hydrogen radical. A maximum temperature of 1715 K is attained in the flame. Shifting of the maxima of hydrogen at this inflow velocity as well as the concave shape   of the flame is consistent with the findings of Pizza et al. (2008b). With an increasing inlet velocity, the flame stabilizes further downstream from the inlet due to a relative increase in the difference between the velocity of the fresh mixture and the flame speed. Furthermore, at higher inlet velocities, more species shift away from the tube centreline. In a 3D microtube, the shifting of the maxima causes the flame to form a ring like structure (Pizza et al. 2010) around the tube centreline. We explore this phenomenon in detail in the corresponding 3D simulations.

Asymmetric stable flames
Starting from the V-shaped flame as an initial condition, we increase the inlet flow velocity gradually to u in = 300 cm/s. After shifting downstream and maintaining its symmetric shape for some time, the flame transitions into an asymmetric stable flame. At the upstream interface between the flame and the unburnt mixture, a flame forming an acute angle with the lower wall is termed as a lower asymmetric flame (Pizza et al. 2008b). Similarly, a flame forming an acute angle with the upper wall is termed an upper asymmetric flame. As shown in Fig. 10, a lower symmetric flame was first encountered in our computation. Interestingly, the asymmetric flame is metastable in this regime. The flame can be made to transition from a lower asymmetric shape to an upper asymmetric shape by heating the lower wall momentarily and then restoring the wall temperature  back to the previous wall temperature of 960 K. A snapshot of the flame during this transition is shown in Fig. 11. An upper asymmetric flame formed as a result of the temperature perturbation is shown in Fig. 12. The resultant flame is also metastabe and remains in its upper asymmetric shape unless perturbed. The heat release rate profile is very similar to the profile of the mass fraction of hydrogen, except for the location of the maxima. The maxima of the heat release rate occurs at the walls in this regime. The distance of the flame from the inlet remains unchanged. In the LBM simulations, the location of the beginning of the flame is 0.24 l from the inlet, which is in good agreement with 0.21 l, as obtained by the DNS of Pizza et al. (2008b). The asymmetric nature of the flame and its metastable behaviour is consistent with the findings of Pizza et al. (2008b) for this regime.

Premixed hydrogen/air flame in microtube
For three-dimensional simulations, we choose the circular tube setup with a width d = 1.5 mm from the simulations in Pizza et al. (2010). Richer dynamics are exhibited by this wide tube as compared to the narrower 1 mm tube. The composition of the incoming mixture is the same as for the two-dimensional simulations. The aspect ratio is also unchanged at l/d = 10. The diameter corresponds to a flame thickness of d = 3.4884δ f . With the viscosity of the mixture at the inlet unburnt composition and at the inlet temperature as reference, the inlet velocity as the reference velocity and the diameter of the tube as the length scale, the Reynolds number is Re = 80.47. We use a computational domain of l × d × d = 540 × 54 × 54 nodes which translates into a spatial resolution of 15.5 computational nodes per flame thickness. As discussed in section (5.1) with the aid of table 1, this resolution was sufficient to produce correct and converged Figure 11. Flame transitioning from a lower asymmetric shape to an upper asymmetric flame: Contours of hydroxide mass fraction (top), hydrogen mass fraction (middle) and heat release rate (bottom). results in the two-dimensional ignition-extinction simulations. In the three-dimensional simulations, we explore the 'open axisymmetric flame' characterized by the maximum of the hydroxide radical being shifted away from the tube centreline. The iso-surfaces of hydroxide therefore form a ring shaped structure around the tube centreline. In the DNS of Pizza et al. (2010), such open flames are observed for an inflow velocity over two disconnected ranges, 60 cm/s u in 100 cm/s and 170 cm/s u in 350 cm/s. In this work, we verify the existence of open flames by performing a simulation at an inflow of u in = 100 cm/s. The bulk of the fluid in the tube is initialized with the inflow velocity u in = 100 cm/s and the temperature is initialized to the wall temperature profile. The composition is initialized with a 1D laminar flame solution computed for the same equivalence ratio as this 3D setup. The initial pressure in the domain is homogeneous at 1 atm. As a consequence of the initial conditions and the inflow velocity, the ignition-extinction regime is not encountered. The incoming fresh mixture enters the tube at a speed which is nearly twice of the flame speed. Therefore, the resulting flame does not propagate  upstream into the inlet and stabilizes at a distance downstream of the inlet. The flame has the location of the maximum of the hydroxide radical and the temperature shifted away from the longitudinal axis of the tube. Therefore, as visible in Fig. 13, iso-surfaces of the mass fraction of hydroxide form ring shaped structures. This type of a ring-like flame is called an 'open flame' in Pizza et al. (2010). At u in = 100 cm/s, the open flame is axisymmetric and maintains a fixed distance of approximately 0.1l from the inlet. The maximum temperature in the flame is 1649 K. Fig. 14 shows iso-surfaces of the mass fraction of the hydrogen radical, which also forms a ring structure similar to the hydroxide radical. Streamlines of the fluid velocity in Fig. 14 show a marked acceleration in the fluid velocity downstream of the flame location. The post combustion fluid in the tube is seen to have attained a maximum velocity which is 9.3 times of the inlet velocity. The maxima of the fluid velocity resides on the tube centreline.

Conclusion
In this paper, we aimed at developing an accurate and robust LB model for reactive flows of practical interest. In Sawant et al. (2021a), we proposed a lattice Boltzmann framework consisting of M +2 lattice Boltzmann equations for multicomponent mixtures of ideal gases. We introduced a new LBM system for the Stefan-Maxwell diffusion along with a reduced, mean-field description of the mixture momentum and energy using a twopopulation approach. Thermodynamic consistency of this model allowed us to naturally account for the temperature and energy changes due to chemical reactions by including the energy of formation, which avoids any ad-hoc modelling for the heat of reaction. The proposed model uses the extended lattice Boltzmann method of Saadat et al. (2021) for the mean fields and a multistep approach for integrating the mass source terms. Furthermore, we introduced novel kinetic boundary conditions for walls, inlets and outlets that are compatible with the underlying reactive flow model.
Our model was validated in detail, starting from a zero-dimensional perfectly stirred reactor and the one-dimensional laminar flame. The accuracy of the boundary conditions was assessed by simulations of premixed hydrogen/air flames in a microtube in both two and three dimensions. In all cases the results were found to be in excellent agreement with reference and DNS solutions that can be found in the literature.
To conclude, the proposed model is not only a viable alternative to traditional reactive computational fluid dynamics but it is also the first model capable of solving reactive flows entirely in the lattice Boltzmann framework. We believe that this model marks a significant stage in the development of LBM by expanding the applicability of LBM to a wide range of setups including but not limited to diffusion, reactive flows and combustion.