Discrete effects on the source term for the lattice Boltzmann modelling of one-dimensional reaction–diffusion equations

This work presents a detailed numerical analysis of one-dimensional, time-dependent (linear) reaction–diffusion type equations modelled with the lattice Boltzmann method (LBM), using the two-relaxation-time (TRT) scheme, for the D1Q3 lattice. The interest behind this study is twofold. First, because it applies to the description of many engineering problems, such as the mass transport in membranes, the heat conduction in fins, or the population growth in biological systems. Second, because this study also permits understanding the general effect of solution-dependent sources in LBM, where this problem offers a simple, yet non-trivial, canonical groundwork. Without recurring to perturbative techniques, such as the Chapman-Enskog expansion, we exactly derive the macroscopic numerical scheme that is solved by the LBM-TRT model with a solution-dependent source and show that it obeys a four-level explicit finite difference structure. In the steady-state limit, this scheme reduces to a second-order finite difference approximation of the stationary reaction–diffusion equation that, due to artefacts from the source term discretization, may operate with an effective diffusion coefficient of negative value, although still remaining stable. Such a surprising result is demonstrated through an exact stability analysis that proves the unconditional stability of the LBM-TRT model with a solution-dependent source, in line with the already proven source-less pure diffusion case (Lin et al., 2021). This proof enlarges the confidence over the LBM-TRT model robustness also for the (linear) reaction–diffusion problem class. Finally, a truncation error analysis is performed to disclose the structure of the leading order errors. From this knowledge, two strategies are proposed to improve the scheme accuracy from second- to fourth-order. One exclusively based on the tuning of the LBM-TRT scheme free-parameters, namely the two relaxation rates and the lattice weight coefficient, and the other based on the redefinition of the structure of the relaxation rates, where the leading order truncation error is absorbed into one of the relaxation rates, liberating the other to improve additional features of the scheme. Numerical tests presented in the last part of the work support the ensemble of theoretical findings.