A lattice Boltzmann model for the fractional advection–diffusion equation coupled with incompressible Navier–Stokes equation

Applied Mathematics Letters 101 (2020) 106074

Contents lists available at ScienceDirect

Applied Mathematics Letters www.elsevier.com/locate/aml

A lattice Boltzmann model for the fractional advection–diffusion equation coupled with incompressible Navier–Stokes equation Rui Du ∗, Zixuan Liu School of Mathematics, Southeast University, Nanjing 210096, PR China

article

info

Article history: Received 10 July 2019 Received in revised form 25 September 2019 Accepted 26 September 2019 Available online 4 October 2019 Keywords: Fractional advection–diffusion equation Caputo derivative Lattice Boltzmann method

abstract The fractional advection–diffusion problem coupled with incompressible Navier– Stokes equations is important in science and engineering. In this paper, a fresh lattice Boltzmann model with double-distribution functions for this problem is presented. One distribution function is used for pressure and velocity. The other is applied for concentration. Through Chapman–Enskog expansion, the macroscopic continuum equations can be recovered from the LB model. Numerical simulations are carried out to verify the numerical accuracy and efficiency. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Fractional advection–diffusion equations (FADEs) [1] have attracted more and more attention because of their frequent appearance in various applications in science and engineering. They are more appropriate for the description of memorial and hereditary properties of various materials and processes than the classical integer-order differential equations. However, most FADEs cannot get the exact solution. So various kinds of numerical methods have been used to seek the approximation solution for the FADEs, such as finite element method [2], Parareal algorithms [3], finite difference method [4–7], and so on. Compared to the traditional methods, the LB method has some distinct advantages, including the simplicity of program, locality of computation and ease in dealing with complex boundary. Nowadays, the LB method has also gained great successes in a variety of fields, for example, fluid flows in porous media, turbulence, multiphase and multicomponent flows, to name but a few [8]. On the other hand, as a numerical solver to nonlinear PDEs, the LB method has also been adopted to solve some special PDEs, such as Korteweg–de Vries equation [9], Cahn–Hilliard and Allen–Cahn equations [10,11], Sine–Cosine–Gordon equation [12], FADEs [13–16] and so on. But there is less research on advection–convection coupled with fluid field. ∗ Corresponding author. E-mail addresses: [email protected] (R. Du), [email protected] (Z. Liu).

https://doi.org/10.1016/j.aml.2019.106074 0893-9659/© 2019 Elsevier Ltd. All rights reserved.

2

R. Du and Z. Liu / Applied Mathematics Letters 101 (2020) 106074

In this paper, the incompressible Navier–Stoles (N–S) equations and time fractional advection–diffusion equation in Caputo sense are considered ∇ · u = 0, ∂u + ∇ · (uu) = −∇p + ν∇2 u, ∂t C α 0 Dt ϕ + ∇ · (uϕ) = D∆ϕ + g(x, t),

(1.1) (1.2) (1.3)

where x ∈ Ω , 0 < t ⩽ T . Ω is the physical domain. u is the velocity of flow and ϕ is the concentration. g(x, t) α is the source term. D is diffusion coefficient and ν is viscosity. The fractional derivative term C 0 Dt ϕ(x, t) is given by ∫ t ⎧ ∂ϕ 1 ⎨ (t − ξ)−α (x, ξ)dξ, 0 ⩽ α < 1, C α Γ (1 − α) 0 ∂ξ 0 Dt ϕ(x, t) = ⎩ ϕt (x, t), α = 1, where Γ is the gamma function. The rest of this paper is arranged as follows. In Section 2, a novel lattice Boltzmann model with double distribution functions for fractional convection–diffusion equation coupled with incompressible N–S equations is proposed and Eqs. (1.1)–(1.3) can be recovered from the LB model through Chapman–Enskog analysis. Detailed numerical simulation is performed to demonstrate the efficiency of the LB model in Section 3. The paper ends with a brief conclusion in Section 4. 2. Lattice Boltzmann model for FADE coupled with incompressible N–S equations In this study, the D2Q9 lattice model is used. The velocity and the corresponding weights are expressed as =

ci

ωi

=

[ 0 1 c 0 0 ⎧ 4 ⎪ ⎪ , ⎪ ⎪ ⎪ ⎨9 1 , ⎪ 9 ⎪ ⎪ ⎪ ⎪ ⎩ 1, 36

0 1

−1 0

1 −1 1 1

0 −1

−1 −1

] 1 , −1

i=0

(2.1)

i = 1, . . . , 4 ,

c2s

2

= c /3.

i = 5, . . . , 8

where c = ∆x ∆t . ∆x and ∆t are the lattice spacing and time step respectively. ci is the lattice velocity. cs is the so-called lattice sound speed. ωi is the weight coefficient. The weight coefficient ωi and discrete ci in D2Q9 model satisfy the following conditions: ∑ ∑ ∑ ωi = 1, cip = 0, ωi cip ciq = δpq c2s , i ∑

i

ωi cip ciq cir = 0,

i

i

∑

cip ciq cir cil = ∆pqrl c2s .

i

where δpq and ∆pqrl = δpq δrl + δpr δql + δpl δqr are both the Kronecker tensors. α At t = tn , the time fractional derivative C 0 Dt ϕ(x, tn ) can be approximated by [14] ( n−1 ∫ ∫ tn ∂ϕ (x, ξ) ) ∑ tm ∂ϕ 1 ∂ξ (x, ξ) ∂ξ C α dξ + dξ 0 Dt ϕ(x, tn ) = Γ (1 − α) m=1 tm −1 (t − ξ)α (t − ξ)α tn −1 ≈

n−1 ∆t−α ∑ [(n − m + 1)1−α − (n − m)1−α ] · [ϕ(x, tm ) − ϕ(x, tm−1 )] Γ (2 − α) m=1

+

∆t2−α ∂ϕ (x, tn ). Γ (2 − α) ∂t

(2.2)

R. Du and Z. Liu / Applied Mathematics Letters 101 (2020) 106074

Let σ =

∆t1−α Γ (2−α) ,

3

and substitute Eq. (2.2) into Eq. (1.3). Then we obtain σ

∂ϕ (x, tn ) + ∇ · (uϕ) = D∇2 ϕ + F (x, tn ), ∂t

(2.3)

where F (x, tn ) = −

n−1 σ ∑ [(n − m + 1)1−α − (n − m)1−α ] · [ϕ(x, tm ) − ϕ(x, tm−1 )] + g(x, tn ). ∆t m=1

(2.4)

Based on the idea by Du et al. [14] and Guo et al. [17], the evolution equations can be written as 1 [hi (x, t) − heq i (x, t)] , τh ) ( 1 1 F (x, t), fi (x + ci ∆t, t + ∆t) − fi (x, t) = − [fi (x, t) − fieq (x, t)] + ∆t ωi 1 − τf 2τf

hi (x + ci ∆t, t + ∆t) − hi (x, t) = −

where hi (x, t) and fi (x, t) are the distribution functions for fluid field and concentration field respectively. τh and τf are the dimensionless relaxation parameters. The equilibrium distribution function heq i is defined by ⎧ p ⎪ i = 0, ⎨ (ω0 − 1) 2 + s0 (u), c s heq (x, t) = i p ⎪ ⎩ ωi 2 + si (u), i = 1 ∼ 8, cs 2

2

i ·u) where si (u) = ωi ( cci 2·u + (c2c − |u| ). 4 2c2 s s s The equilibrium distribution function fieq (x, t) can be given by [ ] M : (ci ci − c2s I) ci · u fieq = ωi ϕ σ + + , 2c2s c2s

M = (β − σ)I,

where I is the unit matrix. β is a new parameter to adjust the relaxation parameter τf for a fixed diffusion coefficient D in order to keep LB model stable, and we should adjust the value of β to make sure fieq (x, t) non-negative. In addition, the equilibrium function satisfies the following conditions which can be verified ∑ eq ∑ ∑ hi = 0, ci heq ci ci heq (2.5) i = u, i = uu + pI, i

∑ i

i

fieq

= σϕ,

i

∑

ci fieq

∑

= uϕ,

i

ci ci fieq = βc2s ϕI.

(2.6)

i

Finally, the macro expression of the velocity u, pressure P and concentration ϕ is given by ( ) ∑ c2s ∑ 1 ∑ 1 u= ci hi , p = [ hi + s0 (u)], ϕ = fi + ∆tF . 1 − ω0 σ 2 i i

(2.7)

i̸=0

Through the Chapman–Enskog analysis, the impressible Navier–Stokes equations [i.e., Eqs. (1.1) and (1.2)] can be recovered with [17] ν = c2s (τh − 0.5)∆t. In the following, we will verify the validity of this model by Chapman–Enskog analysis and get the relationship between D and τf . The distribution function fi (x, t), the derivatives of time and space and the source term F are expanded as (0)

fi = fi

(1)

+ ϵfi

∂ ∂ ∂ =ϵ + ϵ2 , ∂t ∂t1 ∂t2

(2)

+ ϵ2 fi

∇ = ϵ∇1 ,

+ ··· , F = ϵF (1) .

(2.8) (2.9)

R. Du and Z. Liu / Applied Mathematics Letters 101 (2020) 106074

4

By Taylor expansion Eq. (2.5), one can obtain ∆tDi fi + where Di =

∂ ∂t

( ) 1 1 ∆t2 2 Di fi = − (fi − fieq ) + ∆tωi 1 − F, 2 τf 2τf

(2.10)

+ ci · ∇. Substituting Eqs. (2.8) and (2.9) into Eq. (2.10), we have [ ( )2 ] ( ) ∆t (0) (1) (2) 2 ∂ 2 ∂ + ϵDi1 + ϵ fi + ϵfi + ϵ2 fi + · · · ϵDi1 + ϵ ∂t2 2 ∂t2 ) ( ( ) 1 1 (0) (1) (2) fi − fieq + ϵfi + ϵ2 fi + · · · + ϵωi 1 − F (1) , =− τf ∆t 2τf

(2.11)

where Di1 = ∂t∂1 + ci · ∇1 . We will separate the terms in order of ϵ0 , ϵ1 and ϵ2 respectively, and then get the zeroth, first and second-order equations in ϵ, (0)

ϵ0 :

fi

= fieq , (0)

(2.12)

( ) 1 1 (1) f + ωi 1 − =− F (1) , τf ∆t i 2τf

ϵ1 :

Di1 fi

ϵ2 :

1 ∆t 2 (0) ∂fi (2) (1) D f =− f . + Di1 fi + ∂t2 2 i1 i τf ∆t i

(2.13)

(0)

(2.14)

Multiplying Di1 on both sides of Eq. (2.13) and substituting into Eq. (2.14), we can rewrite the second-order equation in ϵ as ] )[ ( (0) ∂fi 1 ∆t 1 (2) (1) 2 (1) ϵ : ωi F =− f . (2.15) + Di1 1 − fi + ∂t2 2τf 2 τf ∆t i (1)

Besides, we can derive the following of non-equilibrium distribution functionfi and (2.14) ∑ (1) ∆t ϵ1 : fi = − F (1) , 2 i ∑ (k) ϵ2 : fi = 0, (k ≥ 2).

from Eqs. (2.6), (2.13)

(2.16) (2.17)

i

From Eq. (2.13), we can obtain ∑

(1)

fi

i

( ) ] ∑[ 1 (0) Di1 fi − 1 − ωi F (1) 2τ [i ( ) ] ∂ 1 (1) = −τf ∆t (σϕ) + ∇1 · (uϕ) − 1 − F . ∂t1 2τf = −τf ∆t

(2.18)

With Eq. (2.16), we can derive ∂ σϕ + ∇1 · (uϕ) = F (1) . ∂t1 Then, we will use Eqs. (2.6) and (2.13) to get ) ] ( ∑ ∑ [ 1 (1) (0) (1) ci fi = −τf ∆t ci Di1 fi − ωi 1 − F 2τf i i ∂ = −τf ∆t (uϕ) − τf ∆tβc2s ∇1 ϕ. ∂t1 After a summation of Eq. (2.15), we can obtain ( ) [ ] ∂ 1 ∂ τf ∆t ∇1 (σϕ) − 1 − (uϕ) + βc2s ∇21 ϕ = 0. ∂t2 2τf ∂t1

(2.19)

(2.20)

(2.21)

R. Du and Z. Liu / Applied Mathematics Letters 101 (2020) 106074

5

It is noted that we neglect τf ∆t∇1 ∂t∂1 (uϕ) in order of ∆t. The recovered equation at second-order of ϵ is as follows ( ) ∂ 1 σϕ − 1 − τf ∆tβc2s ∇21 ϕ = 0. (2.22) ∂t2 2τf Multiplying ϵ and ϵ2 on both sides of Eqs. (2.19) and (2.22), we have recover Eq. (1.3) with ( ) 1 D = βc2s τf − ∆t. 2

(2.23)

3. Numerical simulation In this section, we give some numerical results to demonstrate the accuracy of the LB model proposed in this paper. In the simulation, non-equilibrium extrapolation method [18] is used to deal with all the boundary. Take the region Ω = {(x, y)|0 ⩽ x ⩽ 2π, 0 ⩽ y ⩽ 0.16π}, and consider Eqs. (1.1)–(1.3) with the source term ] [ 625 2 25 2 2−α t + Dt sin(x) sin( y) + ∇ · (uϕ). g(x, y, t) = Γ (3 − α) 4 2 This problem has the analytic solutions, u(x, y, t) = A, v(x, y, t) = B cos(kx − kAt)e−k p(x, y, t) = 3 + 0.01 sin(20πt), 25 ϕ(x, y, t) = t2 sin(x) sin( y), 2

2

νt

, (3.1)

where the constants A, B, k, ν are chosen as A = B = 0.1, k = ν = 1.0. The initial and boundary conditions are implemented according to the analytic solution. In the simulation, D = 1 × 10−4 , β = 0.0001. Table 1 shows the maximum errors with α = 0.3, 0.5, 0.8, 1.0 by different grid sizes. From this table, we can see that the LB model with double distribution functions is efficient for the fractional advection–diffusion problem coupled with incompressible N–S equation. From the table, we can conclude that the maximum errors E = ∥ϕ − ϕLB ∥∞ decrease as the lattice sizes ∆x decrease. Besides, it is obvious that present LB method model is convergent and the convergency rate is approximately 2. Furthermore, we use the fast solver [19] for the fractional derivative in Caputo sense. For the grid size 24 × 300 with α = 0.8 at t = 0.5, the CPU time is about 112 s while the CPU time of the algorithm without fast solver is about 11807s on the computer with iI7-7700HQ CPU. We can see that computation time of the history part is very big and the LB model with fast solver is very efficient and is easy to extend to 3D problem. 4. Conclusion In this paper, we construct a D2Q9 lattice Boltzmann model with double distribution functions for fractional advection–diffusion equations coupled with impressible Navier–Stokes equation. Through the Chapman–Enskog analysis, it has been found out that the macro controlling equations can be recovered correctly from our present LB model. Numerical simulations have been carried out. From the numerical results, it has been found out the LB solutions are agreement with the analytical solutions which demonstrate the LB model is accuracy and efficiency.

R. Du and Z. Liu / Applied Mathematics Letters 101 (2020) 106074

6

Table 1 Maximum errors at t = 0.1 with different lattice sizes and different α. α

Nx × Ny

E

0.3

8 × 100 16 × 200 24 × 300

1.4566 × 10−4 4.8625 × 10−5 2.4746 × 10−5

0.5

8 × 100 16 × 200 24 × 300

1.1358 × 10−4 2.6943 × 10−5 1.2531 × 10−5

0.8

8 × 100 16 × 200 24 × 300

2.7342 × 10−4 5.6372 × 10−5 2.2509 × 10−5

1.0

8 × 100 16 × 200 24 × 300

5.5165 × 10−4 1.7941 × 10−4 6.7785 × 10−5

Acknowledgments This work is financially supported by the National Natural Science Foundation of China (Grant No. 11602057). References [1] K.S. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, John Wiley and Sons Inc., New York, 1993. [2] F. Zeng, C. Li, F. Liu, I. Turner, Numerical algorithms for time-fractional subdiffusion equation with second-order accuracy, SIAM J. Sci. Comput. 37 (2015) 55. [3] S.L. Wu, T. Zhou, Parareal algorithms with local time-integrators for time fractional differential equations, J. Comput. Phys. 358 (2018) 135. [4] M.M. Meerschaert, C. Tadjeran, An explicit finite difference method and a new Von Neumann-type stability analysis for fractional diffusion-wave equations, J. Comput. Appl. Math. 172 (2004) 65. [5] S.B. Yuste, L. Acedo, An explicit finite difference method and a new Von Neumann-type stability analysis for fractional diffusion-wave equations, SIAM J. Numer. Anal. 42 (5) (2005) 1862. [6] C.P. Li, H.F. Ding, Higher order finite difference method for the reaction and anomalous-diffusion equation, Appl. Math. Model. 38 (2014) 3802. [7] R. Du, Z.P. Hao, Z.Z. Sun, Lubich second-order methods for distributed-order time-fractional differential equations with smooth solutions, East Asian J. Appl. Math. 6 (2016) 131. [8] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and beyond, Oxford University Press, Oxford, 2001. [9] W.Q. Hu, S.L. Jia, General propagation lattice Boltzmann model for variable-coefficient non-isospectral KdV equation, Appl. Math. Lett. 91 (2019) 61. [10] Z.H. Chai, N.Z. He, Z.L. Guo, B.C. Shi, Lattice Boltzmann model for high-order nonlinear partial differential equations, Phys. Rev. E 97 (2018) 013304. [11] H. Liang, Y. Li, J.X. Chen, J.R. Xu, Axisymmetric lattice Boltzmann model for multiphase flows with large density ratio, Int. J. Heat Mass Transfer 130 (2019) 1189. [12] H.L. Lai, C.F. Ma, Numerical study of the nonlinear combined sine-cosine-gordon equation with the lattice Boltzmann method, J. Sci. Comput. 53 (2012) 569. [13] Y. Xia, J.C. Wu, Y. Zhang, Numerical study of the nonlinear combined sine-cosine-gordon equation with the lattice Boltzmann method, Eng. Appl. Comput. Fluid Mech. 6 (2012) 581. [14] R. Du, D.K. Sun, B.C. Shi, Z.H. Chai, Lattice Boltzmann model for time sub-diffusion equation in Caputo sense, Appl. Math. Comput. 358 (2019) 80. [15] J.G. Zhou, P.M. Haygarth, P.J.A. Macleod, et al., Lattice Boltzmann method for the fractional advection-diffusion equation, Phys. Rev. E 93 (2016) 043310. [16] J.Y. Zhang, G.W. Yan, Lattice Boltzmann method for the fractional sub-diffusion equation, Internat. J. Numer. Methods Fluids 80 (8) (2015) 490. [17] Z.L. Guo, C.Z. Zheng, B.C. Shi, Lattice BGK model for incompressible Navier–Stokes equation, J. Comput. Phys. 165 (2000) 288. [18] Z. Guo, C.G. Zheng, B.C. Shi, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chin. Phys. 11 (2002) 366. [19] S.D. Jiang, J.W. Zhang, Q. Zhang, Z.M. Zhang, Fast evaluation of the Caputo fractional derivative and its applications to fractional diffusion equations, Commun. Comput. Phys. 21 (2017) 650.