From discrete to continuous Wardrop equilibria

1.   Introduction

In this paper, Backward Difference Lattice Boltzmann (BD-LB) method is used to simulate the Burgers equation with variable coefficient in space and time. It is very difficult to construct its analytical solution directly for the partial differential equation (PDE) with variable coefficient in space and time, whereas finding its numerical solution is an effective way. The traditional macroscopic numerical simulation methods include finite difference (FD) method, finite volume (FV) method, and finite element (FE) method. Compared with these methods, lattice Boltzmann method (LBM) is a mesoscopic numerical approach with many advantages, such as clear physical background, simple algorithm, easy parallel computing, and easy to implement program and handle complex boundary conditions ^[1,2,3].

The Boltzmann equation ^[4], a complex microintegral equation, is the fundamental equation of the gas kinetics theory. The right-end term is called the collision term, which is referred to as $ \Omega \left (f \right) $. The existence of $ \Omega \left (f \right) $ brings great difficulties to solve Eq (1.1).

However, the collision operators used in the LBM are generally based on the much simpler Bhatnagar-Gross-Krook (BGK) collision operator. Chen and Qian et al. ^[5,6] put forward to single-relaxation-time lattice Boltzmann (SRT-LB) model. The model controls the speed of different particles near the equilibrium state by using the same time relaxation coefficient.

where $ c_{j} $ is the discrete lattice velocity, $ \tau $ stands for the relaxation time, and $ f_{j}^ {eq} $ represents the distribution function of local equilibrium state.

Qin et al. ^[7] adds a body force term in Eq (1.2) to simulate the incompressible Navier-Stokes flow and then investigate aqueous humor dynamics in human eye. Hu and Lan et al. ^[8,9] add the modified function $ h_{i} $ to the Eq (1.2) to simulate the Gardner equation and KdV-Burgers equation with time-dependent variable coefficient respectively. Chai et al. ^[10,11] add the auxiliary distribution function $ G_{j}\left (\bf{x}, t \right) $ and a source term $ F_{j} \left (\bf{x}, t \right) $, proposing a multi-relaxed time lattice Boltzmann (MRT-LB) model to solve the Navier-Stokes equation and the convective diffusion equations. Gerasim V. Krivovichev ^[12] uses the Eq (1.2) formula to analyze the parameterized higher-order finite difference schemes of the linear advection equations. These schemes are based on a linear combinations of the spatial approximations of the convective term at the characteristic directions.

We use a FD method to solve the PDE problems in computational mathematics, but the forward difference format is often conditionally stable and the backward difference format is unconditionally stable. In this paper, we do not add correction function and auxiliary function, but we adopt the backward difference format and then rewrite the Eq (1.2) as follow:

The Burgers equation with variable coefficient is generally used to simulate the formation and decay of non-plane shock waves. The Burgers equation with the nonlinear term $ a(t) $ and the dispersion term $ b(t) $ can simulate the propagation of long shock waves in two shallow liquids ^[13].

It is shown in the literature that when $ b(t) $ is equal to 1, the shock wave solution reverses its velocity and collapses after $ a(t) $ changes the critical point of its symbol. In literature, the soliton-type solutions are constructed by using $ B\ddot{a}cklund $ transformation for the given form of $ a(t) $ and $ b(t) $. But these methods are very difficult for solving the variable coefficient partial differential equations, and we study the following Burgers equation with variable coefficient by using BD-LB method ^[14].

Equation (1.5) can represent a variety of physical models widely used in the fields of solid state materials, plasmas and fluid. $ u $ is amplitude function about time $ t $ and space $ x $. $ a(x, t) $ and $ b(x, t) $ are both analytical functions about $ x $ and $ t $. The subscripts represent partial derivatives.

This paper is organized as follows: In Section 2, the LB method for the Burgers equation of variable coefficient is recovered by using both the Chapman-Enskog (CE) analysis and the direct Taylor expansion (DTE) method. In Section 3, we first test Example 3.1 and obtain some reasonable lattice parameters. Then we apply these parameters to other variable coefficient examples and compare the numerical result to the analytical solution. In Section 4, we discuss the obtained results and draw a conclusion.

2.   Lattice Boltzmann model of nonlinear system

The Eq (1.3) is rewrote as follow:

We will use the following notations in this paper:

● $ f_{j}: = f_{j}\left(\boldsymbol{x}, t\right) $;

● $ f_{j, n-1}: = f_{j}\left(\boldsymbol{x}-c_{j} \delta t, t-\delta t\right) $;

● $ f_{j, n-1}^{e q}: = f_{j}^{e q}\left(\boldsymbol{x}-c_{j} \delta t, t-\delta t\right) $;

● $ f_{j, n-1}^{neq}: = \left(1-\frac{1}{\tau}\right)(f_{j, n-1}-f_{j, n-1}^{e q}) $;

● $ \delta $ is the Kronecker delta;

● $ \nabla $ is first-order spatial partial derivative.

In Eq (2.1), $ f_{j, n-1}^{e q} $ is the equilibrium distribution function, $ f_{j, n-1}^{neq} $ is the non-equilibrium distribution function. We define the macroscopic variable $ u $ as the sum of $ \sum_{i} f_{i} $.

Judging from Eq (2.1), we can get the following equation:

According to the equation proposed by Eq (1.2),

and then

Equation (2.4) can be written as follow:

we delimit $ f_{j, n+1} = e^{D_{j}}f_{j} $, where $ D_{j} $ is the difference operator $ D_{j} = \partial_{t}+c_{j} \nabla $.

The Eq (2.1) can be rewritten as follow:

2.1.   The Chapman-Enskog analysis

In the $ CE $ expansion, the LB equation is expanded by a dimensionless parameter $ \epsilon $, which is proportional the Knudsen number ($ Kn = \lambda /L $), $ \lambda $ is the mean free path and $ L $ is the feature length. We choose the local equilibrium distribution function in the following form ^{[2,10,15,16,17,18,20]}:

we adopt the weight family $ w_{j, a} $ selection satisfy

here, $ c_{S, a} $ is the lattice sound speed and the expansion of the function, $ r, s, t $ are expanded by $ CE $ as follows:

Correspondingly, the expansion of the equilibrium state distribution function is as follow:

that is to say, the equilibrium state distribution function varies when it is close to the scaling limit.

We define the second moment of $ u $ in Eq (2.7) as follows:

here, $ j > 0 $ breaks the conservation of $ u $ through the collision process because $ R_{n} \left (u \right) = r_{n} \left (u \right) $ and the $ r_{0} \left (u\right) = u $. In this paper, we only consider Eq (1.5) with a conservation $ u $, so we adopt $ r_{j} = 0 $ when $ j > 0 $. Furthermore, we hypothesize that $ u $ is differentiable in the analysis.

For $ \epsilon $, the temporal and spatial scale expansion as $ \partial_{t} = \sum_{k = 1}^{\infty} \epsilon^{k} \partial_{t_{k}} $, $ \nabla = \epsilon \nabla_{1} $, $ t_{k} $ is expressed as $ k $ time scales, $ \nabla_{k} $ is expressed as $ k $ space scales, respectively.

in which $ D_{j, k}: = \partial_{t_{k}}+\delta_{k, 1} c_{j} \nabla_{1} $. The solution of the Eq (2.6) can be written in the following form:

where

The Burgers equation has a second-order spatial derivative. It needs to sort the $ \epsilon^{2} $ by using the formalized $ u $, the results are summarized in Table 1 and we choose the following equation in order to make the results equivalent to the Burgers equation.

2.2.   The direct Taylor expansion method

We rewrite $ f_{j, n-1}^{neq} $ in Eq (2.1) into Eq (2.16) ^[19]:

Here, we assume $ f_{j}^{e q} $ as follow, and the equilibrium state distribution function $ f_{j}^{eq} $ is an analytic function.

the moments of $ w_{j} $ are shown in Eq (2.17) as in Table 2. Explicit forms for these weights are presented in Eqs (2.28) and (2.30).

In Eq (2.1), for $ j $ sum easy to calculate, on the left of the equation is $ u\left(x, t\right) $, and the equilibrium part to the right of Eq (2.1) changes to

At the same time, for $ f_{j, n-1}^{neq} $ in Eq (2.16), one obtain

where $ \mathcal{T}_{1} = \tau -1, \; \mathcal{T}_{2} = 2\tau^{2}-\tau-1, \; \mathcal{T}_{3} = 6\tau^{3}-6\tau^{2}+\tau-1 $ and $ \tau > 1/2 $ ^[20].

Putting Eqs (2.18) and (2.19) into Eq (2.1), we can get

In order to recover the macroscopic Burgers equation, the parameters are defined as

Compared with Eqs (2.15) and (2.21), the two analytical equations produce the same results for non-linear equations. These two methods are very different, but they recover the macroscopic Burgers equation.

In numerical simulation, we assume the physical space $ X = \delta x\cdot x $ and $ T = \delta t\cdot t $, and then put the $ X, \; T $ into Eq (2.20). Therefore, Eq (2.21) is as follows:

At the same time, the leading truncation error term at $ (\delta t)^{0} $ of Eq (2.20) is the fourth spatial derivative term whose coefficients involve $ \mathcal{K} $. This error term is

For the D1Q5 lattice, the format of $ f_{j}^{e q} $ is Eq (2.17). In order to remove the truncation error, the following $ \delta f_{j}^{eq} $ is added to $ f_{j}^{eq} $ in the D1Q7 lattice,

For the D1Q7 lattice, the following $ f_{j}^{e q} $ is employed:

We adopt the notation standardized in the LBM literature, where DdQq ^[6] refers to the d spatial dimensional model with q kinetic velocity. For the D1Q5 and D1Q7 models, a set of discrete weights having only the the unit $ n $ moment, is provided. These set of weights have the following properties:

In the case of D1Q5, when $ c_{i} = \left \{0, \pm 1, \pm 2 \right \} $, $ w_{i}^{(n)} $ can be obtained by inverting the matrix:

The $ w_{i}^{(n)} $ is

In the case of D1Q7, when $ c_{i} = \left \{0, \pm 1, \pm 2, \pm 3 \right \} $, $ w_{i}^{(n)} $ can be obtained by inverting the matrix:

The $ w_{i}^{(n)} $ is

3.   Numerical experiment

In this paper, the global relative error (GRE) is used to verify the effectiveness of the lattice Boltzmann model and the least-squares fitting is used to calculate the accuracy of the model.

among these, $ u\left (x, t \right) $ and $ u^ {\ast} \left (x, t \right) $ respectively represent the numerical solution and analytical solution.

Example 3.1. For $ a\left (x, t \right) = C_{1}, b\left (x, t \right) = C_{2} $, Eq (1.5) and the initial boundary conditions are taken as follows:

analytical solution ^[21]

We assume $ C_{1} = 1, \; C_{2} = 0.01, \; C_{3} = 2, \; T = 1 $. Figures 1–3 depict the images when $ \delta x = 0.05 $, $ \delta x = 0.1 $, $ \delta x = 0.2 $, and the left-hand plots in each figure show the minimum stability values of $ \tau $ and $ \delta t $ about the D1Q5 and D1Q7 lattices. We compare them with the reference ^[19], and the fitting about the images is very good. The value of $ \delta t $ is shown in Table 3, $ \delta t $ displays the maximum time increments to maintain stability with $ \tau = 1 $.

The plots on the right-hand in Figures 1–3 show the relationship of GRE and $ \delta t $ about D1Q5 and D1Q7 lattices. The crossed points show corresponding results between the D1Q5 and D1Q7 models. In Figure 4, the numerical results agree very well with the analytical solution for any t.

Example 3.2. For $ a\left (x, t \right) = sech^{2} \left (t \right), \; b\left (x, t \right) = C_{1} sech^{2} \left (t \right) $, Eq (1.5) becomes

analytical solution ^[22]

where $ k(x) = C_{2} x, \; c(t) = -C_{3} tanh\left (t \right) $, $ C_{2} $ and $ C_{3} $ are any constant, $ r^{2} \ge 1 $.

Here, macro parameters are assumed to be $ r = 2, \; l = 0, \; C_{1} = 2, \; C_{2} = C_{3} = 1 $ and lattice parameters are assumed $ \delta x = 0.2, \; \tau = 1 $, and the value of $ \delta t $ is shown in Table 3. The calculation area is fixed at $ \left [-10, 10 \right] $. Figure 5 shows the comparison plot of numerical results and analytical solution at different moments. The space-time plots of the numerical result from $ t = 0 $ to $ t = 10 $ are shown in Figure 6.

Remark 3.1. The colormap default in the MATLAB color box represents the solution of the D1Q5 simulation and the colormap jet stands for the solution of the D1Q7 simulation.

Then we do some numerical accuracy experiments. Several simulations are performed at different lattice resolutions $ \delta x = \left \{ 0.05, 0.1, 0.15, 0.2 \right \} $, and the value of $ \delta t $ is $ 1.0 \times 10^{-5} $. Based on the GRE at $ t = 1 $ and $ t = 2 $, the slopes of the fitting lines are very close to 2 in Figure 7, which indicates all of three models have a second-order accuracy in space. When $ \delta t = \left \{ 1.0 \times 10^{-5}, 5.0 \times 10^{-5}, 1.0 \times 10^{-4} \right \} $, $ \delta x = 0.2 $. Based on the GRE at $ x = 1 $ and $ x = 2 $, the slopes of the fitting lines are very close to 1 in Figure 8, which indicates all of three models have a first-order accuracy in time. The results are the same as Eq (2.20).

Example 3.3. For $ a\left (x, t \right) = 2sech^{2} \left (t \right) ln(cosh(\frac{x}{2})), \; b\left (x, t \right) = C_{1} sech^{2} \left (t \right) $. Equation (1.5) becomes

analytical solution ^[22]

where

When the analytical solution is positive, and the parameters are $ a_{0} = 1, \; r = 2 $, $ C_{1} = 2, \; C_{2} = C_{3} = 1 $, $ C_{4} = l = 0 $, lattice parameters are assumed to be $ \delta x = 0.2, \; \tau = 1 $, and the calculation area is fixed at $ \left[-10, 10 \right] $.

We present the comparison between detailed numerical results and analytical solution. Figure 9 shows the two-dimensional visual comparisons at some different times. The space-time evolution graph of the numerical results is shown in Figure 10. The numerical results show that the scheme has good long-time numeric simulation for the Burgers equation with variable coefficient in space and time. All of them clearly show that the numerical results agree with the analytical solutions well.

The calculation cost is given in Table 4 so as to compare the calculation time of the proposed scheme on D1Q7 and the basic scheme on D1Q5. We find that D1Q7 costs less than D1Q5. This improvement is even more significant as spatial step $ \delta x $ is increased.

4.   Conclusions

We derive the nonlinear Burgers equation with variable coefficient from the LB equation by using the CE analysis and the DTE methods. Based on comparative observations in Section 2, we suggest that when we derive LB models for macro equation, it is best to start with using the CE analysis of general equilibrium states. After obtaining some results on the equilibrium distribution function, we can apply DTE method to conduct error analysis to improve the stability and accuracy of the model.

In this study, we derive the LB model of Burgers by using a SRT-LB model and compare the LB solutions of the model with the corresponding analysis, which verifies the accuracy of the model. What's more, the improvement of accuracy by increasing lattice speeds can be regarded as a compensation for the deteriorated precision due to increased $ \tau $. Increasing the spatial step size $ \delta x $ can reduce the computational cost. In the future, we prepare to use LB model to simulate more non-linear PDEs with variable coefficient.

Acknowledgments

This study is supported by National Natural Science Foundation of China (Grant Nos. 11761005, 11861003), the Natural Science Foundation of Ningxia (2021AAC03206), Postgraduate Innovation Project of North Minzu University (YCX21156) and the First-Class Disciplines Foundation of Ningxia (Grant No. NXYLXK2017B09).

Conflict of interest

All authors declare no conflicts of interest in this paper.