Structure preservation in high-order hybrid discretisations of potential-driven advection-diffusion: linear and nonlinear approaches

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, $\mathbf{u}$ is velocity, $\rho$ is the fluid density, $p$ is the hydrodynamic pressure, $\eta$ is the coefficient of dynamic viscosity, $\mathbf{g}$ is the gravity, $n$ and $c$ represent the concentrations of bacteria and oxygen, respectively. $\kappa$ is the oxygen consumption rate, $\chi$ is the chemotactic sensitivity, $D_n$ and $D_c$ are diffusion coefficients for bacteria and oxygen, respectively. $V_b$ and $\rho_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$\acute{e}$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 three-dimensional 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 chemotaxis-fluid 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 Ⅱ, the MRT-LBM with B-W scheme for the chemotaxis-fluid model is proposed. In Section Ⅲ, 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 Ⅳ.

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.

where $\mathbf{z}$ is the unit vector in the vertical direction. $r(c)$ is regularized by

where $c^*$ is a cut-off constant and $\epsilon>0$ is a constant close to zero.

The dimensionless parameters $\alpha$, $\beta$, $\gamma$, $\delta$, 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 (\mathbf{x}, t)$ is the discrete distribution function of particle at site $\mathbf{x}$ and time $t$ moving with speed $\hat{c}$ along the direction $\mathbf{e}_i$. $\{\bf{e}_i, i = 0, 1, \cdots, k-1\}$ is the set of discrete velocity directions, $\hat{c}$ is the sound speed, ${\Lambda _{ij}}$ is the collision matrix, $h_i^{eq}(\mathbf{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 $\varepsilon$ 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. $\mathbf{M}$ is the transformation matrix. $\mathbf{h} = \{h_i, i = 0, 1, \cdots k-1\}$. For the D2Q9 model, $\mathbf{M}$ is defined as

Through multiplying the transformation matrix $\mathbf{M}$, the collision process can be rewritten onto the moment space as

where $\mathbf{\widetilde{S}} = \mathbf{M\Lambda M^{-1}}$ is a diagonal relaxation matrix, given by

The equilibrium moment $\mathbf{m}^{eq}$ is defined as

$\widehat{\bf{H}} = \mathbf{M}\mathbf{H}$ is the corresponding force moment, which is expressed as

where $\mathbf{H} = (H_0, ..., H_8)^T$. $\mathbf{F}_t$ is the external force, which has the following form

In this simulation, $\mathbf{F}_t$ should be expressed as

The first order explicit Euler scheme is used to discrete (12) as

where $\mathbf{m}^+ = \mathbf{M} \mathbf{h}^+$ is the post-collision moments, and $\mathbf{h}^+ = (h^+_0, ..., h^+_8)^T$ is the post-collision distribution function. The non-dimensional relaxation matrix is defined as $\mathbf{S} = \delta t \mathbf{\widetilde{S}}$, which is given by

In simulations, $s_0 = s_3 = s_5$, $s_4 = s_6$ and $s_7 = s_8$.

Next, through using the second-order B-W scheme, Eq. (11) can be solved on a regular lattice with spacing $\delta x$,

where the time step $\delta t$ is related to the CFL number $A$ as $\delta t = A \delta x /c$.

The macroscopic quantities, $\rho$ and $\mathbf{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 $\tau_n$ and $\tau_c$ are the dimensionless relaxation factors, $\mathbf{\hat{c}}_i = \hat{c}\mathbf{e}_i$, $f_i ({\bf{x}}, t)$ and $g_i ({\bf{x}}, t)$ are the distribution functions of particle moving with velocity $\mathbf{\hat{c}}_i$ at position $\mathbf{x}$ and time $t$. $S_i ({\bf{x}}, t)$ and $\hat{S}_i({\bf{x}}, t)$ are source distribution functions, which are defined as

Here, $B$ is a parameter, which satisfies $\alpha = B c_s^2{\tau _n}\delta t$. $f_i^{eq} (\mathbf{x}, t)$ and $g_i^{eq} (\mathbf{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(\mathbf{x}, t)$ and oxygen $c(\mathbf{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.

(27)-(28) can be regarded as a special case of the generalized K-S model. $\mu$ and $\chi$ represent the diffusion constants. Without loss of generality, the values of these parameters are set as $\mu = \chi = 1$. The computational domain is $\Omega = [-0.5, 0.5]\times[-0.5, 0.5]$ and the initial conditions are

The boundary conditions are

where $\mathbf{v}$ is unit normal vector and $\partial \Omega$ 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\times256$ uniform grid in space and time step $\delta t = 5.0 \times 10^{-7}$. Fig. 1 shows numerical solutions of bacteria concentrations $n(\mathbf{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 $\delta t = 2.5 \times 10^{-7}$. Numerical solutions can be found in Fig. 2, which show that the positivity of bacteria concentrations $n(\mathbf{x}, t)$ is well preserved. In addition, one-dimensional profiles of $n(\mathbf{x}, t)$ along $y = 0$ at $t = 1.15\times10^{-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.

3.2.   Numerical study of the coupled NS-KS equations

In this section, the incompressible N-S equations coupled with K-S equation is solved for the convergence test of the proposed LBM with B-W scheme. The mathematical model has the following form,

In the numerical simulation, $f = \kappa = 1$. The boundary conditions are

The computational domain is $\Omega = [0, 1]\times[0, 1]$. The above problem has the following exact solutions [22],

To measure the accuracy, the following relative error is applied:

where $\phi^*$ and $\phi$ 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, \mathbf{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 $\gamma = 418$, $S_c = 7700$, $c^* = 0.3$, the grid spacing and time step are $h = 1/128$ and $\delta t = 2\times10^{-5}$, respectively, and $\epsilon = h$. The initial conditions are defined as:

At the top boundary $\Omega_{top}$, the boundary conditions are given as

where $\mathbf{u} = (u, v)$. At the bottom of the domain $\partial \Omega_{bot}$, the boundary conditions are:

The mixed boundary condition in Eq. (31) can be regrouped as $\alpha c_y = (\ln n)_y$, where $\alpha = \chi r(c)/D_n$. By integrating $\alpha 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 $\alpha = 5$, $\beta = 5$ and $\delta = 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\leq 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 $\beta$ is investigated for $\alpha = 5$ and $\delta = 1$. Numerical results for $\beta = 7.23$, 10, 20, 40 are shown in Fig. 5. When the value of $\beta$ 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 $\beta$, while the cell up-swimming increases (see Fig. 5(b)). The effect of $\alpha$ is also investigated for the fixed value of $\beta = 10$ and $\delta = 1$. The cases of $\alpha = 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 $\alpha$ 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.   Numerical study of bacterial bioconvection

In this subsection, the system (5)-(8) on the rectangular domain is numerically studied. In the following numerical simulations, the coefficients are set as $\alpha = 10$, $\beta = 10$, $\gamma = 1000$, $\delta = 5$, $S_c = 500$, and the cut-off $c^* = 0.3$. The computational domain is $\Omega = [0, 6]\times[0, 1]$. The mesh grid and time step are $h = 1/256$ and $\delta t = 10^{-5}$, respectively, and we use $\epsilon = h$. The boundary conditions (31) and (32) are applied.

3.4.1.   The case of randomly perturbed homogeneous initial data

As in [17], the following initial data is considered:

where rand() is a random number uniformly distributed in the interval [0, 1].

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 $\delta$

In this subsection, the effect of the parameter $\delta$ on evolution of the cell density is investigated. In numerical simulations, the following initial data is used.

The computational domain is $\Omega = [0, 4]\times[0, 1]$ and the values of $\delta$ are set as: $\delta = 5$, 25, 50, while the value of $D_n$ is fixed. The other parameters are $\alpha = 10$, $\beta = 10$, $\gamma = 1000$, and $S_c = 500$.

Numerical results of $\delta = 5,$ 25, and 50 are shown in Figs. 8-10, which are consistent with the results in [17,15]. For $\delta = 5$ (see Fig. 8), profiles of plumes are similar to those in Fig. 7. However, for $\delta = 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.

4.   Conclusions

An accurate and stable MRT-LBM with B-W scheme for a chemotaxis-fluid model is developed. Through introducing B-W scheme in the evolution process of the LBM, the numerical stability is improved. The numerical study of the classical K-S model shows that the proposed LBM with B-W scheme could preserve the positivity of bacteria concentrations. Then, the second order accuracy is tested by simulating the coupled NS-KS model. The nonlinear dynamics of the chemotaxis-fluid model was investigated numerically by the proposed MRT-LBM with B-W scheme. Our numerical results agree well with those in the literature under different settings. In the numerical simulation of bioconvection phenomenon, the evolution of the instability of falling plumes and the convergence towards numerically stable stationary plumes are observed.

Appendix A.   From MRT-LBE to the Hydrodynamic equations: Multi-scale Chapman-Enskog expansion

Rewritten the evolution equation (17) up to $O(\delta x^3)$, one can obtain the following equation

Multiplying Eq. (16) by inverse of the transformation matrix $M$, and substituting it into the above equation and expanding the variables around $(\mathbf{x}, t)$ up to $O(\delta x^2)$ and $O(\delta t^2)$, we obtain

where $\Omega_i = - {\Lambda _{ij}}({h_j} - h_j^{eq}) + {H_i}$.

$\partial _t^2{h_i}$ can be expressed as

Thus, Eq. (34) can be rewritten as

where $D_i = {\partial _t} + c{\mathbf{e}_i} \cdot \nabla$.

In addition, from Eq. (34), we can see that ${\Omega _i} = {D_i}{h_i} + O(\delta t)$. Thus, Eq. (36) is also equivalent to

Then we introduce the following expansions:

where $\varepsilon$ is a small parameter. With the expansions, Eq. (37) can be rewritten in consecutive orders of $\varepsilon$ as

where ${D_{1i}} = {\partial _{{t_1}}} + {{\mathbf{c}}_i} \cdot {\nabla _1}$.

Multiplying the transformation Matrix $\mathbf{M}$ on both side of Eq. (39), we can obtain the following moment equations:

where $\mathbf{D}_1 = \partial_{t_1}\mathbf{I}+\mathbf{C}_\alpha \partial_{1\alpha}$, $\mathbf{C}_\alpha$ 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+3\mathbf{u}^2$, $B_y = -1+6u^2+3\mathbf{u}^2$.

Similarly, From Eq. (40b), the scale equations of conserved quantities $\rho$, $j_x$ and $j_y$ on the $t_2$ time scale can be rewritten as

To close the hydrodynamic equations at the second order of $\varepsilon$, the terms of $e^{(1)}$, $p_{xx}^{(1)}$ and $p_{xy}^{(1)}$ in Eqs. (44) and (45) should be estimated. Under the low Mach number assumption, these terms can be evaluated as

With Eqs. (46), (47) and (48), we can obtain the hydrodynamic equations at $t_1$ and $t_2$ scales,

Continuity equations

Momentum equations

where $\nu = \rho \eta$ and $\zeta = \rho \xi$ are the kinematic and bulk viscosities, respectively. In the present MRT-LB model, we enforce $\nu = \frac{1}{3}(\frac{1}{s_7}-\frac{1}{2})\delta t$ and $\zeta = \frac{1}{3}(\frac{1}{s_1}-\frac{1}{2})\delta t$.

Combining the above equations on $t_0$ and $t_1$ scale, the following hydrodynamic equations can be obtained.