A Multilevel Finite Difference Scheme for One-Dimensional Burgers Equation Derived from the Lattice Boltzmann Method

Abstract

An explicit finite difference scheme for one-dimensional Burgers equation is derived from the lattice Boltzmann method. The system of the lattice Boltzmann equations for the distribution of the fictitious particles is rewritten as a three-level finite difference equation. The scheme is monotonic and satisfies maximum value principle; therefore, the stability is proved. Numerical solutions have been compared with the exact solutions reported in previous studies. The 𝐿₂, 𝐿_∞ and Root-Mean-Square (RMS) errors in the solutions show that the scheme is accurate and effective.

1. Introduction

The lattice Boltzmann method (LBM) has been introduced as a new computational tool for the study of fluid dynamics and systems governed by partial differential equations. It has made a rapid development in theory and application over the last couple of decades since its inception [1–4]. This method can be either regarded as an extension of the lattice gas automaton [5] or as a special discrete form of the Boltzmann equation for kinetic theory [6]. The lattice Boltzmann models can also be used as partial differential equation (PDE) solvers. By choosing appropriate collision operator or equilibrium distribution, the lattice Boltzmann model is able to recover the PDE of interest. Recently, it has been developed to simulate linear and nonlinear PDE such as Laplace equation [7], Poisson equation [8, 9], the shallow water equation [10], Burgers equation [11], Korteweg-de Vires equation [12], Wave equation [13, 14], reaction-diffusion equation [15, 16], and convection-diffusion equation [17, 18].

The numerical schemes based on the LBM are given as a system of two-level explicit difference equations composed of the distribution functions of fictitious particles for each direction in which the particles move. For one-dimensional advection-diffusion problems, Ancona [19] showed that the LB schemes with the velocity model D1Q2 which includes two velocities with speed 1 in opposite directions to each other can be rewritten as the DuFort-Frankel scheme [20] which is a second-order three-level difference scheme. This shows that the accuracy of the LB schemes based on the model D1Q2 is identical to that of the DuFort-Frankel scheme. Suga [21] have proposed a four-level explicit finite difference scheme for 1D diffusion equation which is derived from the lattice Boltzmann method with rest particles. The consistency analysis of the scheme shows that the two parameters which appear in the scheme, the relaxation parameter and the amount of rest particles, can be determined such that the scheme has the truncation error of fourth order. In spite of the vast and successful applications, the numerical stability of the method has not been well understood. For certain specific class of lattice Boltzmann methods, for example, solving for linear and nonlinear convective-diffusive equation, there are some convergence and stability results given by Elton et al. [22].

Many works have been developed on lattice Boltzmann method to the Burgers equation in one or higher dimension [23–25]. In those papers, the standard lattice Boltzmann method was used and the macroscopic quantities were computed by the distribution function. However, those models are suffered from the stability. In this paper, we derive a three-level difference scheme for 1D Burgers equation based on the model D1Q2 from the LB schemes. It is generally recognized that LBM is a finite difference scheme of Boltzmann equation that has higher-order discretization error. We develop this method with the point of view above, but, at the same time, we also regard the LBM with BGK model as finite difference method for macroscopic equation. We find such LB scheme is a three-level finite difference one, which is monotonic and satisfies maximum value principle; therefore, we complete the proof of stability.

The rest of the paper is organized as follows. Section 2 describes the LB scheme with the velocity model D1Q2 and derives the three-level finite difference scheme which is equivalent to the LB scheme. A stability analysis of the scheme is given in Section 3. In Section 4, numerical solutions are compared with exact solutions reported in previous studies. And the conclusions are given in the end.

2. The Three-Level Finite Difference Scheme for 1D Burgers Equation Based on the LB Schemes

The one-dimensional Burgers equation take the following form:𝜕𝑢𝜕𝑡+𝑢𝜕𝜕𝑥𝑢=𝜈𝜕²𝑢𝜕𝑥²,(2.1) with the initial condition 𝑢(𝑥,0)=𝑢₀(𝑥). Here, the viscous coefficient 𝜈=1/Re, Re is the Reynolds number. Historically, (2.1) was first introduced by Bateman [26] who gave its steady solutions. It was later treated by Burgers [27] as a mathematical model for turbulence and after whom such an equation is widely referred to as Burgers equation. For a small value of 𝜈, Burgers equation behaves merely as hyperbolic partial differential equation and the problem becomes very difficult to solve as a steep shock-like wave fronts developed.

2.1. The Lattice Boltzmann Scheme

According to the theory of the LBM, it consists of two steps: (1) streaming, where each particle moves to the nearest node in the direction of its velocity; (2) colliding, which occurs when particles arriving at a node interact and possibly change their velocity directions according to scattering rules. Fictitious particles are introduced at each of the mesh points 𝑥=𝑗Δ𝑥(𝑗=…,−2,−1,0,1,2,…), and they move with the velocity 𝑐_𝑖 determined by the D1Q2 model from 𝑥 to the neighboring mesh point which was shown in Figure 1. The lattice Boltzmann schemes are established on grids with two directions𝑐₁,𝑐₋₁=[−𝑐,𝑐],(2.2) where 𝑐=Δ𝑥/Δ𝑡 is the speed in the system. Let 𝑓_𝑖(𝑥,𝑡) denote the distribution function of the particles moving with velocity 𝑐_𝑖. So the time evolution of the distribution function 𝑓_𝑖(𝑥,𝑡) is given by the following lattice Boltzmann equation (LBE) based on the Bhatnagar-Gross-Krook (BGK) model:𝑓_𝑖𝑥+𝑐_𝑖Δ𝑡,𝑡+Δ𝑡=𝑓_𝑖(𝑥,𝑡)−1𝜏𝑓^𝑖(𝑥,𝑡)−𝑓^eq_𝑖(𝑥,𝑡),(2.3) where 𝑓^eq_𝑖(𝑥,𝑡) is the local equilibrium distribution function of particles and 𝜏 is the dimensionless relaxation time which controls the rate of approach to equilibrium. The change in the distribution function produced by the collision of particles is approximated by the second term on the right-hand side of (2.3). The macroscopic velocity 𝑢(𝑥,𝑡) is defined in terms of the distribution function as𝑢(𝑥,𝑡)=_𝑖𝑓_𝑖(𝑥,𝑡)=_𝑖𝑓^eq_𝑖(𝑥,𝑡).(2.4)

Figure 1 

D1Q2 model with two velocities in one dimension.

In this paper, 𝑓^eq_𝑖(𝑥,𝑡) are determined as to satisfy (2.4) and the following conditions:_𝑖𝑐_𝑖𝑓^eq_𝑖(𝑥,𝑡)=𝑢²(𝑥,𝑡)_𝑖2,𝑐_𝑖𝑐_𝑖𝑓^eq_𝑖(𝑥,𝑡)=𝑐²𝑢(𝑥,𝑡).(2.5)

Solving these equations determines the equilibrium distribution functions𝑓^eq₁(𝑥,𝑡)=𝑢(𝑥,𝑡)2+𝑢²(𝑥,𝑡)Δ𝑡𝑓4Δ𝑥_eq−1,(𝑥,𝑡)=𝑢(𝑥,𝑡)2−𝑢²(𝑥,𝑡)Δ𝑡4Δ𝑥.(2.6) Applying the Chapman-Enskog expansion [24] yields the above Burgers equation (2.1) from the LBE and the equilibrium distribution functions given by (2.3) and (2.6), respectively. The viscosity 𝜈 is defined by 𝜈=(𝜏−1/2)(Δ𝑥²/Δ𝑡).

2.2. The Multilevel Finite Difference Scheme

Now, we let 𝑓^{𝑛𝑖,𝑗} denote 𝑓_𝑖(𝑗Δ𝑥,𝑛Δ𝑡) and let 𝑢^𝑛𝑗 denote 𝑢(𝑗Δ𝑥,𝑛Δ𝑡). We note that the subscript 𝑖,𝑗 combines information about the channel or direction of propagation (𝑖=1,−1) and location (𝑗 denotes a grid node). Using the equilibrium distribution function (2.6), the lattice Boltzmann equation (2.3) can be rewritten by classical finite different notation𝑓^{𝑛+11,𝑗+1}=1−1𝜏𝑓^𝑛1,𝑗+12𝜏𝑢^𝑛𝑗+Δ𝑡4𝜏Δ𝑥𝑢^𝑛𝑗²𝑓,(2.7)_{𝑛+1−1,𝑗−1}=1−1𝜏𝑓^{𝑛−1,𝑗}+12𝜏𝑢^𝑛𝑗−Δ𝑡4𝜏Δ𝑥𝑢^𝑛𝑗².(2.8)

According to (2.4), the macroscopic velocity can be computed by𝑢^𝑗𝑛+1=𝑓^𝑛+11,𝑗+𝑓^{𝑛+1−1,𝑗}=1−1𝜏𝑓^{𝑛1,𝑗−1}+1+𝑓^{𝑛−1,𝑗+1}2𝜏𝑢^𝑛𝑗−1+𝑢^𝑛𝑗+1+Δ𝑡4𝜏Δ𝑥𝑢^𝑛𝑗−1²−𝑢^𝑛𝑗+1²=𝐻𝑓_{𝑛1,𝑗−1},𝑓^{𝑛−1,𝑗+1},𝑢^𝑛𝑗−1,𝑢^𝑛𝑗+1.(2.9)

In addition,𝑓^{𝑛1,𝑗−1}+𝑓^{𝑛−1,𝑗+1}==𝑢𝑢^𝑛𝑗−1_𝑛𝑗−1−𝑓^{𝑛−1,𝑗−1}+𝑢^𝑛𝑗+1−𝑓^𝑛1,𝑗+1+𝑢^𝑛𝑗+1−𝑓^𝑛1,𝑗+1+𝑓^{𝑛−1,𝑗−1},(2.10) while𝑓^𝑛1,𝑗+1+𝑓^{𝑛−1,𝑗−1}=1−1𝜏𝑓^{𝑛−11,𝑗}+12𝜏𝑢^𝑗𝑛−1+Δ𝑡4𝜏Δ𝑥𝑢^𝑗𝑛−1²+1−1𝜏𝑓^{𝑛−1−1,𝑗}+12𝜏𝑢^𝑗𝑛−1−Δ𝑡4𝜏Δ𝑥𝑢^𝑗𝑛−1²=1−1𝜏𝑓^{𝑛−11,𝑗}+𝑓^{𝑛−1−1,𝑗}+1𝜏𝑢^𝑗𝑛−1=1−1𝜏𝑢^𝑗𝑛−1+1𝜏𝑢^𝑗𝑛−1=𝑢^𝑗𝑛−1.(2.11)

Then, (2.10) becomes𝑓^{𝑛1,𝑗−1}+𝑓^{𝑛−1,𝑗+1}=𝑢^𝑛𝑗−1+𝑢^𝑛𝑗+1−𝑢^𝑗𝑛−1.(2.12)

Substitute (2.12) to (2.9), we finally obtain the following three-level explicit finite difference scheme𝑢^𝑗𝑛+1=1−1𝜏𝑢^𝑛𝑗−1+𝑢^𝑛𝑗+1−𝑢^𝑗𝑛−1+12𝜏𝑢^𝑛𝑗−1+𝑢^𝑛𝑗+1+Δ𝑡4𝜏Δ𝑥𝑢^𝑛𝑗−1²−𝑢^𝑛𝑗+1².(2.13)

3. Stability Analysis

In this section, assumed the initial value 𝑢₀(𝑥) is bounded and smooth enough, we will prove the multilevel finite difference scheme is stable in 𝐿¹⋂𝐿^∞ space. Suppose𝑢₀(𝑥)∈𝐿¹,||𝑢₀(𝑥)||≤1.(3.1)

It is not difficult to see that, if |𝑢^𝑛𝑗|≤1 and𝜏≥1,Δ𝑡Δ𝑥≤1,(3.2) then the scheme (2.9) is monotonic increase. 𝜏≥1 means𝜈Δ𝑡Δ𝑥²≥12.(3.3)

Now, we will point out that the solution of the scheme (2.13) satisfies the maximum value principle.

Lemma 3.1 (maximum value principle). If initial value |𝑢₀(𝑥)|≤1 and the restrictions (3.2) hold, then, for all 𝑗∈𝑍, there are min_𝑙𝑢^0𝑙≤𝑢^𝑗𝑛+1≤max_𝑙𝑢^0𝑙,𝑛≥0.(3.4)

Proof. It is known that if we take 𝑓^01,𝑗=𝑢^0𝑗/2,𝑓^0−1,𝑗=𝑢^0𝑗/2, and 𝑢^𝑛𝐿=max_𝑗𝑢^𝑛𝑗,𝑢^𝑛𝑆=min_𝑗𝑢^𝑛𝑗𝑗∈𝑍, then, for all 𝑗,𝑘∈𝑍, 𝑓^11,𝑗+𝑓^1−1,𝑘==𝐻𝐻𝑓𝑢^01,𝑗−1_0𝑗−1,𝑓^0−1,𝑘+1,𝑢^0𝑗−1,𝑢^0𝑘+12,𝑢2^0𝑘+1,𝑢^0𝑗−1,𝑢^0𝑘+1≤𝐻𝑢_0𝐿2,𝑢2^0𝐿,𝑢^0𝐿,𝑢^0𝐿=1−1𝜏𝑢2^0𝐿+𝑢2^0𝐿+21𝜏𝑢^0𝐿+𝑢^0𝐿+Δ𝑡4𝜏Δ𝑥𝑢^0𝐿²=𝑢_0𝐿−𝑢^0𝐿²,(3.5) and similarly 𝑓^11,𝑗+𝑓^1−1,𝑘=𝐻𝑢^0𝑗−12,𝑢2^0𝑘+1,𝑢^0𝑗−1,𝑢^0𝑘+1≥𝐻𝑢_0𝑆2,𝑢2^0𝑆,𝑢^0𝑆,𝑢^0𝑆=𝑢_0𝑆.(3.6)
If we suppose 𝑢^0𝑆≤𝑓^𝑛1,𝑗+𝑓^{𝑛−1,𝑘}≤𝑢^0𝐿 is also correct. Particularly 𝑗=𝑘, we have 𝑢^0𝑆≤𝑢^𝑛𝑗≤𝑢^0𝐿, then 𝑓^𝑛+11,𝑗+𝑓^{𝑛+1−1,𝑘}≤=𝐻𝐻𝑓𝑓^{𝑛1,𝑗−1}_{𝑛1,𝑗−1},𝑓^{𝑛−1,𝑘+1},𝑢^𝑛𝑗−1,𝑢^𝑛𝑘+1=1−1,𝑓^{𝑛−1,𝑘+1},𝑢^0𝐿,𝑢^0𝐿𝜏𝑓^{𝑛1,𝑗−1}+𝑓^{𝑛−1,𝑘+1}+1𝜏𝑢^0𝐿≤𝑢^0𝐿.(3.7)
Similarly, we get 𝑓^𝑛+11,𝑗+𝑓^{𝑛+1−1,𝑘}≥𝑢^0𝑆.(3.8) Let 𝑗=𝑘, we can get min_𝑙𝑢^0𝑙≤𝑢^𝑗𝑛+1≤max_𝑙𝑢^0𝑙,𝑛≥0.(3.9)
Assume that 𝑢(̃𝑥,𝑡) is another solution of (2.1) with subject to initial condition 𝑢(̃𝑥,0)=𝑢₀̃(𝑥), and the initial condition satisfies |𝑢₀̃(𝑥)|≤1. Using the same scheme (2.13) and same restriction condition (3.2), we have the following.

Lemma 3.2. If the conditions of Lemma 3.1 are fulfilled, there are inequalities _𝑗_𝑗max𝑢^𝑗𝑛+1,𝑢̃^𝑗𝑛+1≤_𝑗max𝑢^0𝑗,𝑢^0𝑗̃,min𝑢^𝑗𝑛+1,𝑢̃^𝑗𝑛+1≥_𝑗min𝑢^0𝑗,𝑢^0𝑗̃.(3.10)
Denote that 𝑢^𝑛Δ𝑥={𝑢^𝑛𝑗,𝑗∈𝑍} is the discrete solution of LBE (2.7)–(2.9) at time 𝑛Δ𝑡, and ‖𝑢^𝑛Δ𝑥‖_𝐿¹=∑_𝑗|𝑢^𝑛𝑗|Δ𝑥 is the 𝐿¹ norm of discrete function 𝑢^𝑛Δ𝑥. Then, the solution is stable in the meaning of 𝐿¹.

Theorem 3.3. If 𝑢^𝑛Δ𝑥,𝑢̃^𝑛Δ𝑥 are the solutions of (2.13), 𝑢^0Δ𝑥,𝑢̃^0Δ𝑥∈𝐿¹(𝑅²) with subject to the corresponding initial conditions (3.1) and restrictions (3.2), then there are ‖‖𝑢^𝑛Δ𝑥−‖‖𝑢𝑢̃^𝑛Δ𝑥‖‖_𝐿1≤‖‖𝑢^0Δ𝑥−𝑢̃^0Δ𝑥‖‖_𝐿1,(3.11)_𝑛Δ𝑥‖‖_𝐿1≤‖‖𝑢^0Δ𝑥‖‖_𝐿1.(3.12)

Proof. Consider||𝑢^𝑗𝑛+1−𝑢̃^𝑗𝑛+1||=max𝑢^𝑗𝑛+1,𝑢̃^𝑗𝑛+1−min𝑢^𝑗𝑛+1,𝑢̃^𝑗𝑛+1.(3.13) Summing the absolute value to all 𝑗, by Lemma 3.2, we have _𝑗||𝑢^𝑗𝑛+1−𝑢̃^𝑗𝑛+1||=_𝑗max𝑢^𝑗𝑛+1,𝑢̃^𝑗𝑛+1≤_𝑗−_𝑗min𝑢^𝑗𝑛+1,𝑢̃^𝑗𝑛+1max𝑢^0𝑗,𝑢^0𝑗̃−_𝑗min𝑢^0𝑗,𝑢^0𝑗̃=_𝑗||𝑢^0𝑗−𝑢^0𝑗̃||.(3.14) If we let 𝑢̃_Δ𝑥(𝑥,𝑡)=0 in (3.11), we can get (3.12).

Remark 3.4. The restriction (3.2) is sufficient but not necessary.

4. Numerical Experiments

Example 4.1. We investigate the accuracy of the scheme by solving (2.1) on the domain (𝑡,𝑥)∈(0,𝑇]×[0,1]. The initial condition is 𝑢(𝑥,0)=sin(2𝜋𝑥),0≤𝑥≤1, and the homogenous boundary condition is 𝑢(0,𝑡)=𝑢(1,𝑡)=0. In this case, the exact Fourier solution is given by [28] 𝑢(𝑥,𝑡)=2𝜋𝜈∑^∞𝑛=1𝑎_𝑛exp−𝑛²𝜋²𝜈𝑡𝑛sin(𝑛𝜋𝑥)𝑎₀+∑^∞𝑛=1𝑎_𝑛exp−𝑛²𝜋²𝜈𝑡cos(𝑛𝜋𝑥),(4.1) where 𝑎𝑎₀=¹⁰exp−(2𝜋𝜈)⁻¹(1−cos(𝜋𝑥))𝑑𝑥,_𝑛=2¹⁰exp−(2𝜋𝜈)⁻¹(1−cos(𝜋𝑥))cos(𝑛𝜋𝑥)𝑑𝑥,𝑛=1,2,….(4.2)
In comparison with the analytical solutions, the efficiency of proposed model is validated. The following error norms are used to measure the accuracy:(1)𝐿₂-error ‖𝑒‖_𝐿2=^𝑛_𝑖=1𝑒^2𝑖^1/2,(4.3)(2)𝐿_∞-error ‖𝑒‖_𝐿∞=Max||𝑒_𝑖||,1≤𝑖≤𝑛,(4.4)(3) The root mean square (RMS) error ‖𝑒‖_RMS=^𝑛_𝑖=1𝑒^2𝑖𝑛^1/2.(4.5)
The numerical solutions of (2.1), which are computed by using different step size at time 𝑇=0.1 for 𝜈=1, are given in Table 1. The above error norms are given in Table 2 for different mesh size.
From Table 2, we find that the accuracy measured in 𝐿₂,𝐿_∞ and RMS norm errors increases as the step size decrease. The numerical solutions are in the symmetric pattern as the exact solutions are. Table 3 and Figure 1 show a comparison between numerical and exact solutions at different times for 𝜈=0.005. The curves for distribution of absolute errors at different times are also shown in Figure 2. It is known that the Fourier solutions for 𝜈≤0.001 fail to converge because of the slow convergence of the infinite series [28]. The numerical solution cures for 𝜈=0.001 at different time are drawn in Figure 3, which shows the correct physical behavior.

Table 1 

Comparison of the LB finite difference solutions with exact solution at 𝑇=0.1 for 𝜈=1 with 𝜏=1.

Table 2 

Error norms for 𝜈=1 at 𝑇=0.1 with different step size.

Table 3 

Comparison of the LB finite difference solutions with exact solution for 𝜈=0.005 with 𝑑𝑥=0.005,𝑑𝑡=0.003, and 𝜏=1.1 at different times.

(a) Numerical solutions

(b) Absolute errors

Figure 2 

Numerical solutions (a) and distribution of absolute errors (b) for 𝜈=0.005 at different times with 𝑑𝑥=0.005,𝜏=1.1, and 𝑑𝑡=0.003.

Figure 3 

Numerical solutions for 𝜈=0.001, at different times with 𝑑𝑥=0.001,𝜏=1 and 𝑑𝑡=0.0005.

Example 4.2. Consider Burgers equation with the following forms: 𝜕𝑢𝜕𝑡+𝑢𝜕𝜕𝑥𝑢=R1e𝜕²𝑢𝜕𝑥²,1𝑢(𝑥,0)=12≤𝑥≤23,𝑡>0,𝑢1Re𝑥+tan𝑥2,12≤𝑥≤32,𝑢32,𝑡=Re1+𝑡12+tan4(RRee+𝑡),𝑡>0,2,𝑡=Re1+𝑡32+tan4(R3eRe+𝑡),𝑡>0.(4.6)
It possesses the exact solution [23] 𝑢(𝑥,𝑡)=1Re+𝑡𝑥+tan2(R𝑥eRe+𝑡).(4.7)
In the computation, we compare the result with the D1Q2 and D1Q3 lattice Boltzmann model whose equilibrium distribution functions are taken as 𝑓^eq₁(𝑥,𝑡)=𝑢(𝑥,𝑡)2+𝑢²(𝑥,𝑡)𝑓4𝑐_eq2,(𝑥,𝑡)=𝑢(𝑥,𝑡)2−𝑢²(𝑥,𝑡)𝑓4𝑐_eq0,(𝑥,𝑡)=2𝑓3𝑢(𝑥,𝑡),_eq1(𝑥,𝑡)=𝑢(𝑥,𝑡)6+𝑢²(𝑥,𝑡)𝑓4𝑐_eq2,(𝑥,𝑡)=𝑢(𝑥,𝑡)6−𝑢²(𝑥,𝑡)4𝑐.(4.8)
Let Re=500, we give the results of our model, and exact solution as Figure 4 at 𝑡=0.4. Table 4 shows the results of the D1Q2, D1Q3, our model and the exact solution at different lattice at time 𝑡=0.4. The global relative errors GRE=∑^𝑖||𝑢^𝐸𝑥_𝑖,𝑡−𝑢^𝑁𝑥_𝑖,𝑡||∑_𝑖||𝑢^𝑁𝑥_𝑖,𝑡||,(4.9) which are used to measure the accuracy are presented in Table 5.
From Figure 4 and Table 4, we find that the D1Q2,D1Q3, and our model are all in excellent agreement with the exact solutions. The accuracy of the multilevel finite difference model is even higher than the D1Q2 and D1Q3 model. It should be pointed out that in order to attain better accuracy, the LB model requires a relatively small time step Δ𝑡 but the multilevel finite difference model does not have this restriction.

Table 4 

Comparison of the results with D1Q2, D1Q3, our model, and exact solution.

Table 5 

Global relative errors with different models.

Figure 4 

Comparison of the exact solution and our model. Parameters are: Re=500,𝑑𝑥=0.01,𝑑𝑡=0.002,𝜏=1.

5. Conclusion

In the current study, a three-level explicit finite difference scheme for 1D Burgers equation is derived by rewriting the LB scheme. Furthermore, it is proved that the scheme is conditionally stable. The efficiency and accuracy of the proposed scheme are validated through detail numerical simulation. It can be found that the numerical solutions are in excellent agreement with the analytical solutions. In order to derive LB scheme in a higher dimension, the LBM with the multispeed velocity model will be useful, in which different free parameters will be assigned for different values of the speed. Application of our method to 2D and 3D equations is left for future work.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (no. 51174236) and National Basic Research Program of China (2011CB606306).