A Non-Isothermal Chemical Lattice Boltzmann Model Incorporating Thermal Reaction Kinetics and Enthalpy Changes

The lattice Boltzmann method is an efficient computational fluid dynamics technique that can accurately model a broad range of complex systems. As well as single-phase fluids, it can simulate thermohydrodynamic systems and passive scalar advection. In recent years, it also gained attention as a means of simulating chemical phenomena, as interest in self-organization processes increased. This paper will present a widely-used and versatile lattice Boltzmann model that can simultaneously incorporate fluid dynamics, heat transfer, buoyancy-driven convection, passive scalar advection, chemical reactions and enthalpy changes. All of these effects interact in a physically accurate framework that is simple to code and readily parallelizable. As well as a complete description of the model equations, several example systems will be presented in order to demonstrate the accuracy and versatility of the method. New simulations, which analyzed the effect of a reversible reaction on the transport properties of a convecting fluid, will also be described in detail. This extra chemical degree of freedom was utilized by the system to augment its net heat flux. The numerical method outlined in this paper can be readily deployed for a vast range of complex flow problems, spanning a variety of scientific disciplines.


Introduction
The lattice Boltzmann method (LBM) is a kinetic-based computational fluid dynamics (CFD) technique that was traditionally viewed as a somewhat esoteric, alternative paradigm for the simulation of hydrodynamic systems.It originated from the lattice gas cellular automata (LGCA) [1][2][3], which is a discrete fluid model involving the movement and collision of particles on a lattice.The LGCA attracted interest because it was based on a fictitious mode of molecular interactions, but was nonetheless shown to reproduce the standard Navier-Stokes equations.However, due to its particulate nature, it suffered from considerable stochastic noise and required coarse-graining and other modifications to be mapped onto many systems of interest.It remains useful for exploring concepts of statistical mechanics and dynamical systems theory [4,5], but is less popular for conventional CFD problems.
When the LGCA is coarse-grained, the propagation and collision of individual particles become the propagation and modification of probability distribution functions for a discrete set of particle velocities [2,3].Fluid mass is encoded in these velocity distribution functions, which move under their own characteristic velocity vectors and exchange mass at discrete grid points due to inter-particle collisions [6][7][8][9][10].Thus, the LBM is a kinetic-based model, which does not solve the equations of motion for a fluid in the continuum limit, but instead solves for the evolution of velocity distribution functions.With a discrete set of velocities and a suitably-discretized spatial and time domain, the LBM can accurately model a range of flow systems, and a suitable multi-scale analysis can be used to derive the Navier-Stokes equations from the governing equations of the LBM [10,11].
Given its kinetic nature, the LBM can also be expanded to include the advection and diffusion of internal energy [12][13][14][15][16][17][18], alongside the advection and diffusion of momentum.Heat is included as a passively advected scalar quantity, with an extra set of distribution functions.By extension of this logic, further passive scalars can be added.This allows multi-phase and multi-component fluids to be simulated [19][20][21][22].Versions in recent years have developed yet further and incorporated a range of interactions between these additional components, advancing the LBM into the realm of chemistry [23][24][25][26][27], including thermal chemical problems such as combustion [28][29][30][31].
There are two objectives of the present work: (1) to provide a pedagogical mini review of a general LBM framework that has become a commonly-deployed model for thermal and chemical flow systems; (2) to present new results from an investigation of the transport properties of a convective, reacting fluid system.The key features and equations of the numerical method will be described, and several example systems will be presented for the purpose of validating the method's efficacy.
The first sections will present the single-phase LBM, followed by the two-population thermal version, and example applications will be discussed.The introduction of additional components and chemical reactions will then be described.A simple, analytically-tractable decay reaction will be simulated, and it will be shown that the reactive LBM (RLBM) is second order accurate, in line with its single-phase counterpart.A simple method for extending the RLBM to include thermal reaction rates and enthalpy changes will then be outlined.The penultimate section explores a system that addresses a simple question in the field of transport phenomena: to what extent can the presence of chemical species and reactions alter the heat flux dynamics of a convecting fluid?The paper concludes with a brief discussion, and suggestions are given for the many applications to which this form of LBM could be amenable.

Single-Phase Fluid
For simulating the hydrodynamics of single-phase fluids, the most basic LBM is a relatively simple CFD method.As a kinetic-based scheme, it can either be thought of as deriving directly from the non-equilibrium Boltzmann equation: or as the result of a coarse-grained averaging of the LGCA [1][2][3].Equation (1) describes changes in the velocity distribution function f as a result of self-advection and particle collisions (represented by the function Ω( f )).The expression above can be recast into a system with a finite set of discrete velocities (e.g., the D2Q9 system in two dimensions, see below), which can then be discretized in space and time.
The exact discretisation method depends on the type of LBM being used (see, e.g., Dellar [32]); in this work, a hidden, semi-implicit scheme is employed.Assuming a simplification of the collision term in the above equation, wherein the velocity distribution functions relax towards their respective equilibria (given by a suitable discretization of the Boltzmann distribution) at a single characteristic rate τ ν , we arrive at the governing equation for the LBM [3,7,10,16,33]: where i indexes the discrete velocity set, f i is the velocity distribution for velocity i (a measure of the fluid mass possessing velocity i at that point in space and time), e i is the velocity vector, ∆t is the time step, τ ν is the relaxation time scale and f eq i is the equilibrium distribution: where ω i is the lattice weight for velocity i, c = ∆x/∆t is the lattice speed and C s = √ RT = c/ √ 3 is the speed of sound.The use of a single-time relaxation step is commonly known as the lattice Bhatnagar, Gross and Krook (LBGK) method [34].
Fluid properties are calculated from the appropriate moments of the distribution function, including density and velocity (higher order moments can also be used to calculate internal energy density and viscous stresses, etc.), All modeling in the present work is carried out using the two-dimensional D2Q9 velocity set, which utilizes a square lattice with eight velocities and rest particles.For this set, the velocity weights are ω 0 = 4/9, ω i = 1/9 for i = 1, 2, 3, 4 and ω i = 1/36 for i = 5, 6, 7, 8.The velocity vectors are thus: e 0 = (0, 0) e 1,3 , e 2,4 = (±c, 0), (0, ±c) e 5,6,7,8 = (±c, ±c) Using the Chapman-Enskog expansion, it can be shown that at the macro-level, a fluid obeying the above equations satisfies the Navier-Stokes equations (see, e.g., Chen and Doolen [7], Succi [10]) with a kinematic viscosity given by Many years of analysis have revealed that the above-described model can reproduce incompressible flow dynamics to second order accuracy [6][7][8][9][10]35].Furthermore, the LBM can be applied to turbulent flows [6,10,36], and recent novel lattice schemes [37], entropic multi-speed approaches [38] and cumulant methods [39] have further broadened the range of applications amenable to simulation by the LBM.

Thermohydrodynamics
As the LBM developed further, researchers began probing its ability to simulate thermohydrodynamic, multi-phase and multi-component systems.What emerged were two main methods for simulating the internal energy density of a modeled fluid.Given the intrinsic connection between local temperature and particle velocities, one can derive the temperature from the appropriate moment of the velocity distribution function [14,16,40]: where D is the spatial dimension, R is the molar gas constant, f is the distribution function before discretisation, v is the particle velocity, u is the net fluid velocity and ĝ is the distribution function for internal energy density.The integral is taken over the entire (continuous) velocity space.
After discretisation, this integral becomes: It is possible to construct simulations in which only the evolution of the velocity distribution function needs to be modeled.However, one must include more than one characteristic relaxation time and the method suffered from various issues including instability and the presence of spurious physical effects [41][42][43].Note that this approach was more recently developed into a fully-consistent method using the entropic lattice Boltzmann approach [13].
The alternative to multi-speed models is to track the internal energy distribution using an additional set of distribution functions.He et al. [14] developed a thermal LBM wherein the function ĝ is discretized in velocity space, physical space and time in the same way that f is discretized in the standard LBM.This method was then simplified by Peng et al. [18] into the following form.
The discretized internal energy distribution function evolves analogously to the velocity distribution function: where τ c is the characteristic relaxation time for internal energy.This parameter controls the thermal diffusivity: The equilibrium distribution functions for the internal energy in the D2Q9 velocity system are given by [18]: Many thermohydrodynamic problems are related to natural convection, in which buoyancy forces produce fluid accelerations due to temperature (and hence density) differences.Such an effect can easily be incorporated into the LBM framework using the Boussinesq approximation, in which it is assumed that density differences are negligible except with regard to the generation of buoyancy forces.The body force is represented by an additional term in the evolution equation for the velocity distribution function, which subsequently becomes: The body force term is given by: where T 0 is the mean temperature, G = βg 0 (T − T 0 ) ĵ is the vector representing the buoyancy force, β is the thermal expansion coefficient and g 0 is gravitational acceleration [14,18].As with the single-phase system, this model has been validated over many years by repeated analyses and comparisons with analytical and experimental results [14,16,18,44,45].It has been successfully applied to a variety of thermal flow problems including natural convection [14,18,44,[46][47][48][49][50], turbulent convection [51][52][53][54], thermal channel flows [14,40,55,56] and more complex systems involving multiple phases and phase change [57][58][59][60][61][62][63][64][65].
There have also been efforts to understand the relationship between natural convection and the principle of maximum entropy production (MEP) [66][67][68] (see Dyke and Kleidon [69], Martyushev and Seleznev [70] and references therein for details of the MEP principle).It was found that while previous authors proposed that natural convection was a prime example of MEP [71,72], in fact, this is not the case.The MEP prediction does not take proper account of the physical properties of the system in question (hence its generality), but convective heat flux is a strong function of the parameters of the system and is far more constrained than atmospheric flows (in which MEP predictions can be quite accurate [73][74][75]).A variety of counter examples, generated using LBM simulations, showed that MEP does not in fact predict the steady states of a convective fluid [66,67,76].

Dissolved Chemical Species
As well as internal energy, the LBM can readily include further additional components (see, e.g., [19,21]) (note that the original work in this direction was carried out by Gunstensen et al. [20]).There are a variety of approaches to this problem of multi-component fluids [22,23,25,29,30,[77][78][79], but arguably the simplest and most commonly-used involves a simple addition of distribution functions corresponding to each extra component in the fluid.The physical assumption underlying this approach is that the additional components are low in concentration such that they do not explicitly influence the fluid flow.They are classed as passive scalars and are advected by any net fluid velocity.Under this assumption, any number of extra species can be added, and they follow analogous evolution equations to the internal energy: ), where σ is the component index and τ σ is the relaxation time.The equilibrium distributions take on the following form: where ψ σ represents the concentration of species σ.Note that unlike Equation (3), higher order terms can be neglected from these distributions, since they were shown to be unnecessary by Ayodele et al. [24] (assuming the Mach number is kept low).Those authors found that both the linear and quadratic forms of the equilibrium distributions were second order accurate.The diffusivity of species σ is given by: and the concentration by:

Isothermal Reactive LBM
The RLBM described above requires further modifications to allow interactions between the components.The simplest way to introduce chemical reactions is to add extra terms to Equation (15).It will be assumed that the chemical kinetics can be described using a standard mass action rate law.

Single Reaction Benchmark Test
As an example test case, the irreversible decay of a reactant A → B in a closed system with no fluid motion is particularly favorable.Such a system can be described by a pair of coupled differential equations: which can be solved analytically using Fourier transforms.This yields an expression for the concentration field of species A: where it has been assumed that the mass of species A is initially concentrated at the point (x 0 , y 0 ), is the Thiele modulus, a dimensionless parameter that measures the ratio of the characteristic reaction time to that of mass diffusion.
For systems that are very large, have a high rate constant or a low mass diffusivity, φ 2 1, and the reaction will be the dominant process.If φ 2  1, because of a combination of small characteristic length, low reaction rate and/or high diffusivity, the dynamics of the system are dominated by the process of diffusion, and reaction plays only a minor role.In this work, It is straightforward to construct the RLBM equations for modeling this system.The two passive scalar species are governed by the following equations: ) ) where k is the rate constant.The equilibria for this system are simplified due to the lack of fluid motion (u = 0 and ρ = 1 everywhere): The RLBM results can now be compared to the exact solution of Equation ( 21), as was done by Ayodele et al. [23].The present assessment differs slightly in that all simulations were initialized at dimensionless time t 0 = 1 × 10 −4 and run through to t = 2 × 10 −4 .This was motivated by the desire to avoid initial conditions with large discontinuities.Such discontinuities are inevitable if the simulation is initiated simply with unit concentration of the reactant at the central grid node, i.e., when . The numerical error is calculated using: and is plotted for a range of relaxation parameters τ A and rate constants k in Figure 1.
Note that in order to compare the same problem for different parameter combinations, the dimensionless group φ 2 and the dimensionless end time t must be kept constant.Thus, the characteristic length L (grid size) was appropriately adjusted for each [τ A , k] pair.
As one might expect, smaller reaction rates and larger diffusivities produce lower errors.The errors are very similar in magnitude to those found by Ayodele et al. [23]; however, there are no minima in the error curves here.Using a third-order Chapman-Enskog expansion, Ayodele et al. [23] derived an expression for E as a function of various system parameters and found crucially that it is a quadratic function of the relaxation parameter τ A , but also contains several derivatives of the density ψ A .They concluded that since their error curves were not parabolic, errors were clearly entering into the system due to derivatives of the density.The different shapes of the curves in Figure 1 compared to those of Ayodele et al. [23] are probably related to the different initial conditions used in the present work.However, in agreement with Ayodele et al. [23], it was found that there is a linear dependence between the error and the rate constant k.In conclusion, it is clear that the RLBM is capable of simulating simple reactive systems.The appropriate Chapman-Enskog expansions have been shown to recover the relevant reaction-diffusion equations, and the method is second order accurate, in line with the single-phase version [23,77].Variants of this RLBM have been used to model a range of chemical flow systems including the velocity and chemical species distribution around a reacting block in a channel flow [80], reactive dissolution and precipitation of solid and porous materials [81][82][83][84], as well as heterogeneous catalysis in microflows [85,86].

Pattern Formation in the Gray-Scott System
This section will present a more detailed application of the RLBM to a non-linear, pattern-forming chemical system known as the Gray-Scott reaction diffusion system (GSRDS) [87][88][89].Its defining characteristic is an autocatalytic reaction that allows stable or oscillatory structures to form in the concentration fields of the chemical species.As the two main parameters of the system are varied, a vast suite of different patterns is exhibited, including plane waves, spiral waves, lamellar phases and circular solitons, or spots [76,[88][89][90][91][92][93]].Such a characteristic example of non-equilibrium pattern formation is a prime candidate for testing hypotheses regarding the thermodynamics of dissipative structures and the relationship of such structures to the self-organization processes that contributed to the origin of life [94,95].Understanding such self-organization processes could also guide future work in the fields of microfluidics and nanotechnology [96].
The majority of past investigations of the GSRDS were analytical or numerical, due to its interesting mathematical properties as a non-linear dynamical system, capable of exhibiting spatio-temporal chaos and complex pattern formation.These features made it an ideal exemplar of the original ideas of Turing [97].Despite the simplicity of the GSRDS, experimental work has shown that it can be readily mapped to the ferrocyanide-iodate-sulfite reaction, which exhibits a similar range of oscillatory and pattern-forming behaviors [90,98,99].
The GSRDS is described by the following equations: where the second RHS terms represent the autocatalytic reaction and the third terms implement porous wall BCs.These BCs act to maintain the concentration of species A close to unity and the concentration of species B close to zero.The strength of these supply and removal terms is controlled by the two parameters F and R. Note that the pattern formation process relies upon the following diffusivity constraint: D B = D A /2.One can vary F and R freely and observe the aforementioned range of stable and dynamic structures.
In the case of non-zero solvent velocity, the two equations above gain an extra advection term and become: The dynamics of the GSRDS can be reproduced by the RLBM with the following evolution equations: ) with equilibria defined by Equation (16).Such an RLBM system was explored for a still fluid by Ayodele et al. [23], under several parameter combinations, and subsequently for a sheared fluid [24].As a comprehensive assessment of the RLBM's ability in the context of the GSRDS, the entire phase space of chemical patterns can be reproduced by simulating a system in which the two parameters F and R vary as a function of the vertical and horizontal coordinate, respectively.An instantaneous snapshot from such a system is shown in Figure 2.
This phase diagram agrees well with analogous diagrams in previous works [89,91] and highlights the complex emergent properties of this simple system.As was found by Ayodele et al. [23], the characteristic size of the chemical structures can be modified by changing the diffusivity D A (keeping the ratio D B = D A /2 constant).This is achieved in the RLBM simply by changing the relaxation parameters τ A and τ B .To highlight this effect, Figure 3 shows four different phase portraits, each with a different value of τ A .
These phase diagrams highlight the coarse-graining effect of increasing the diffusivities of the two species, which affects all patterns in an essentially equivalent manner.Thus, we now see how the three key parameters of the system can be tuned to produce any desired pattern from the system's repertoire.Changes in F and R take the system through different structural phases, and changes in D A allow a simple increase or decrease of the characteristic pattern length scale.
Such observations can be leveraged for the exploration of more complex interactions in systems with large numbers of interacting chemical species.If there were mutual catalytic enhancements and inhibitions, as well as physical interactions such as modulations of viscosity and/or passive scalar diffusivities, one can imagine a vast range of potential life-like dynamics, including competition, symbiosis and arms races.

Thermal Reactive LBM
The isothermal constraint of the RLBM is somewhat limiting.A vast range of important systems involves the coupled dynamics of reactions, heat transfer and convection (e.g., Andres and Cardoso [100], Rogers and Morris [101]).One of the most common applications of thermal RLBMs (TRLBMs) has been for problems of combustion [28][29][30][77][78][79].This section will present a simple approach to relaxing the isothermal assumption of the RLBM, beginning with thermal rate constants and moving on to enthalpy changes.
Chemical reaction rates are functions of many variables (pressure, cross-section of reactant molecules, etc.), but one of the most influential is temperature.The temperature dependence of simple reactions is encapsulated in the Arrhenius equation, which expresses the rate in terms of the frequency factor A, activation energy E and temperature T: which can simply replace the constant rate parameter used to scale reactions in the isothermal RLBM.
Another crucial thermal feature of reactions is changes in enthalpy from the reactant to product state.If we assume that no pressure-volume work is done during a reaction, any free energy difference will be released as a local change in internal energy, i.e., a gain or loss of heat.In the majority of LBMs, the fluid is assumed to be incompressible, which automatically constrains the model to cases of negligible compression work from reactions.If we follow the standard convention of a negative enthalpy change implying an exothermic reaction, the internal energy distributions of a TRLBM are governed by: where R is the number of reactions, A k , E k and ∆H k are the frequency factor, activation energy and enthalpy change of reaction k, N is the number of chemical species and α kl is the stoichiometric coefficient for reaction k and chemical species l.
In addition to combustion, variants of the TRLBM have been applied to several complex fluid problems including microfluidics [102] and fuel cells [103].It has also been used to explore thermal effects in the GSRDS.The addition of thermal rate constants and finite reaction enthalpies introduced a completely new set of physical effects.A strongly endothermic reaction in a convecting fluid produced a competitive dynamic, in which chemical structures and convection plumes were observed to compete for a common heat source.This behavior was analogous to ecological effects, allowing compelling parallels to be drawn between non-living dissipative structures and organisms [76,104].
It was also found that two sets of GSRDSs embedded in the same system could exhibit spontaneous temperature regulation, if one set had an exothermic reaction and the other an endothermic reaction.The composite system was found to be robust to large changes of boundary temperature, even when such changes would normally eliminate the conditions for stable pattern formation [105].This was a further example of an effect strongly reminiscent of the living world (homeostasis), exhibited by a purely physicochemical system.

Reversibility and Heat Transport Enhancement
This final section will present new results from a study in which the TRLBM was used to address a simple thermodynamic transport problem.The system in question will be a two-dimensional fluid with two solute species undergoing a reversible reaction A B. If a steady vertical temperature gradient is imposed on the system, it will respond with a net heat flux in the same direction.This heat flux is normally constrained by the physical properties of the fluid (viscosity, thermal diffusivity and coefficient of thermal expansion) and buoyancy force (strength of gravity), as well as the strength of the temperature gradient.If that gradient is strong enough, the fluid will undergo convective motion.With the two chemical species and reaction present, the net heat flux of the system will be affected (assuming non-zero enthalpy change ∆H).We can ask whether it will be enhanced or diminished, and if so, how does this scale with the amount of chemical species present.They represent an extra degree of freedom that the system can potentially make use of.
To address this problem, a series of simulations of a single size were carried out, with two different sets of fluid parameters.The first set had a Rayleigh number of Ra = 5 × 10 3 , and the second had a value of Ra = 5 × 10 4 .For each of the two Rayleigh numbers, simulations were performed with fixed temperature BCs, and also with fixed flux BCs.The total mass of chemical species present in the system was varied for each run (but held constant during the run, there was no flux of mass in or out of the system).It was raised from zero to a value that gave mean concentrations of ψ A = ψ B = 4.Each simulation was initialized with a small amount of noise and concluded when a steady state was reached.Figure 4 displays the steady state configuration of a simulation with Ra = 5 × 10 4 , fixed temperature BCs, and a mean initial chemical species concentration of The system configured itself such that the endothermic reverse reaction was more prevalent at the lower, hotter end of the domain, leading to a higher steady state concentration ψ A in that region.The forward, exothermic reaction dominated at the cooler upper boundary, yielding higher concentrations of B near the top wall.Indeed, it appears that the system maximized its heat flux by arranging in this way.Having the endothermic reaction dominate at the lower wall means that the system invests heat energy in converting B to A; this energy is carried by A as it is advected by the flow before being released gradually as the reaction changes to A → B at the upper wall.The two chemical species and their reaction have been utilized as an additional heat delivery mechanism in parallel with the advection and diffusion of heat alone.Note, however, that the velocity field is identical to that of a purely convective, Rayleigh-Bénard system (see Figure 4a).
It is straightforward to quantify the heat transport increase exhibited by these systems.Figure 5 shows the different components of heat flux as a function of average chemical species concentration.The heat flux values are normalized by total heat flux; thus, the curves show fractional contributions to the three different modes of heat transport: diffusion, convection and advection by chemical species.Similar trends can be seen for both types of BC.With increases in chemical species concentration, an increasing heat transport function is exhibited by the chemical aspect of the system, and the relative roles of the other heat transport mechanisms are reduced.It seems the system is gradually switching over to a configuration in which advection and reaction provide an increasingly larger fraction of the total heat transport.This makes sense because higher quantities of chemical species mean that the heat exchange as a result of reactions will increase, and the amount of heat that can be invested in the enthalpy of chemical species also increases.
However, one would expect that this trend cannot continue ad infinitum.There is likely to be a lower limit for the diffusive heat transport (in the fixed boundary temperature case, the absolute diffusive heat flux is constant since the temperature difference and thermal diffusivity are constant), since all heat must enter and leave the system by diffusion.The blue curves of Figure 5 appear to have a decreasing curvature with a gradient that seems to be tending to zero.
Considering the role of convection, it should be kept in mind that the advection of chemical species relies on there being a sustained fluid flow.That fluid flow is induced by the temperature gradient, and so, in the case of fixed flux BCs, as the system becomes "more efficient" at transporting heat, it gradually diminishes the driving force for the generation of buoyant motion.Thus, for these BCs, it is likely that there will be an eventual saturation of the relative magnitude of chemical advective heat transport.
In the fixed temperature case, this can never be an issue because the driving force is constant (the boundary heat flux adjusts itself to maintain the boundary temperatures at their relevant values).One observation that is not shown in Figure 5 is that for this BC, the total heat flux increased with ψ.The presence of the chemical species enhanced the system's heat transport abilities, and it used those new abilities to augment its total heat flux above that which would occur without the action of the chemistry.

Conclusions
This paper described a powerful version of the LBM, in which fluid dynamics, heat transfer, buoyancy-driven convection, passive scalar advection, chemical reactions and enthalpy changes can all be incorporated simultaneously.The method is a natural extension of previous versions, and it retains the second order accuracy of the standard LBM.It is simple to code and easily parallelizable.
The single-phase and thermohydrodynamic forms of the model were described, and previous applications to natural convection and the validity of the principle of MEP were presented.It was shown that isothermal chemical reactions can be accurately simulated, and the complex, pattern-forming behavior of the GSRDS was comprehensively reproduced.Thermally-resolved systems were then explored, by introducing a temperature-dependent rate factor and reaction-induced internal energy changes due to enthalpy differences.Previous results including ecological phenomena (competition and temperature regulation) of dissipative structures in the thermal GSRDS were discussed.
Finally, an investigation of heat flux enhancement by the reaction and advection of chemical species was presented.It was observed that convecting fluid systems with a reversible reaction between solute species are able to transport more heat than their non-chemical counterparts.The presence of the chemical species provided a parallel channel for heat flow, the spatial structure of which was optimized by the system to maximize total heat flux.
There is an almost limitless range of systems that can be investigated using this TRLBM, since the number of chemical species and reactions can be increased until computational power becomes a binding constraint.Porous media can also be incorporated as they are in the standard LBM, and geochemical systems such as hydrothermal vents can be readily modeled.A range of fields from biophysics to astrobiology to chemical engineering could make use of this method.

Further Work
To further validate the present method, it would be useful to take a range of complex thermal chemical systems and compare their results with those of mainstream commercial multi-physics solvers, and perhaps also with laboratory experiments.It would be interesting to explore the limits at which the TRLBM begins to lose accuracy, and the computational efficiency of the method should be compared to more conventional CFD models.Accuracy limits might also be amenable to analytical calculation.

Figure 1 .
Figure 1.Numerical error of reaction diffusion simulations as a function of relaxation parameter τ A for several reaction rates k.

FFigure 2 .
Figure 2. RLBM simulation of a GSRDS, with spatially-varying parameters.The supply rate F and removal rate R are linear functions of the vertical and horizontal coordinate, respectively.The color map shows the order parameter φ(x, t) = ψ A (x, t) − ψ B (x, t), with red corresponding to values of φ = 1 and the deepest blue corresponding to φ ≈ −0.5.Also shown is the Hopf bifurcation (dotted line) and the saddle node bifurcation (solid line).

Figure 3 .
Figure 3. Phase portraits of RLBM simulations of a GSRDS with different relaxation parameters τ A .

Figure 4 .
Figure 4. Steady state flow and concentration fields of a TRLBM simulation in which chemical species A reversibly transforms into species B: A B, which has a lower enthalpy value (the reaction is exothermic in the forward direction).(a) Temperature field with velocity streamlines; (b) concentration field of the first component ψ A ; and (c) concentration field ψ B .The initial concentrations were ψ A0 = ψ B 0 = 2.

Figure 5 .
Figure 5. Heat flux components as a function of average chemical species concentration for TRLBM chemical convection simulations at two different Rayleigh numbers.The heat flux values (corresponding to diffusion in blue, convection in red and advection by passive scalar chemical species in green) are normalized by the total heat flux.