A MULTIPLE-RELAXATION-TIME LATTICE BOLTZMANN METHOD WITH BEAM-WARMING SCHEME FOR A COUPLED CHEMOTAXIS-FLUID MODEL

. In this work, a multiple-relaxation-time (MRT) lattice Boltzmann method (LBM) is proposed to solve a coupled chemotaxis-ﬂuid model. In the evolution equation of the proposed LBM, Beam-Warming (B-W) scheme is used to enhance the numerical stability. In numerical experiments, at ﬁrst, the comparison between the classical LBM and the present LBM with B-W scheme is carried out by simulating blow up phenomenon of the Keller-Segel (K-S) model. Numerical results show that the stability of the present LBM with B-W scheme is better than the classical one. Then, the second order convergence rate of the proposed B-W scheme is veriﬁed in the numerical study of the coupled Navier-Stokes (N-S) K-S model. Finally, through solving the coupled chemotaxis-ﬂuid model, the formation of falling bacterial plumes is numerically investigated. Numerical results agree well with existing ones in the literature.

1. Introduction. The chemotaxis refers to the collective motion of bacteria, cells or an organisms in response to an attractant gradient. There are numerous examples of chemotaxis in animal and insect ecology, biological and biomedical sciences [9]. Deep understanding of the behavior of chemotaxis phenomenon is of great significance to the biological and medical applications.
Various kinds of mathematical models have been developed by both experimentalists and theoreticians to describe the chemotactic phenomenon. In addition, to further describe and interpret the bioconvection phenomenon [14], a variety of coupled chemotaxis-fluid models have been proposed and studied in [5,8,24]. The coupled model is expressed as follows: Here, u is velocity, ρ is the fluid density, p is the hydrodynamic pressure, η is the coefficient of dynamic viscosity, g is the gravity, n and c represent the concentrations of bacteria and oxygen, respectively. κ is the oxygen consumption rate, χ is the chemotactic sensitivity, D n and D c are diffusion coefficients for bacteria and oxygen, respectively. V b and ρ b are the bacterial volume and density, respectively. r(c) is the dimensionless cut-off function. It can be seen from the above system that the chemotaxis-fluid model involves highly nonlinearity coupled with chemotaxis, diffusion and convection, which brings great difficulties to get the analytical solution. Recently, several numerical methods have been developed to study the chemotaxis phenomenon, such as finite element, finite volume schemes [13,21], upwind-difference potentials method [10], fractional step algorithms [20,25], interior penalty/discontinuous Galerkin methods [11,12], and cell-overcrowding prevention models [1,7], etc. In the work by Chertock et al. [5], a coupled chemotaxis-fluid system has been numerically studied by a high-resolution vorticity-based hybrid finite-volume finite-difference scheme. In [22], the chemotaxis phenomenon in incompressible viscous fluid has been numerically studied by a compact finite difference scheme. The plume patterns were studied by solving a chemotaxis-diffusion-convection coupling system in [6], where the comparison of different systems of chemotaxis-diffusion-convection, Rayleigh-Bénard convection and the double diffusive convection were investigated. And it was indicated that although physical mechanisms are different, the convection patterns and the dimensionless differential systems are similar. In addition, a threedimensional chemotaxis-fluid model was numerically investigated in [17]. The formation of falling bacterial plumes and its convergence to a steady state have been observed.
In this work, we will develop an efficient LBM to solve the coupled chemotaxisfluid model. LBM originates from kinetic Boltzmann equation, which has been successfully applied in modeling complex fluid flows [18,4,2]. In addition, LBM also shows its excellent capability in solving the nonlinear systems [23,3,28,26]. For the chemotaxis problem, a simplified LBM was developed for a bacterial chemotaxis model by Hilpert [16]. While, Yan et al. [27] developed a MRT-LBM to investigate the traveling bacterial bands of chemotactic bacteria. However, due to the highly nonlinearity of the coupled chemotaxis-fluid model, we need to develop more accurate and stable LBM to solve this complex problem. In [19], a novel MRT-LBM was proposed to simulate multi-phase fluids with Peng-Robinson (P-R) equation of state. To enhance the stability of the MRT-LBM, B-W scheme is introduced in the discretization of the evolution equation. Numerical results showed that the P-R free energy model was solved precisely and spurious currents were eliminated effectively. Inspired by this, we will apply the MRT-LBM with B-W scheme to solve the coupled chemotaxis-fluid model. The rest of this article is organized as follows. In Section II, the MRT-LBM with B-W scheme for the chemotaxis-fluid model is proposed. In Section III, several numerical experiments are carried out on different chemotaxis problems, including the classical K-S model, the coupled Navier-Stokes-Keller-Segel (NS-KS) model, and the bacterial biocovection problem. The paper ends with some conclusions in Section IV.
2. Numerical method. In this section, we will introduce the MRT-LBM with B-W scheme for the coupled chemotaxis-fluid model (1)- (4). To better investigate the coupled chemotaxis-fluid model, the following dimensionless form is taken: where n r is the characteristic cell density, c air is the oxygen concentration in the air, and L is a characteristic length. Through introducing the above scaled variables and dropping the prime notation, the following system can be obtained.
∇ · u = 0, where z is the unit vector in the vertical direction. r(c) is regularized by where c * is a cut-off constant and > 0 is a constant close to zero. The dimensionless parameters α, β, γ, δ, and the Schmidt number S c are given below 2.1. Lattice Boltzmann method for the incompressible Navier-Stokes equations. In this section, LBM with MRT collision operator is applied to solve the incompressible N-S equations (7)- (8). To capture the chemotaxis phenomenon more accurately, B-W scheme is introduced in the discretization of the evolution equation of MRT-LBM. The discrete velocity Boltzmann equation with MRT collision operator is expressed as where, h i (x, t) is the discrete distribution function of particle at site x and time t moving with speedĉ along the direction e i . {e i , i = 0, 1, · · · , k − 1} is the set of discrete velocity directions,ĉ is the sound speed, Λ ij is the collision matrix, h eq i (x, t) is the equilibrium distribution function (EDF), and H i is the force distribution function.
The discrete velocity Boltzmann equation (9) is solved by a time-splitting scheme, which is decomposed into two sub-processes, i.e., the collision process, and the streaming process, In the MRT model, the collision subprocess is carried out in the moment space. We take the generally used D2Q9 (two-dimensional space with nine discrete velocities) model as an example. The distribution functions h i in the moment space are defined as where e and ε are related to total energy and the function of energy square; j x and j y are relevant to the momentum; q x and q y are related to the x and y components of the energy; p xx and p xy are the corresponding diagonal and off-diagonal components of the stress tensor, respectively. M is the transformation matrix.
Through multiplying the transformation matrix M, the collision process can be rewritten onto the moment space as where S = MΛM −1 is a diagonal relaxation matrix, given by The equilibrium moment m eq is defined as H = MH is the corresponding force moment, which is expressed as where H = (H 0 , ..., H 8 ) T . F t is the external force, which has the following form In this simulation, F t should be expressed as The first order explicit Euler scheme is used to discrete (12) as where m + = Mh + is the post-collision moments, and h + = (h + 0 , ..., h + 8 ) T is the post-collision distribution function. The non-dimensional relaxation matrix is defined as S = δt S, which is given by S = diag{s 0 , s 1 , s 2 , s 3 , s 4 , s 5 , s 6 , s 7 , s 8 }.
Next, through using the second-order B-W scheme, Eq. (11) can be solved on a regular lattice with spacing δx, where the time step δt is related to the CFL number A as δt = Aδx/c. The macroscopic quantities, ρ and u are calculated as 2.2. Lattice Boltzmann method for the transport equation of concentrations of bacteria and oxygen. The evolution equations of the LBM for (5) and (6) are given by where τ n and τ c are the dimensionless relaxation factors,ĉ i =ĉe i , f i (x, t) and g i (x, t) are the distribution functions of particle moving with velocityĉ i at position x and time t. S i (x, t) andŜ i (x, t) are source distribution functions, which are defined as Here, B is a parameter, which satisfies α = Bc 2 s τ n δt. f eq i (x, t) and g eq i (x, t) are local equilibrium distribution functions, and read as ). (19) and (20) are decomposed into two sub-processes: the collision process and the streaming process, The first order explicit Euler scheme is used to solve the collision process (21)- (22), which leads to B-W scheme is used to solve (23) and (24), The concentrations of the bacteria n(x, t) and oxygen c(x, t) are computed by 3. Numerical experiments. To validate that the stability is improved by the present LBM with B-W scheme, and also to verify the capability of the proposed method, a series of chemotaxis problems are numerically investigated. In the following numerical experiments, the general bounce-back scheme in reference [29] is used to treat the macroscopic boundary conditions in the LBM with B-W scheme. 3.1. Validation of the LBM with B-W scheme. In this section, the following classical K-S chemotaxis problem is considered. The boundary conditions are where v is unit normal vector and ∂Ω represents the domain boundary. Under above initial and boundary conditions, the solution of (27)- (28) will blow up at the origin in finite time, which brings a great challenge in the numerical simulation. To illustrate that the stability is improved by B-W scheme, the comparison between standard LBM and the present LBM with B-W scheme is carried out. In the numerical experiment, the D2Q9 lattice model is used with a 256 × 256 uniform grid in space and time step δt = 5.0 × 10 −7 . Fig. 1 shows numerical solutions of bacteria concentrations n(x, t) at different time points by the standard LBM, i.e., A = 0. It can be observed that numerical solutions appear negative values as time evolves. Then, the LBM with B-W scheme is applied to solve the above initial-value problem. In the simulation, the same grid is used with the CFL number A = 0.5, which gives δt = 2.5×10 −7 . Numerical solutions can be found in Fig. 2, which show that the positivity of bacteria concentrations n(x, t) is well preserved. In addition, one-dimensional profiles of n(x, t) along y = 0 at t = 1.15×10 −4 are depicted in Fig.  3. Negative values can be observed clearly near the origin in Fig. 3(a), while not in Fig. 3(b), which verifies that the numerical stability is improved by the proposed LBM with B-W scheme.
To measure the accuracy, the following relative error is applied: where φ * and φ represent the analytical and numerical solutions, respectively, and the summation is taken over all grid points. It can be seen from Table 1 that numerical solutions of n, c, u by the proposed LBM with B-W scheme has second order convergence.

3.3.
Numerical study of the coupled chemotaxis-fluid system. The chemotaxis response of bacterial suspensions can be described by the chemotaxis-fluid model (5)- (8). The balance between chemotaxis, diffusion, and convection of bacteria leads to the formation and stability of plumes. There still needs to a deep understanding of the particular impact of each mechanism. In the simulations, the parameters of (5)-(8) are set as γ = 418, S c = 7700, c * = 0.3, the grid spacing and time step are h = 1/128 and δt = 2 × 10 −5 , respectively, and = h. The initial conditions are defined as: At the top boundary Ω top , the boundary conditions are given as where u = (u, v). At the bottom of the domain ∂Ω bot , the boundary conditions are: The mixed boundary condition in Eq. (31) can be regrouped as αc y = (ln n) y , where α = χr(c)/D n . By integrating αc y = (ln n) y with respect to y, the mixed boundary condition is transformed into Dirichlet boundary condition, then the general bounce-back scheme in reference [29] is used.
First of all, as the numerical experiments in [17], the case of α = 5, β = 5 and δ = 0.25 is considered. Numerical results of the vertical profiles of c and cell densities n at t=0.22 are shown in Fig. 4(a), which agree well with the available data in [17]. As a result of the cut-off of the chemotactic convection for oxygen levels below c ≤ c * , an increase of cells towards the bottom can be found from the vertical profile of the cell density. The streamlines in Fig. 4(b) shows that the bioconvection phenomenon is well captured. Next, the effect of β is investigated for α = 5 and δ = 1. Numerical results for β = 7.23, 10, 20, 40 are shown in Fig.  5. When the value of β increases, the oxygen consumption increases relative to the oxygen diffusion. It can be seen from Fig. 5(a) that the oxygen density decreases with the increases of the value of β, while the cell up-swimming increases (see Fig.  5(b)). The effect of α is also investigated for the fixed value of β = 10 and δ = 1. The cases of α = 1, 2, 4, 5.952 are considered. Numerical results are shown in Fig.  6, which indicate that the cell density near the surface increases when the value of α increases (see Fig. 6(b)), and less overall oxygen consumption occurs in these regions (see Fig. 6(a)). Similar results as in Fig. 5 and Fig. 6 can be found in [17].
3.4.1. The case of randomly perturbed homogeneous initial data. As in [17], the following initial data is considered: n 0 (x, y) = 0.8 + 0.2rand(), c 0 (x, y) = 1, u 0 (x, y) = 0, The evolution of the cell density n at different time is depicted in Fig. 7. It can be observed that the cell concentration increases with time as a result of the cells accumulating at the top of the chamber, while the distinct phenomenon of the cell concentration can be found at the bottom of the chamber. As a consequence, several instabilities at a high-concentration layer is triggered below the fluid-air surface at t = 0.32. At t = 0.37, two falling plumes is formulated due to the instabilities and all plumes hit the bottom of the chamber at t = 0.42. At the bottom of the plumes, the cells move onto the edge of the plumes, and are advected by the flow field near the fluid-air surface. Thus the smaller plumes and thicker cell-rich layer at t = 0.47 than those at t = 0.42 can be observed. The velocity components u(x, y, t) and v(x, y, t) are depicted as arrows in Fig. 7. It can be seen that the fluid velocity is downward in regions of large concentrations. The above phenomenon was also observed in [5,17,6].

3.4.2.
Effect of increased value of δ. In this subsection, the effect of the parameter δ on evolution of the cell density is investigated. In numerical simulations, the following initial data is used. Numerical results of δ = 5, 25, and 50 are shown in Figs. 8-10, which are consistent with the results in [17,15]. For δ = 5 (see Fig. 8), profiles of plumes are similar to those in Fig. 7. However, for δ = 25 and 50 (see Figs. 9 and 10, respectively), profiles of plumes are significantly different, where we could observe that the two instabilities are approaching each other and merge into one instability.  and O(δt 2 ), we obtain where Thus, Eq. (34) can be rewritten as where D i = ∂ t + ce i · ∇. In addition, from Eq. (34), we can see that Ω i = D i h i + O(δt). Thus, Eq. (36) is also equivalent to Then we introduce the following expansions: where ε is a small parameter. With the expansions, Eq. (37) can be rewritten in consecutive orders of ε as where Multiplying the transformation Matrix M on both side of Eq. (39), we can obtain the following moment equations: where D 1 = ∂ t1 I + C α ∂ 1α , C α is the discrete velocity matrix. In addition, from Eqs. (18) and (40a), we derive On the t 1 time scale, Eq. (40b) can be rewritten as follows: where B x = −1 + 6v 2 + 3u 2 , B y = −1 + 6u 2 + 3u 2 .