A lattice Boltzmann model for the coupled cross-diffusion-fluid system
aa r X i v : . [ phy s i c s . f l u - dyn ] J a n A lattice Boltzmann model for the coupled cross-di ff usion-fluidsystem Chengjie Zhan a , Zhenhua Chai a,b , Baochang Shi a,b, ∗ a School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan, 430074, China b Hubei Key Laboratory of Engineering Modeling and Scientific Computing, Huazhong University of Science andTechnology, Wuhan 430074, China
Abstract
In this paper, we propose a lattice Boltzmann (LB) model for the generalized coupled cross-di ff usion-fluid system. Through the direct Taylor expansion method, the proposed LB model cancorrectly recover the macroscopic equations. The cross di ff usion terms in the coupled system aremodeled by introducing additional collision operators, which can be used to avoid special treat-ments for the gradient terms. In addition, the auxiliary source terms are constructed properlysuch that the numerical di ff usion caused by the convection can be eliminated. We adopt the de-veloped LB model to study two important systems, i.e., the coupled chemotaxis-fluid system andthe double-di ff usive convection system with Soret and Dufour e ff ects. We first test the presentLB model through considering a steady-state case of coupled chemotaxis-fluid system, then weanalyze the influences of some physical parameters on the formation of sinking plumes. Finally,the double-di ff usive natural convection system with Soret and Dufour e ff ects is also studied, andthe numerical results agree well with some previous works. Keywords:
Cross-di ff usion-fluid system, lattice Boltzmann model, chemotaxis-fluid,double-di ff usive natural convection, Soret and Dufour e ff ects
1. Introduction
Cross-di ff usion systems have been widely used to describe multi-species interaction in manyfields, for instance, chemotactic cell migration[1, 2], population dynamic in biological systems[3],pedestrian dynamics[4], multicomponent di ff usion [5] and so on. Such systems can also inducesome interesting and multiple phenomena, which are significant in both science and engineer-ing. However, most of these systems usually couple with the fluid field, and one also needs toconsider the e ff ect of fluid flow. In this work, we consider a general cross-di ff usion-fluid (CDF)system in d dimensional space, ∂φ α ∂ t + ∇ · φ α u = ∇ · D αβ ∇ φ β + S α , (1a) ∇ · u = , (1b) ∗ Corresponding author.
Email address: [email protected] (Baochang Shi )
Preprint submitted to Elsevier January 11, 2021 u ∂ t + ∇ · uu = −∇ p + ∇ · ν ∇ u + F , (1c)where Φ = ( φ α ) and S = ( S α ) are m × m is the number of scalar variables. D = ( D αβ )is a m × m di ff usive matrix, and D αβ ∇ φ β means P m β = D αβ ∇ φ β with α = , , · · · , m . u and F are the d dimensional fluid velocity and external force, respectively. p is the pressure and ν is the viscosity. We note that there are two important cases of such a general system, i.e., thecoupled chemotaxis-fluid (CF) system [6–8] and the double-di ff usive convection (DDC) systemwith Soret and Dufour e ff ects [9, 10], which can be obtained through taking the specific formsof Φ , D and F .The CDF system is nonlinear and coupled, and usually it is di ffi cult to obtain its analyticalsolution. For this reason, some numerical methods have been developed for this type of system[11–15]. For CF system, Chertock et al. [16] developed a high-resolution vorticity-based hybridfinite-volume finite-di ff erence scheme to understand the interplay of gravity and chemotaxis inthe formation of two-dimensional plumes. Sheu and Chiang [17] proposed a combined com-pact di ff erence scheme of fifth-order spatial accuracy in a three-point grid stencil to investigatethe flow convection and di ff usion e ff ects on the distributions of the bacteria and oxygen con-centrations. To avoid a strong restriction on the time step, Lee and Kim [18] used an operatorsplitting-type Navier-Stokes solver to study the nonlinear dynamics of a three-dimensional CFsystem. Recently, an upwind finite element method was developed to investigate the pattern for-mation and hydrodynamical stability of the CF system, and through a comparison of chemotaxis-di ff usion, double di ff usive, and Rayleigh-B´enard convection, some similarities among them wereobtained in [19]. For DDC system with Soret and Dufour e ff ects, Nithyadevi and Yang [20], us-ing SIMPLE algorithm with QUICK scheme, studied the DDC of water with Soret and Dufourfactors in a partially heated enclosure around the density maximum. B´eg et al. [21] analyzed thesteady free convection heat and mass transfer from a spherical body in a micropolar fluid withSoret and Dufour e ff ects by the Keller-box implicit finite di ff erence method. The free convectionboundary layer flow over an arbitrarily inclined heated plate in a porous medium with Soret andDufour e ff ects was studied by transforming the governing equations into a universal form in [22].Wang et al. [23, 24] studied the oscillation of double-di ff usive B´enard convection with Soret andDufour e ff ects in a horizontal cavity by the SIMPLE algorithm. In this work, we will considerthe lattice Boltzmann (LB) method for its advantages in the study of complex physical systems.The LB method, as a mesoscopic numerical approach, has achieved great success in model-ing complex flows [25–30] and nonlinear systems, such as reaction-di ff usion equation [31, 32],convection-di ff usion equation [33–35] and so on. In addition, the LB method can describe thecoupling interaction among di ff erent physical fields well for its mesoscopic kinetic background[36, 37], and it is also suitable to study coupled problems like the present CDF system. How-ever, there is a very di ffi cult issue in the study of the CDF system in the framework of the LBmethod, i.e., how to treat the cross-di ff usion terms? Hilpert [38] presented a strategy to put thecross-di ff usion term into the equilibrium distribution function, but the finite-di ff erence scheme isneeded to compute the gradient operator. Yu et al. [39] handled the cross-di ff usion terms throughredesigning the second-moments of equilibrium distribution functions, thus no gradient opera-tors included. The cross-di ff usion term was put into the evolution equation as a source term suchthat the gradient operator can be computed by a local computational scheme in Refs. [40, 41].Ren and Chan [42] diagonalized the coupling di ff usivities matrix, thus the coupled system canbe transformed to the uncoupled convection-di ff usion equations.In this work, we will develop a LB model for the CDF system. Inspired by Refs.[43, 44]where the source terms and fluid field have not been included, to avoid the special treatments2or the gradient terms, as mentioned above, the cross-di ff usion terms in the coupled system aremodeled by some extra collision operators, which are added in the evolution of LB model.This paper is organized as follows. In Section 2, the coupled LB model for general CDF sys-tem is proposed. Through the direct Taylor expansion, the governing equations can be correctlyrecovered from present LB model. In Section 3, we simulate the coupled CF system and theDDC system with Soret and Dufour e ff ects. For the former system, after testing the LB model,we analyze the influences of some parameters on the formation of plume structures. We also con-sider the e ff ects of some physical parameters for the latter problem. Finally, some conclusionsare summarized in Section 4.
2. The coupled lattice Boltzmann model
In this section we will develop a LB model for the CDF system (1), and write the evolutionequations of present LB model as f i ,α ( x + c i ∆ t , t + ∆ t ) = f i ,α ( x , t ) − ω αβ (cid:16) f i ,β ( x , t ) − f eqi ,β ( x , t ) (cid:17) + ∆ tS i ,α ( x , t ) + ∆ tG i ,α ( x , t ) + ∆ t D i ,αβ S i ,β ( x , t ) , (2a) h i ( x + c i ∆ t , t + ∆ t ) = h i ( x , t ) − ω (cid:16) h i ( x , t ) − h eqi ( x , t ) (cid:17) + ∆ t (1 − ω F i ( x , t ) , (2b)where f i ,α ( x , t ) ( i = , , · · · , q − q represents the number of discrete velocity directions) and h i ( x , t ) are the distribution functions of scalar variable φ α and fluid field at position x and time t , f eqi α ( x , t ) and h eqi ( x , t ) are the corresponding equilibrium distribution functions. c i is the discretevelocity, ∆ t is the time step. ( ω αβ ) is a invertible m × m matrix, ω αβ and ω are relaxation factors.¯ D i ,αβ = δ αβ ∂ t + γ αβ c i · ∇ with ( δ αβ ) and ( γ αβ ) representing m × m identity and diagonal matrices,respectively.To obtain the macroscopic equation (1), the equilibrium distribution functions are defined by f eqi ,α = W i " φ α + c i · φ α u c s + [( λ αβ φ β − φ α ) c s I + ϑφ α uu ] : ( c i c i − c s I )2 c s , (3a) h eqi = σ i + W i " c i · u c s + uu : ( c i c i − c s I )2 c s , (3b)where W i is the weight coe ffi cient, and c s is the lattice sound speed in LB method. ( λ αβ ) is a m × m invertible matrix that can be used to adjust the relaxation matrix ( ω αβ ), ϑ is an adjustableparameter. σ = ( W − p / c s + ρ with ρ being a constant, σ i = W i p / c s ( i ,
0) [45].The force term F i ( x , t ) is given by F i = W i " c i · F c s + ϕ ( uF + Fu ) : ( c i c i − c s I )2 c s , (4)where ϕ is another adjustable parameter. Here the source term S i ,α ( x , t ) and auxiliary source term G i ,α ( x , t ) to be determined later.The macroscopic scalar variable, fluid velocity and pressure can be computed by φ α = X i f i ,α , (5a)3 = X i c i h i + ∆ t F , (5b) p = c s − W X i , h i − W u · u c s + ϕ ∆ t ( 1 ω −
12 ) u · F c s , (5c)where the details on the computation of pressure can be found in Appendix A.We noted that the present LB model can be used for d dimensional problems. Here somecommonly used lattice velocity models are listed below:D1Q3: c = ˆ c [0 , , − , (6a) c s = ˆ c / , W = / , W − = / , (6b)D2Q5: c = ˆ c " − − , (7a) c s = ˆ c / , W = / , W − = / , (7b)D2Q9: c = ˆ c " − − − − − − , (8a) c s = ˆ c / , W = / , W − = / , W − = / , (8b)D3Q7: c = ˆ c − − − , (9a) c s = ˆ c / , W = / , W − = / , (9b)D3Q15: c = ˆ c − − − − − − − − − −
10 0 0 1 0 0 − − − − − , (10a) c s = ˆ c / , W = / , W − = / , W − = / , (10b)D3Q19: c = ˆ c − − − − − − − − − − −
10 0 0 0 1 − − − − , (11a) c s = ˆ c / , W = / , W − = / , W − = / , (11b)where ˆ c = ∆ x / ∆ t is the particle speed, and ∆ x is the lattice spacing.4 .1. The direct Taylor expansion of present lattice Boltzmann model We now perform a direct Taylor expansion analysis of LB model for convection cross-di ff usion equation, while the derivation process for incompressible Navier-Stokes equations isshown in Appendix A. In the direct Taylor expansion [46–49], the time step ∆ t is used as thesmall expansion parameter.Applying the Taylor expansion to Eq. (2a) and based on f i ,α = f eqi ,α + f nei ,α ( f nei ,α is the non-equilibrium part of distribution function f i ,α ), we have [49] N X l = ∆ t l l ! D li f i ,α + O ( ∆ t N + ) = − ω αβ f nei ,β + ∆ tS i ,α + ∆ tG i ,α + ∆ t D i ,αβ S i ,β , (12)where D i = ∂ t + c i · ∇ . From above equation one can obtain f nei ,α = O ( ∆ t ) , (13a) N − X l = ∆ t l l ! D li ( f eqi ,α + f nei ,α ) + ∆ t N N ! D Ni f eqi ,α = − ω αβ f nei ,β + ∆ tS i ,α + ∆ tG i ,α + ∆ t D i ,αβ S i ,β + O ( ∆ t N + ) . (13b)According to Eq. (13), we can derive the equations at di ff erent orders of ∆ t , D i f eqi ,α = − ω αβ ∆ t f nei ,β + S i ,α + G i ,α + O ( ∆ t ) , (14a) D i ( f eqi ,α + f nei ,α ) + ∆ t D i f eqi ,α = − ω αβ ∆ t f nei ,β + S i ,α + G i ,α + ∆ t D i ,αβ S i ,β + O ( ∆ t ) . (14b)From Eq. (14a), we can get ∆ t D i f eqi ,α = − D i ω αβ f nei ,β + ∆ t D i ( S i ,α + G i ,α ) + O ( ∆ t ) . (15)Substituting Eq. (15) into Eq. (14b) yields D i f eqi ,α + D i ( δ αβ − ω αβ f nei ,β + ∆ t D i ( S i ,α + G i ,α ) = − ω αβ ∆ t f nei ,β + S i ,α + G i ,α + ∆ t D i ,αβ S i ,β + O ( ∆ t ) . (16)To recover the macroscopic equation (1a), the distribution functions f i ,α , f eqi ,α , S i ,α and G i ,α shouldsatisfy following conditions, X i f i α = X i f eqi α = φ α , X i c i f eqi α = φ α u , X i c i c i f eqi α = λ αβ φ β c s I + ϑφ α uu , (17a) X i S i ,α = S α , X i c i S i ,α = M S ,α , (17b) X i G i ,α = , X i c i G i ,α = M G ,α , (17c)5hich can also be used to derive the following equation, X i f nei ,α = X i f i ,α − X i f eqi ,α = . (18)Summing Eqs. (14a) and (16) over i and using above relations, one can obtain ∂ t φ α + ∇ · φ α u = S α + O ( ∆ t ) , (19a) ∂ t φ α + ∇ · φ α u + ∇ · ( δ αθ − ω αθ X i c i f nei ,θ = S α + ∆ t ∇ · h ( γ αβ − δ αβ ) M S ,β − M G ,α i + O ( ∆ t ) , (19b)where the term P i c i f nei ,θ can be derived from Eq. (14a), X i c i f nei ,θ = − ω − θβ ∆ t X i c i h D i f eqi ,β − S i ,β − G i ,β i + O ( ∆ t ) = − ω − θβ ∆ t h ∂ t φ β u + ∇ · ϑφ β uu + c s λ βκ ∇ φ κ − M S ,β − M G ,β i + O ( ∆ t ) , (20)where ( ω − θβ ) is defined as the inverse matrix of ( ω θβ ), i.e., P θ ω αθ ω − θβ = δ αβ .Substituting Eq. (20) into Eq. (19b), we get ∂ t φ α + ∇ · φ α u = ∇ · ( ω − αβ − δ αβ λ βθ c s ∆ t ∇ φ θ + S α + ∆ t ∇ · RH α + O ( ∆ t ) , (21)where RH α = (2 ω − αβ − δ αβ ) (cid:16) ∂ t φ β u + ∇ · ϑφ β uu (cid:17) + ( γ αβ − ω − αβ ) M S ,β − ω − αβ M G ,β . (22)If RH α = or RH α = O ( ∆ t ), we can correctly recover the convection cross-di ff usion equation(1a) at the order of ∆ t , ∂ t φ α + ∇ · φ α u = ∇ · D αβ ∇ φ β + S α + O ( ∆ t ) , (23)where D αβ = ( ω − αθ − δ αθ λ θβ c s ∆ t . (24)When taking M S ,α =
0, the term S i ,α can be defined by S i ,α = W i S α , (25)and ( γ αβ ) can be set as zero matrix for simplicity, thus ¯ D i ,αβ = δ αβ ∂ t , and a first-order explicitfinite-di ff erence scheme can be applied, i.e., ∂ t S i ,α ( x , t ) = [ S i ,α ( x , t ) − S i ,α ( x , t − ∆ t )] / ∆ t . In thiscase, the term RH α can be written as RH α = (2 ω − αβ − δ αβ ) (cid:16) ∂ t φ β u + ∇ · ϑφ β uu (cid:17) − ω − αβ M G ,β . (26)In the following, two cases are considered according to whether the flow field is coupled or not.6 ase 1: If we only consider the convection cross-di ff usion equation (1a) without including thefluid field, one can set ϑ = λ αβ = δ αβ , thus the equilibrium distribution function f eqi ,α is linear,and DdQ (2 d +
1) velocity models can be used for simplicity. Under condition of RH α =
0, weget M G ,α = ( δ αβ − ω αβ ∂ t φ β u . (27)In this case the term G i ,α can be given by G i ,α = ( δ αβ − ω αβ W i c i · ∂ t φ β u c s , (28)where the time derivative can be discretized by a first-order explicit finite-di ff erence scheme [50]. Case 2:
When the incompressible fluid field is considered, the value of ϑ can be arbitrary since uu in equilibrium distribution function f eqi ,α is the order of O ( Ma ), and can be neglected, thus f eqi ,α can also become linear with a special value λ αβ = δ αβ . In addition, if we take ϑ =
1, and with thehelp of Eqs. (19a) and (A.6b), RH α can be written as RH α = (2 ω − αβ − δ αβ ) h u ∂ t φ β + φ β ∂ t u + u ∇ · φ β u + φ β ∇ · uu i − ω − αβ M G ,β = (2 ω − αβ − δ αβ ) h u S β + φ β ( F − ∇ p ) i − ω − αβ M G ,β + O ( ∆ t ) . (29)It is clear that the term RH α is the order of ∆ t , if M G ,α is taken as M G ,α = ( δ αβ − ω αβ h u S β + φ β ( F − ∇ p ) i . (30)Under the incompressible condition, ∇ p is also the order of O ( Ma ). Keep this in mind, the term G i ,α can also be given by G i ,α = ( δ αβ − ω αβ W i c i · ( u S β + φ β F ) c s . (31)From above discussion, one can find that there are no special treatments needed for the cross-di ff usion terms in the present LB model, and the discretization of time derivative in auxiliarysource term also can be avoided when the fluid field is considered.
3. Numerical experiments
In this section, we will simulate the CF system and DDC system with Soret and Dufour e ff ectsin two-dimensional space. we first test proposed LB model, and then discuss the e ff ects of somephysical parameters in this two systems. In the simulations below, we set ( λ αβ ) as identity matrixand ϕ = ff usion equationsand the D2Q9 model for the Navier-Stokes equations. The relaxation factor ω αβ is computedby Eq. (24), the anti-bounce-back scheme [51] is used to treat Dirichlet boundary conditions ofconvection cross-di ff usion equations, and the half-way bounce-back scheme [52, 53] is appliedfor the other no-flux and velocity boundary conditions.7 .1. Numerical results and discussion on the chemotaxis-fluid system From system (1), we can obtain the coupled CF system through taking Φ = ( n , c ) T , D = " D n − µ r ( c ) n D c , S = (0 , − κ r ( c ) n ) T , F = g nV b ( ρ b − ρ ) /ρ and p = p /ρ , n t + ∇ · n u = ∇ · [ D n ∇ n − µ r ( c ) n ∇ c ] , (32a) c t + ∇ · c u = ∇ · D c ∇ c − κ r ( c ) n , (32b) ∇ · u = , (32c) u t + u · ∇ u = −∇ p /ρ + ∇ · ν ∇ u + g nV b ( ρ b − ρ ) /ρ , (32d)where n and c are the concentrations of bacteria and oxygen, respectively. µ is the chemotacticsensitivity, κ is the consumption rate of oxygen, D n and D c are the di ff usion coe ffi cients ofbacteria and oxygen. ρ is the pure fluid density, g is the gravitation acceleration, V b and ρ b are the volume and density of bacteria, respectively. r ( c ) is a dimensionless truncated function,which is relevant to the chemotaxis cut-o ff value c ∗ .Consider this system in two-dimensional rectangular domain Ω = [ − a , a ] × [0 , d ], the bound-ary conditions are given as follows. The top of the domain ∂ Ω top is fluid-air surface where theflux of bacteria is zero, the value of oxygen concentration equals to the air oxygen concentration c air and no fluid flows through the fluid-air surface,[ D n ∇ n − µ r ( c ) n ∇ c ] · ˆ n = , c = c air , ∂ u ∂ y = , v = , ∀ x ∈ ∂ Ω top , (33)where ˆ n is the unit outer normal vector, c air is the air oxygen concentration, x = ( x , y ) and u = ( u , v ). At the bottom of the domain ∂ Ω bot , the fluxes of bacteria and oxygen, and the fluidvelocity are zero, ∇ n · ˆ n = ∇ c · ˆ n = , u = , ∀ x ∈ ∂ Ω bot . (34)Periodic boundary condition is used for the left and right sides of the domain.Before performing any simulations, we first rewrite the system (32) in a dimensionless formthrough introducing the following variables [8, 16]: x ′ = x L , t ′ = D n L t , c ′ = cc air , n ′ = nn r , u ′ = LD n u , p ′ = L νρ D n p , (35)where L and n r are characteristic length and characteristic bacteria density, respectively. Afterdropping the prime notation in the rescaled variables, we can obtain the following dimensionlesssystem: n t + ∇ · n u = ∇ · [ ∇ n − α r ( c ) n ∇ c ] , (36a) c t + ∇ · c u = ∇ · δ ∇ c − β r ( c ) n , (36b) ∇ · u = , (36c) u t + ∇ · uu = − S c ∇ p + ∇ · S c ∇ u − S c γ n ˆ z , (36d)where ˆ z is the upwards unit vector, α , β , δ , γ and Schmidt number S c are dimensionless param-eters, which are defined as α = µ c air D n , β = κ n r L c air D n , δ = D c D n , γ = V b n r g ( ρ b − ρ ) L νρ D n , S c = ν D n . (37)8 y n Numerical solutionAnalytical solution (a) y c Numerical solutionAnalytical solution (b)Figure 1: A comparison of the numerical and analytical solutions [(a) n and (b) c ]. Similarly, we can also derive the dimensionless boundary conditions,[ ∇ n − α r ( c ) n ∇ c ] · ˆ n = , c = , u · ˆ n = , ∀ x ∈ ∂ Ω top , (38) ∇ n · ˆ n = ∇ c · ˆ n = , u = , ∀ x ∈ ∂ Ω bot , (39)the periodic boundary condition is still used for the left and right boundaries.Now we can simulate the dimensionless CF system with proposed LB model by setting Φ = ( n , c ) T , D = " − α r ( c ) n δ , S = (0 , − β r ( c ) n ) T , ν = S c and F = − S c γ n ˆ z . Based on the results in Ref. [6], under some suitable parameters, the solutions of Eqs. (36a)-(36d) converge to homogeneous-in- x steady-state solutions of the following time-independentsystem: ∇ · [ ∇ n − α r ( c ) n ∇ c ] = , ∇ · δ ∇ c − β r ( c ) n = . (40)The analytical solutions n s ( y ) and c s ( y ) can be explicitly derived if c > c ∗ (i.e., r ( c ) = n s ( y ) = δ A β α ( α Ay ) , c s ( y ) = − α ln cos( α Ay )cos( α A ) ! , (41)where A is a positive constant and satisfies β Z d n s ( y ) dy = δ A tan( α A ) . (42)We now numerically study the system (36) in the domain Ω = [ − , × [0 ,
1] with α = β = δ = γ = S c =
500 and the following initial conditions, n ( x , y ) = π , c ( x , y ) = , u ( x , y ) = . (43)According to Eq. (42), one can determine A = π/
20. In our simulations, we take ∆ x = ∆ y = . ∆ t = . × − such that the relaxation factor is in a proper range. In order to determine9 n , c y Present ( n )Lee and KimPresent ( c )Lee and Kim (a) x y (b)Figure 2: (a) Vertical profiles of the bacteria and oxygen concentrations n and c , (b) velocity vector. whether the result reaches a stable state, the following criterion is adopted, P i , j (cid:16) n N + i , j − n Ni , j (cid:17)P i , j n N + i , j < × − , P i , j (cid:16) c N + i , j − c Ni , j (cid:17)P i , j c N + i , j < × − , (44)where n Ni , j and c Ni , j denote the bacteria and oxygen concentrations at position ( i ∆ x , j ∆ y ) and time N ∆ t . From Fig. 1, one can find the numerical results are in good agreement with the analyticalsolutions, and the relative errors of bacteria and oxygen concentrations are less than 1 × − . In this part, we consider an example with α = β = δ = . γ =
418 and
S c = ∆ t = × − due to the larger value of S c . The cut-o ff value c ∗ = .
3, while the truncated function r ( c ) is regularized as a continuously form, r ( c ) = h + c − c ∗ p ( c − c ∗ ) + ǫ i , (45)where ǫ is a positive constant close to zero, and is set as ∆ x below. The initial conditions aregiven by n ( x , y ) = y > . − .
01 sin(( x − . π ) , . c ( x , y ) = , u ( x , y ) = . (46)This problem is used to validate the present LB model through a comparison with the resultsreported by Lee and Kim [18]. To this end, we present the quasi-homogeneous-in- x verticalprofiles of bacteria and oxygen concentrations at t = .
22 in Fig. 2(a). The results in this figureshows that the bacteria increase towards the bottom of domain, this is because the chemotacticconvection is cut-o ff for oxygen levels below c ∗ . In addition, the convection structure can becaptured by the velocity field in Fig. 2(b).Now we discuss the influence of parameter α to the system. Increasing α denotes the increaseof bacteria directional swimming relative to their di ff usion ( β and δ are fixed). For this purpose,we fixed the other parameters to be β =
10 and δ =
1, and consider di ff erent values of α ( α = , , , . n y Present ( α = 1)Lee and KimPresent ( α = 2)Lee and KimPresent ( α = 4)Lee and KimPresent ( α = 5 . (a) c y Present ( α = 1)Lee and KimPresent ( α = 2)Lee and KimPresent ( α = 4)Lee and KimPresent ( α = 5 . (b)Figure 3: Vertical profiles of the bacteria and oxygen concentrations [(a) n and (b) c ] at t = .
22 for α = , , , . n y Present ( β = 7 . β = 10)Lee and KimPresent ( β = 20)Lee and KimPresent ( β = 40)Lee and Kim (a) c y Present ( β = 7 . β = 10)Lee and KimPresent ( β = 20)Lee and KimPresent ( β = 40)Lee and Kim (b)Figure 4: Vertical profiles of the bacteria and oxygen concentrations [(a) n and (b) c ] at t = .
22 for β = . , , , t = .
22 are shown, one can find that with the increase of α , the bacteria concentration increasesnear the surface of the domain, bacteria leave the lower part of the domain faster, and consumelittle oxygen.Then we consider the e ff ect of parameter β . The increase of β indicates the increase ofoxygen consumption compared with oxygen di ff usion ( α and δ are fixed). As shown in Fig. 4where α = δ = β = . , ,
20 and 40 at t = .
22, with the increase of β , the oxygenconcentration decreases at the same height, and the up-swimming bacteria increase. We note thatthe present results of di ff erent cases are in agreement with the previous work [18].Finally, we focus on the influence of parameter δ . Actually, the increase of δ indicates thatoxygen di ff usion increases and oxygen adds to entire domain faster ( α and β are fixed). Wetake α = β = δ = , , ,
15, and present the vertical profiles of bacteria, oxygenconcentrations at t = .
22 in Fig. 5, from which we can see that with the increase of δ , oxygendi ff uses into the entire domain faster, while most of bacteria stay at the lower part of the domainand the up-swimming bacteria decrease. In this part, we investigate the formation of plume structures caused by the high bacteriaconcentration near the surface under some certain parameters. In the following simulations, wetake the bacteria-typical parameters as α = δ = S c = ∆ t = . × − .11 n y δ = 1 δ = 5 δ = 10 δ = 15 (a) c y δ = 1 δ = 5 δ = 10 δ = 15 (b)Figure 5: Vertical profiles of the bacteria and oxygen concentrations [(a) n and (b) c ] at t = .
22 for δ = , , , The truncated function is given by r ( c ) = c > . , c < . . (47)We first consider the case with β =
20 and γ = n ( x , y ) = . + . ξ, c ( x , y ) = , u ( x , y ) = , (48)where ξ is a random number uniformly distributed in the interval [0 , n and oxygen concentration c in time. From this figure, we canobserve the instability phenomena at t = .
16 due to the large amount of bacteria. Then the insta-bility amplifies the random irregularity of initial data, and develops into four plumes at t = . c in most of the bottom half of the domain is less than the chemotaxiscut-o ff c ∗ = . t = .
9, the left two plumes are close to each other, and merge intoa large plume at t = .
44. Meanwhile, the two plumes at the right side also approach to eachother, and combine together at t =
4, this structure no longer changes in time.Then we cut characteristic bacteria density n r by half which means β =
10 and γ = n ( x , y ) = y > . − .
01 sin(( x − . π ) , . c ( x , y ) = , u ( x , y ) = . (49)We carry out some simulations, and present the results at di ff erent times in Fig. 7. From Fig.7(a), one can observe that the instability occurs at the lower edge of the high concentration layerat t = .
2. When the time is increased to t = .
3, three plumes are formed due to the modulationsin initial conditions (49). It is interesting that those plumes bounce upwards slightly at t = . t = .
4, which can be seen from Fig. 7(c). Theplumes have some minor changes, and continue dropping to the bottom of the domain slightlyfrom t = . t = .
0. In addition, we also plot the oxygen concentration at t = . ff value c ∗ = . (a) t=0.16 -2 0 200.51 (b) t=0.16 -2 0 200.51 t=0.54 -2 0 200.51 t=0.54 . -2 0 200.51 t=0.9 -2 0 200.51 t=0.9 . -2 0 200.51 t=2.44 -2 0 200.51 t=2.44 . -2 0 200.51 t=4 -2 0 200.51 t=4 . Figure 6: Evolution of (a) bacteria concentration n and (b) oxygen concentration c in time. (a) t=0.2 -3 -2 -1 0 1 2 300.51 t=0.3 -3 -2 -1 0 1 2 300.51 t=0.4 -3 -2 -1 0 1 2 300.51 t=0.5 -3 -2 -1 0 1 2 300.51 t=1 -3 -2 -1 0 1 2 300.51 (b) t=0.5 t=0.4 . . . . . . . -3 -2 -1 0 1 2 300.51 (c) -3 -2 -1 0 1 2 300.51 t=0.5 . . . . . . . Figure 7: (a) Evolution of bacteria concentration n in time, (b) oxygen concentration c at t = .
5, (c) level sets of bacteriaconcentration n at t = . , . .2. Numerical results and discussion on the double-di ff usive convection system Similar to above discussion, we can derive the DDC system with Soret and Dufour e ff ectsfrom CDF system by setting Φ = ( T , C ) T , D = " α k CT k TC D , S = , F = g [1 − β T ( T − T ) − β C ( C − C )] and p = p /ρ , T t + ∇ · T u = ∇ · [ α ∇ T + k CT ∇ C ] , (50a) C t + ∇ · C u = ∇ · [ k TC ∇ T + D ∇ C ] , (50b) ∇ · u = , (50c) u t + u · ∇ u = −∇ p /ρ + ∇ · ν ∇ u + g [1 − β T ( T − T ) − β C ( C − C )] , (50d)where T is the temperature and C is the concentration. α and D are the thermal di ff usivity andmass di ff usivity. k CT and k TC are the Dufour and Soret coe ffi cients, respectively. In Eq. (50d),the Boussinesq approximation is adopted here to consider the density in the buoyancy term with T and C being the reference temperature and concentration, β T and β C being the coe ffi cients ofthermal and solute expansion.The boundary conditions of this system in 2D cavity [0 , L ] × [0 , H ] are given by ∇ T · ˆ n = , ∇ C · ˆ n = , u = , at y = y = H , (51a) T = T h , C = C h , u = , at x = , (51b) T = T l , C = C l , u = , at x = L , (51c)where T h and C h are the higher temperature and concentration, while T l and C l are the lowerones. Additionally, the initial conditions are T = T l , C = C l , u = , at t = . (52)The governing equations (50a)-(50d) can be expressed as a dimensionless form through in-troducing the following variables [42, 54]: x ′ = x L , t ′ = α L t , u ′ = L u α , p ′ = L α ρ p , T ′ = T − T T h − T l , C ′ = C − C C h − C l . (53)After dropping the prime notation in the rescaled variables, one can obtain the following dimen-sionless system: T t + ∇ · T u = ∇ · [ ∇ T + D CT ∇ C ] , (54a) C t + ∇ · C u = ∇ · [ S TC Le ∇ T + Le ∇ C ] , (54b) ∇ · u = , (54c) u t + u · ∇ u = −∇ p + ∇ · Pr ∇ u + RaPr ( T + NcC )ˆ z , (54d)where Prandtl number Pr , Rayleigh number Ra , buoyancy ratio Nc , Lewis number Le , Dufourfactor D CT , Soret factor S TC , and aspect ratio A are defined as Pr = να , Ra = g β T ( T h − T l ) L να , Nc = β C ( C h − C l ) β T ( T h − T l ) , Le = α D , D CT = k CT ( C h − C l ) α ( T h − T l ) , S TC = k TC ( T h − T l ) D ( C h − C l ) , A = HL . (55)15 able 1: A grid-independence study at Ra = , Pr = . Le = . A = . Nc = − . S TC = . D CT = . Grid number Nu Deviation (%)
S h
Deviation (%)250 × A × A × A Table 2: A comparison of Nu and S h between the present work and Refs. [42, 54] at Pr = . Le = . Nc = − . S TC = . D CT = . Present Ref. [42] Ref. [54]Rayleigh number
Ra Nu S h Nu S h Nu S h ∇ T · ˆ n , ∇ C · ˆ n = , u = , at y = y = A , (56a) T = , C = , u = , at x = , (56b) T = , C = , u = , at x = . (56c) T = , C = , u = , at t = . (57)The present LB model can also simulate above dimensionless system with Φ = ( T , C ) T , D = " D CT S TC / Le / Le , S = , ν = Pr and F = RaPr ( T + NcC )ˆ z .To characterize the heat and mass transfer in DDC system, the local Nusselt number Nu ( y )and Sherwood number S h ( y ), and the average Nusselt number Nu and Sherwood number S h onhigh temperature and concentration wall are defined as Nu ( y ) = ∂ T (0 , y ) ∂ x + D CT ∂ C (0 , y ) ∂ x , Nu = A Z A Nu ( y ) dy , (58a) S h ( y ) = S TC ∂ T (0 , y ) ∂ x + ∂ C (0 , y ) ∂ x , S h = A Z A S h ( y ) dy , (58b)where the space derivative terms in Nu ( y ) and S h ( y ) are calculated with the second-order upwinddi ff erence schemes. We first conduct a grid-independence study on DDC system in the cavity [0 , × [0 , A ] at Ra = , Pr = . Le = . A = . Nc = − . S TC = . D CT = .
1. The results of theaverage Nusselt number Nu and Sherwood number S h at di ff erent grid numbers of 250 × A ,16 able 3: The average Nusselt and Sherwood numbers of the DDC system at Ra = , Pr = . Le = . Nc = − . S TC = . D CT = . Present Ref. [42] Deviation (%)Aspect ratio
A Nu S h Nu S h Nu S h / / / × A , and 350 × A are presented in Table 1. From this table, one can find that therelative deviations of Nu and S h are less than 0.08%. Based on these results, one can find thatthe grid size 300 × A is fine enough, and is also used in the following simulations. ff ects of some physical parameters We first consider the e ff ect of Rayleigh number Ra . To this end, some other parameters arefixed as Pr = . Le = . Nc = − . A = . S TC = . D CT = .
1. As we can seefrom Table 2, the average Nusselt and Sherwood numbers increase with the increase of Ra . At Le = .
0, the concentration di ff usion is weaker than the thermal di ff usion, which gives a largerconcentration gradient at boundary layer with a high concentration, thus the average Sherwoodnumber is larger than the average Nusselt number. We note that present results are in goodagreement with those in the previous works [42, 54], as shown in Table 2 and Fig. 8, where therelative errors are no more than 2%.Then the e ff ect of aspect ratio A is also investigated. We perform some simulations at Ra = , Pr = . Le = . Nc = − . S TC = . D CT = .
1, and show the results in Table 3where the aspect ratio A is changed from 1 / A is increased up to 1 .
5, and decreases when the aspect ratio isfurther increased, this means that there is a critical aspect ratio A = .
5. Moreover, throughchanging the Rayleigh number Ra , we can obtain scaling relations between average Nusselt andSherwood numbers and Rayleigh number through fitting the data, as shown in Fig. 9. As seenfrom this figure, when A < .
0, the increases of Nu and S h are faster than the cases with A > . A = . Ra = are time-averaged quantities.Next the influence of buoyancy ratio Nc is also studied through fixing other parameters as Ra = , Pr = . Le = . S TC = . D CT = .
1. One can see from Fig. 10 that the Nu and S h at all aspect ratios are close to their minimum values at Nc = − .
0, where the thermal andsolutal buoyancy forces are o ff set, and the results are more accurate for the horizontal cavities.Here it is also worth noting that the values of average Nusselt and Sherwood numbers are time-averaged at Nc = − . A = . ff ects of the Soret and Dufour factors at Ra = , Pr = . Le = . Nc = − . A = .
0. When Soret and Dufour factors are equal to each other, Nu and S h increase with the increase of Soret and Dufour factors, as presented in Table 4 where the17 emperature . . . . . . . . . . x y Concentration . . . . . . . . . . x y x y Streamlines (a)(b)Figure 8: Temperature, concentration and streamlines at Ra = , Pr = . Le = . Nc = − . S TC = . D CT = .
1. [(a) present results, (b) results in Ref. [42]]. Ra N u Nu = 0 . Ra . − − − → Nu = 0 . Ra . − − − → Nu = 0 . Ra . − − − → A = 2 . A = 1 . A = 0 . (a) Ra S h Sh = 0 . Ra . − − − → Sh = 0 . Ra . − − − → Sh = 0 . Ra . − − − → A = 2 . A = 1 . A = 0 . (b)Figure 9: The average Nusselt number (a) and average Sherwood number (b) at di ff erent aspect ratios. -4 -3 -2 -1 0 1 2 Nc N u A = 2 . A = 1 . A = 0 . (a) -4 -3 -2 -1 0 1 2 Nc S h A = 2 . A = 1 . A = 0 . (b)Figure 10: The average Nusselt number (a) and average Sherwood number (b) at di ff erent aspect ratios. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 S TC N u / S h NuSh (a) -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 D CT N u / S h NuSh (b)Figure 11: The average Nusselt and Sherwood numbers at di ff erent Soret and Dufour factors, [(a) Soret factor S TC at D CT = D CT at S TC = able 4: The average Nusselt and Sherwood numbers of the DDC system at Ra = , Pr = . Le = . Nc = − . A = . Present Ref. [42]Soret factor S TC Dufour factor D CT Nu S h Nu S h -0.2 -0.2 2.5609 4.3020 - --0.1 -0.1 2.7079 4.4842 - -0 0 2.8265 4.6144 2.8275 4.62320.05 0.05 2.8747 4.6597 2.8753 4.66810.1 0.1 2.9151 4.6914 2.9153 4.69930.15 0.15 2.9475 4.7091 2.9473 4.71640.2 0.2 2.9715 4.7123 2.9709 4.7191relative errors are within 0.2%, compared with [42]. In addition, we take D CT =
0, and presentthe e ff ect of Soret factor in Fig. 11(a). From this figure one can see that the average Sherwoodnumber increases in Soret factor, while the average Nusselt number first decreases when S TC is less than − .
5, and then increases in Soret factor when S TC is larger than − .
5. However,as shown in Fig. 11(b), when we fix S TC = D CT but decreases in a small region near D CT =
4. Conclusions
In this paper, a LB model is proposed for the general CDF system in d dimensional space.Through the direct Taylor expansion, the governing equations of the general system can be recov-ered correctly from the proposed LB model. In the present LB model, the cross-di ff usion termsare modeled through introducing some extra collision operators, and the proper auxiliary sourceterms are also constructed, which can be used to avoid the discretization of some derivative terms.Additionally, the computational scheme of pressure actually contains the external force but canbe simplified under the incompressible condition. We validate the present LB model by two im-portant CDF systems, i.e., CF system and DDC system with Soret and Dufour e ff ects, and findthe results are in good agreement with the analytical solutions and the previous works. Moreover,in the DDC system, the influence of Rayleigh number at di ff erent aspect ratios is considered, andthe power-law relations between the average Nusselt and Sherwood numbers and the Rayleighnumber are observed, which indicate that the heat and mass transfer rates increase faster withthe increase of the Rayleigh number in horizontal cavities. In addition, the average Nusselt andSherwood numbers are close to their minimum values when buoyancy ratio approaches to − ff ects of Soret and Dufour factorsare also investigated, and the results show that the average Nusselt number increases in a largeinterval of Soret and Dufour factors, while the average Sherwood number increases with Soretfactor but decreases with the Dufour factor. Acknowledgments
This work is supported by the National Natural Science Foundation of China (Grant No.51836003), and the National Key Research and Development Program of China (Grant No.20017YFE0100100).
Appendix A. The lattice Boltzmann model for incompressible Navier-Stokes equations
In this appendix, we will adopt the direct Taylor expansion method to recover the incompress-ible Navier-Stokes equations from the present LB model. We first apply the Taylor expansion toEq. (2b), similar to Eq. (13), and can obtain the following equations at di ff erent orders of ∆ t , D i h eqi = − ω ∆ t h nei + (1 − ω F i + O ( ∆ t ) , (A.1a) D i ( h eqi + h nei ) + ∆ t D i h eqi = − ω ∆ t h nei + (1 − ω F i + O ( ∆ t ) , (A.1b)where h nei is the non-equilibrium part of h i . According to Eq. (A.1a), we have ∆ t D i h eqi = − D i ω h nei + ∆ t D i (1 − ω F i + O ( ∆ t ) . (A.2)Then substituting Eq. (A.2) into Eq. (A.1b), one can obtain D i h eqi + D i (1 − ω h nei + ∆ t D i (1 − ω F i = − ω ∆ t h nei + (1 − ω F i + O ( ∆ t ) . (A.3)To give the correct Navier-Stokes equations, the distribution functions h i , h eqi and F i satisfy thefollowing conditions, X i h i = X i h eqi = ρ , X i c i h eqi = u , X i c i c i h eqi = p I + uu , X i c i c i c i h eqi = ∆ · u , (A.4a) X i F i = , X i c i F i = F , X i c i c i F i = ϕ ( uF + Fu ) , (A.4b)where ∆ is a fourth-order tensor giving by δ αβ δ θγ + δ βθ δ αγ + δ αθ δ βγ . From Eqs. (A.4a) and (5b),we can get X i c i h nei = X i c i h i − X i c i h eqi = − ∆ t F . (A.5)Using above relations, we can derive the zeroth and first order moments of Eqs. (A.1a) and (A.3)at the orders of O ( ∆ t ) and O ( ∆ t ), ∇ · u = O ( ∆ t ) , (A.6a) ∂ t u + ∇ · ( p I + uu ) = F + O ( ∆ t ) , (A.6b) ∇ · u = O ( ∆ t ) , (A.7a) ∂ t u + ∇ · ( p I + uu ) + ∇ · (1 − ω X i c i c i h nei + ∆ t ϕ ( uF + Fu ) = F + O ( ∆ t ) , (A.7b)21here the term P i c i c i h nei can be evaluated by Eq. (A.1a), X i c i c i h nei = − ∆ t ω X i c i c i (cid:20) D i h eqi − (1 − ω F i (cid:21) + O ( ∆ t ) = − ∆ t ω ∂ t X i c i c i h eqi + ∇ · X i c i c i c i h eqi − (1 − ω X i c i c i F i + O ( ∆ t ) = − ∆ t ω (cid:20) ∂ t ( p I + uu ) + ∇ · ( ∆ · u ) − (1 − ω ϕ ( uF + Fu ) (cid:21) + O ( ∆ t ) = − ∆ t ω (cid:20) ∂ t ( p I + uu ) + c s [ ∇ u + ( ∇ u ) T ] − (1 − ω ϕ ( uF + Fu ) (cid:21) + O ( ∆ t ) . (A.8)Based on Eq. (A.6), we get ∂ t ( u α u β ) = u α ( F β − ∇ β p − ∇ γ u β u γ ) + ( F α − ∇ α p − ∇ γ u α u γ ) u β + O ( Ma ∆ t ) = u α F β + F α u β − u α ∇ β p − ∇ α pu β − ∇ γ u α u β u γ + O ( Ma ∆ t ) = u α F β + F α u β + O ( Ma ∆ t + Ma ) , (A.9)where Ma is the Mach number. Substituting Eq. (A.9) into Eq. (A.8) and with the help of ∂ t p = O ( Ma ), one can obtain X i c i c i h nei + ∆ t ϕ ( uF + Fu ) = ∆ t ϕ − ω ( uF + Fu ) − c s ∆ t ω [ ∇ u + ( ∇ u ) T ] + O ( Ma ∆ t + ∆ t ) . (A.10)According to the value of ϕ , we have the following two cases. Case 1: ϕ =
0. Due to the fact uF = O ( Ma ), the first term on the right hand side of Eq. (A.10)can be absorbed into the truncation term O ( Ma ∆ t + ∆ t ). Actually, the derivation process ofEq. (A.9) is not necessary since ∂ t uu = O ( Ma ). In this case, the force term F i can be simplified. Case 2: ϕ =
1. Under this condition, the first term on the right hand side of Eq. (A.10) is zero.In both cases, if we substitute Eq. (A.10) into Eq. (A.7b), and omit the truncation terms O ( ∆ t ) in Eq. (A.7a) and O ( Ma ∆ t + ∆ t ) in Eq. (A.7b), one can get the incompressible Navier-Stokes equations: ∇ · u = , (A.11a) u t + ∇ · uu = −∇ p + ∇ · ν [ ∇ u + ( ∇ u ) T ] + F , (A.11b)where ν = ( ω − ) c s ∆ t .Now let focus on the computation of pressure. From Eq. (A.1a) one can obtain h nei = − ∆ t ω (cid:20) D i h eqi − (1 − ω F i (cid:21) + O ( ∆ t ) . (A.12)Taking the zeroth-direction of Eqs. (3b) and (A.12), we have h eq = ( W − pc s + ρ − W u · u c s , (A.13) h ne = − ∆ t ω (cid:20) ∂ t h eq − (1 − ω F (cid:21) + O ( ∆ t ) . (A.14)22wing to the fact that ∂ t h eq is order of O ( Ma ), Eq. (A.14) can be simplified by h ne = ∆ t ω (1 − ω F + O ( Ma ∆ t + ∆ t ) . (A.15)According to Eqs. (A.13) and (A.15), and based on h i = h eqi + h nei , we have(1 − W ) pc s = ρ − ( h − h ne ) − W u · u c s = ρ − X i h i − X i , h i − W u · u c s + ∆ t ω (1 − ω F + O ( Ma ∆ t + ∆ t ) = X i , h i − W u · u c s + ∆ t ω (1 − ω F + O ( Ma ∆ t + ∆ t ) . (A.16)Ignoring the truncation error terms of O ( Ma ∆ t + ∆ t ), we can obtain the computational schemefor pressure, p = c s − W X i , h i − W u · u c s + ∆ t ω (1 − ω F = c s − W X i , h i − W u · u c s − ϕ ∆ t ( 1 ω −
12 ) u · F c s . (A.17) References [1] E. F. Keller, L. A. Segel, Initiation of slime mold aggregation viewed as an instability, Journal of TheoreticalBiology 26 (1970) 399–415.[2] E. F. Keller, L. A. Segel, Model for chemotaxis, Journal of Theoretical Biology 30 (1971) 225–234.[3] N. Shigesada, K. Kawasaki, E. Teramoto, Spatial segregation of interacting species, Journal of Theoretical Biology79 (1) (1979) 83–99.[4] S. Hittmeir, H. Ranetbauer, C. Schmeiser, M.-T. Wolfram, Derivation and analysis of continuum models for cross-ing pedestrian tra ffi c, Mathematical Models and Methods in Applied Sciences 27 (07) (2017) 1301–1325.[5] C. F. Curtiss, R. B. Bird, Multicomponent di ff usion, Industrial & Engineering Chemistry Research 38 (7) (1999)2515–2522.[6] A. J. Hillesdon, T. J. Pedley, J. O. Kessler, Bioconvection in suspensions of oxytactic beacteria: linear theory,Bulletin of Mathematical Biology 57 (1995) 299–303.[7] C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Goldstein, J. O. Kessler, Self-concentration and large-scale co-herence in bacterial dynamics, Physical Review Letters 93 (2004) 098103.[8] I. Tuval, L. Cisneros, C. Dombrowski, C. W. Wolgemuth, J. O. Kessler, R. E. Goldstein, Bacterial swimmingand oxygen transport near contact lines, Proceeding of the National Academy of Sciences of the United States ofAmerica 102 (7) (2005) 2277–2282.[9] O. V. Trevisan, A. Bejan, Combined heat and mass transfer by natural convection in a vertical enclosure, Journalof Heat Transfer 109 (1) (1987) 104–112.[10] S. N. Gaikwad, M. S. Malashetty, K. R. Prasad, An analytical study of linear and non-linear double di ff usiveconvection with Soret and Dufour e ff ects in couple stress fluid, International Journal of Non-Linear Mechanics42 (7) (2007) 903–913.[11] M. A. Budroni, Cross-di ff usion-driven hydrodynamic instabilities in a double-layer system: General classificationand nonlinear simulations, Physical Review E 92 (6) (2015) 063007.[12] M. Bendahmane, F. Karami, M. Zagour, Kinetic-fluid derivation and mathematical analysis of the cross-di ff usion–Brinkman system, Mathematical Methods in the Applied Sciences 41 (16) 6288–6311.[13] M. C. Kim, K. H. Song, Theoretical and numerical analyses of the e ff ect of cross-di ff usion on the gravitationalinstability in ternary mixtures, International Journal of Heat and Mass Transfer 43 (2019) 118511.
14] K. R. Raghunatha, I. S. Shivakumara, M. S. Swamy, E ff ect of cross-di ff usion on the stability of a triple-di ff usiveOldroyd-B fluid layer, Zeitschrift f¨ur angewandte Mathematik und Physik 70 (100) (2019).[15] A. Atlas, M. Bendahmane, F. Karami, D. Meskine, Kinetic-fluid derivation and mathematical analysis of anonlocalcross-di ff usion–fluid system, Applied Mathematical Modelling 82 (2020) 379–408.[16] A. Chertock, K. Fellner, A. Kurganov, A. Lorz, P. A. Markowich, Sinking, merging and stationary plumes in acoupled chemotaxis-fluid model: a high-resolution numerical approach, Journal of Fluid Mechanics 694 (2012)155–190.[17] T. W. Sheu, C. Y. Chiang, Numerical investigation of chemotaxic phenomenon in incompressible viscous fluid flow,Computer and Fluids 103 (2014) 290–306.[18] H. G. Lee, J. Kim, Numerical investigation of falling bacterial plumes caused by bioconvection in a three-dimensional chamber, European Journal of Mechanics B / Fluids 52 (2015) 120–130.[19] Y. Deleuze, C.-Y. Chiang, M. Thiriet, T. W. Sheu, Numerical study of plume patterns in a chemotaxis-di ff usion-convection coupling system, Computer and Fluids 126 (2016) 58–70.[20] N. Nithyadevi, R.-J. Yang, Double di ff usive natural convection in a partially heated enclosure with Soret and Dufoure ff ects, International Journal of Heat and Fluid Flow 30 (2009) 902–910.[21] O. A. B´eg, V. R. Prasad, B. Vasu, N. B. Reddy, Q. Li, R. Bhargava, Free convection heat and mass transfer froman isothermal sphere to a micropolar regime with Soret / Dufour e ff ects, International Journal of Heat and MassTransfer 54 (2011) 9–18.[22] C.-Y. Cheng, Soret and Dufour e ff ects on free convection heat and mass transfer from an arbitrarily inclined platein a porous medium with constant wall temperature and concentration, International Communications in Heat andMass Transfer 39 (2012) 72–77.[23] J. Wang, M. Yang, Y. Zhang, Onset of double-di ff usive convection in horizontal cavity with Soret and Dufoure ff ects, International Journal of Heat and Mass Transfer 78 (2014) 1023–1031.[24] J. Wang, M. Yang, Y.-L. He, Y. Zhang, Oscillatory double-di ff usive convection in a horizontal cavity with Soretand Dufour e ff ects, International Journal of Thermal Sciences 106 (2016) 57–69.[25] F. J. Higuera, S. Succi, R. Benzi, Lattice gas dynamics with enhanced collisions, Europhysics Letters 9 (4) (1989)345–349.[26] R. Benzi, S. Succi, M. Vergassloa, The lattice Boltzmann equation: theory and applications, Physics Reports 222(1992) 145–197.[27] Y. Qian, S. Succi, S. A. Orszag, Recent advances in lattice Boltzmann computing, Annual Reviews of Computa-tional Physics 3 (1995) 195–242.[28] S. Chen, G. D. Doolen, Lattice Boltzmann method for fluid flows, Annual Reviews of Fluid Mechanics 30 (1998)329–364.[29] C. K. Aidun, J. R. Clausen, Lattice-Boltzmann method for complex flows, Annual Review of Fluid Mechanics42 (1) (2010) 439–472.[30] H. Wang, X. Yuan, H. Liang, Z. Chai, B. Shi, A brief review of the phase-field-based lattice Boltzmann method formultiphase flows, Capillarity 2 (2) (2019) 33–52.[31] S. P. Dawson, S. Chen, G. D. Doolen, Lattice Boltzmann computations for reaction-di ff usion equations, The Journalof Chemical Physics 98 (2) (1993) 1514.[32] R. Blaak, P. M. Sloot, Lattice dependence of reaction-di ff usion in lattice Boltzmann modeling, Computer PhysicsCommunications 129 (2000) 256–266.[33] X. He, N. Li, B. Goldstein, Lattice Boltzmann simulation of di ff usion-convection systems with surface chemicalreaction, Molecular Simulation 25 (3) (2000) 145–156.[34] B. Shi, Z. Guo, Lattice Boltzmann model for nonlinear convection-di ff usion equation, Physical Review E 79 (2009)016701.[35] Z. Chai, T. S. Zhao, Lattice Boltzmann model for the convection-di ff usion equation, Physical Review E 87 (2013)063309.[36] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, 2001.[37] T. Kr¨uger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. M. Viggen, The Lattice Boltzmann Method:Principles and Practice, Oxford University Press, 2017.[38] M. Hilpert, Lattice Boltzmann model for bacteria chemotaxis, Journal of Mathematical Biology 51 (2005) 302–332.[39] X. Yu, Z. Guo, B. Shi, Numerical study of cross di ff usion e ff ects on double di ff usion convection with latticeBoltzmann method, International Conference on Computational Science 4487 (2007) 810–817.[40] X. Yang, B. Shi, Z. Chai, Coupled lattice Boltzmann method for generalized Keller-Segel chemotaxis model,Computers and Mathematics with Applications 68 (2014) 1653–1670.[41] Z. Chai, H. Liang, R. Du, B. Shi, A lattice Boltzmann model for two-phase flow in porous media, SIAM Journalon Scientific Computing 41 (4) (2019) B746–B772.[42] Q. Ren, C. L. Chan, Numerical study of double-di ff usive convection in a vertical cavity with soret and dufoure ff ects by lattice Boltzmann method on GPU, International Journal of Heat and Mass Transfer 93 (2016) 538–553.
43] C. Huber, B. Chopard, M. Manga, A lattice Boltzmann model for coupled di ff usion, Journal of ComputationalPhysics 229 (20) (2010) 7956–7976.[44] Z. Chai, X. Guo, L. Wang, B. Shi, Maxwell-stefan-theory-based lattice Boltzmann model for di ff usion in multi-component mixtures, Physical Review E 99 (2019) 023312.[45] N. He, N. Wang, B. Shi, Z. Guo, A unified incompressible lattice BGK model and its application to three-dimensional lid-driven cavity flow, Chinese Physics 13 (01) (2004) 40–46.[46] D. J. Holdych, D. R. Noble, J. G. Georgiadis, R. O. Buckius, Truncation error analysis of lattice Boltzmann meth-ods, Journal of Computational Physics 193 (2) (2004) 595–619.[47] A. J. Wagner, Thermodynamic consistency of liquid-gas lattice Boltzmann simulations, Physical Review E 74(2006) 056703.[48] G. Kaehler, A. J. Wagner, Derivation of hydrodynamics for multi-relaxation time lattice Boltzmann using themoment approach, Communications in Computational Physics 13 (3) (2013) 614–628.[49] Z. Chai, B. Shi, Multiple-relaxation-time lattice Boltzmann method for the navier-stokes and nonlinear convection-di ff usion equations: Modeling, analysis and elements, Physical Review E 102 (2020) 023306.[50] Z. Chai, B. Shi, Z. Guo, A multiple-relaxation-time lattice Boltzmann model for general nonlinear anisotropicconvection–di ff usion equations, Journal of Scientific Computing 69 (2016) 355–390.[51] T. Zhang, B. Shi, Z. Guo, Z. Chai, J. Lu, General bounce-back scheme for concentration boundary condition in thelattice Boltzmann method, Physical Review E 85 (2012) 016701.[52] A. J. C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1.Theoretical foundation, Journal of Fluid Mechanics 271 (1994) 285–309.[53] A. J. C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 2.Numerical results, Journal of Fluid Mechanics 271 (1994) 311–339.[54] H. Yu, Z. Luo, Q. Lou, S. Zhang, J. Wang, Lattice Boltzmann simulations of the double-di ff usive natural convectionand oscillation characteristics in an enclosure with Soret and Dufour e ff ects, International Journal of ThermalSciences 136 (2019) 159–171.ects, International Journal of ThermalSciences 136 (2019) 159–171.