A splitting lattice Boltzmann scheme for (2+1)-dimensional soliton solutions of the Kadomtsev-Petviashvili equation

1.   Introduction

The lattice Boltzmann (LB) method, originating from lattice gas automata (LGA), is a powerful computational fluid technique that has achieved much progress in the past three decades. In an LGA model, the flow of fluids is simulated by the movement and collision of a large number of indistinguishable particles on a regular spatial lattice. The Boolean variables are used to describe the particle distribution in the system. The density and velocity of the fluids are calculated through statistical processing of these Boolean variables. Although the LGA appears promising in simulating fluid problems, it has several disadvantages due to the iteration of Boolean variables. To improve on the LGA, in a general LB model, the ensemble-averaged particle distribution functions and the linearized Bhatnagar-Gross-Krook (BGK) approximation is applied to overcome the statistical noise and exponential complexity from which the LGA suffers ^[1,2,3]. Among the various numerical methods, the LB model has emerged as an impressive candidate because of its efficiency and simplicity in simulating multiphase flows ^[4,5,6,7,8], compressible flows ^[9,10], turbulent flows ^[11,12], combustion ^[13,14], non-Newtonian flows ^[15,16] and more.

In addition, by selecting specific forms of distribution functions, the LB method has been widely employed to solve partial differential equations (PDEs) numerically, such as wave equations ^[17,18,19], Korteweg-de Vries (KdV) equations ^[20,21,22], convection-diffusion equations ^[23,24,25], reaction-diffusion equations ^[26,27], and nonlinear Schrödinger equations ^[28,29].

In this research, the following (2+1)-dimensional Kadomtsev-Petviashvili (KP) equation will be investigated:

If the parameter is $\gamma = -3$, Eq (1) is called the Kadomtsev-Petviashvili-Ⅰ type (KPI) equation. Conversely, in the case $\gamma = 3$, Equation (1) is called the Kadomtsev-Petviashvili-Ⅱ type (KPII) equation. Both the KPI and KPII equations are fully integrable. The equation was first introduced by Kadomtsev and Petviashvili to analyze the stability of one-dimensional solitary pulses under transverse perturbations ^[30]. Regarded as an extension of the KdV equation in two-dimensional space, Eq (1) was developed for many application areas. For example, under the assumptions of long wavelengths and micro amplitudes, the water waves are governed by KP equations ^[30]. Bryant found that the propagation of weakly dispersive waves can be described with KP equations ^[31]. The two-dimensional shallow water waves with finite amplitude can be modeled from the KP equation ^[32]. Moreover, Xue found that several dust-acoustic solitons can be obtained by KP equations ^[33].

Given a wide range of applications in nonlinear science, KP equations have attracted much attention. To date, considerable theoretical studies and numerical simulations have been performed on KP equations. For example, Kricher et al. investigated the explicit solutions of KP equations ^[34,35]. Zhao proposed the extended mapping method to derive several exact solutions of KP equations and analyzed the interactions among different solitary waves ^[36]. Qin et al. introduced the extended homoclinic test method to study the breather wave and rogue wave solutions of the KP equations, which are important for describing several phenomena in plasma physics and fluid mechanics ^[37]. The Hirota bilinear method, which can be used to find the analytical solutions of PDEs ^[38,39], has been widely used to develop both rational and rogue wave solutions of KP equations ^{[40,41,42,43,44]}. The variable-coefficient symbolic computation method is carried out to obtain several multi-rogue wave solutions of the KP equation ^[45]. Using the generalized extended tanh method combined with the F-expansion method, several solitary wave solutions of the KP equations were exploited by Seadawy et al. ^[46]. The explicit finite difference model for the initial or boundary-value problems of KP equations were developed by Bratsos and Twizell ^[47]. To improve the stability of the algorithm, Feng and Mitsui utilized a linearized implicit finite difference method based on the Crank-Nicolson scheme to solve KP equations ^[48]. The pseudo-spectral method was carried out by Minzoni and Smyth to analyze the evolution of lump solitons of KP equations ^[49]. A computational approach based on the Adomian decomposition method was developed by Wazwaz to study the movement of solitons in KP equations ^[50]. A higher-order moment LB method was constructed by Wang to simulate single-line solitons and lump solitons in KP equations ^[51]. Cai et al. proposed a linearized local energy-preserving method to analyze the evolution of solitons in KP equations ^[52].

Notably, although the higher-order moment LB model presented by Wang successfully simulated the evolution of KP equations ^[51], the error of the model was unsatisfactory. The macroscopic quantity of the model was selected as Eq (2):

Thus, the true macroscopic quantity u must be obtained by the central difference in the x direction. This could cause error accumulation, which reduces the accuracy of the model. In this research, a more reliable numerical scheme for the KPI equation based on the LB model will be developed.

Inspired by ^[49], to achieve the aim of solving KPI equations more effectively, a splitting LB scheme is proposed in this paper. First, Eq (1) is translated as follows:

where $\delta = -\gamma$. We can see that Eq (3) is a KdV equation with the source term $Kw$.

It is reasonable to rewrite the KPI equation into Eqs (3) and (4). On one hand, as an extension of the KdV equation in two-dimensional space, the (2+1)-dimensional KP equation is derived from one-dimensional solitary pulses under transverse perturbations ^[30]. Equation (4) describes the second-order derivative term in the y-direction as transverse perturbations, which becomes the source term for the movement of solitons governed by the KdV equation in the x-direction. On the other hand, Beji ^[53] noted that in the nonlinear water wave equation regarding surface displacement, the KP equation can be obtained by taking the x-direction as the main propagation direction of waves, discarding the nonlinear and dispersive terms in the y-direction and retaining only the non-dispersive term in the y-direction. This point of view also verifies the rationality of Eqs (3) and (4).

The KdV equations with the source term were first calculated by Fornberg and Whitham ^[54]. In their work, the pseudo-spectral method combined with the leap-frog scheme was used. Meanwhile, the form of Eq (4) is similar to a diffusion-type equation. Although the characteristic of Eq (4) is different from the diffusion equation because the quantities $u$ and $w$ are not linearly related, Eq (4) can be integrated by using the Crank-Nicolson scheme.

In this paper, we focus on the LB solutions for these two equations. Our strategy is that two sets of distribution functions will be employed to recover these equations. A one-dimensional LB model is used to simulate the KdV equation with a source term, and a spatial evolution LB model is applied to develop (4).

The remainder of this research is constructed as follows. In Section 2, a splitting LB scheme for KPI equations is proposed. In Section 3, some detailed numerical simulations are conducted to verify the scheme. Finally, conclusions are provided in the last section.

2.   A splitting LB scheme for the KPI equation

2.1.   The LB model for the KdV equation with a source term.

A one-dimensional LB equation with a BGK collision term is considered as follows:

where the relaxation time is denoted as $\tau$. ${f}_{\alpha }\left(x, t\right)$ and ${f}_{\alpha }^{eq}\left(x, t\right)$ represent the discrete distribution function and equilibrium distribution function in the $\alpha$ direction, respectively. They should satisfy the following conservation law:

The additional term ${\omega }_{\alpha }\left(x, t\right)$ is used to simulate the source term of the macroscopic equation.

Now, applying the Chapman-Enskog expansion ^[55] to ${f}_{\alpha }\left(x, t\right)$, it yields:

Here, in Eq (7), ${f}_{\alpha }^{\left(0\right)}\left(x, t\right)\equiv {f}_{\alpha }^{eq}\left(x, t\right)$. The parameter $\varepsilon$ is called the Knudsen number. It is reasonable to assume that $\varepsilon$ is a small number. Therefore, we can take it as the time step: $\varepsilon = \delta t$ ^[56]. Thus, the LB equation becomes:

Since the Knudsen number $\varepsilon$ is small, we can employ the Taylor expansion to the left-hand side of the LB equation:

Moreover, the multiscale expansion to the additional term ${\omega }_{\alpha }\left(x, t\right)$ is obtained as follows:

To quantify changes in the LB equation at different time scales, the following is introduced:

and

Next, Eqs (7), (9), (10) and (12) are substituted into Eq (8), and the terms are separated based on the orders of $\varepsilon$:

In the above equations, $\Delta = \partial /\partial {t}_{0}+{e}_{\alpha }\cdot \partial /\partial x$ is a kind of partial differential operator. ${C}_{i}$ are several polynomials related to $\tau$:

To recover the KdV equation with the source term, $u\left(x, t\right)$ is considered as:

Due to Eq (6), the distribution functions take the following forms:

and

Now, the following functions of ${f}_{\alpha }^{\left(0\right)}\left(x, t\right)$ are selected:

In Eqs (23) and (24), the parameter $\chi = 1/{\varepsilon }^{2}{C}_{3}$ can be determined. These functions are equivalent to those in ^[20].

By taking a summation $\left(13\right)+\left(14\right)\times \varepsilon +\left(15\right)\times {\varepsilon }^{2}+\left(16\right)\times {\varepsilon }^{3}$ with respect to $\alpha$, we can obtain the following macroscopic equation:

This is the fourth-order accuracy KdV equation with a source term, in which the recovery of the source term $Kw$ depends on the selection of additional distribution functions ${\omega }_{\alpha }^{\left(i\right)}\left(x, t\right)$. More detailed information regarding this aspect can be found in Appendix A.

Equations (21)–(24) can be simply used to determine ${f}_{\alpha }^{eq}\left(x, t\right)$ by solving linear algebraic equations. For a one-dimensional five-bit lattice (D1Q5) where the discrete velocities are selected as ${e}_{\alpha } = \left\{{e}_{0}, {e}_{1}, {e}_{2}, {e}_{3}, {e}_{4}\right\} = \left\{0, c, -c, 2c, -2c\right\}$, the equilibrium distribution functions are:

2.2.   A spatial evolution LB model for Equation (4)

For the purpose of simulating Eq (4), we consider the following spatial evolution LB equation with the BGK collision term:

In Eq (31), $\widetilde{\tau }$ is a new relaxation time. $\delta x$ is the spatial step in the x-direction. ${g}_{\beta }\left(y, x\right)$ and ${g}_{\beta }^{eq}\left(y, x\right)$ are the distribution function and equilibrium distribution function at position $\left(y, x\right)$ with velocity ${e}_{\beta }$, respectively. They satisfy the following conservation condition:

Now, introducing a new dimensionless parameter $\widetilde{\varepsilon }$ that numerically equals $\delta x$. Then, Eq (31) can be rewritten as:

Next, the Taylor expansion

and the Chapman–Enskog expansion ^[55]:

can be adopted when $\widetilde{\varepsilon }$ is a small parameter. In (35), ${g}_{\beta }^{\left(0\right)}\left(y, x\right)\equiv {g}_{\beta }^{eq}\left(y, x\right)$. By denoting different spatial scales ${x}_{0}$, ${x}_{1}$, ..., they are

Substitute Eqs (34), (35) and (37) into Eq (33) and separate terms based on different orders of $\widetilde{\varepsilon }$:

where ${\widetilde{C}}_{2} = 0.5-\widetilde{\tau }$ and the partial differential operator is $\widetilde{\Delta } = \partial /\partial {x}_{0}+{e}_{\beta }\cdot \partial /\partial y$.

The macroscopic quantity $w\left(y, x\right)$ in Eq (4) is defined as

By reverting back to Eqs (32) and (35), we can obtain

and

Furthermore, several functions related to ${g}_{\beta }^{\left(0\right)}\left(y, x\right)$ are selected as:

By taking $\left(38\right)+\left(39\right)\times \widetilde{\varepsilon }$ and summing about $\beta$, we can obtain the following:

Due to Functions (43), (44) and Eq (38), Eq (4) with second-order accuracy of truncation error can be recovered:

For the D1Q5 lattice where the discrete velocities are selected as ${e}_{\beta } = \left\{{\widetilde{e}}_{0}, {\widetilde{e}}_{1}, {\widetilde{e}}_{2}, {\widetilde{e}}_{3}, {\widetilde{e}}_{4}\right\} = \left\{0, \widetilde{c}, -\widetilde{c}, 2\widetilde{c}, -2\widetilde{c}\right\}$, the equilibrium distribution functions can be derived by combining Eqs (41), (43) and (44). They are represented as follows:

Here, the parameter is $\lambda = -\delta /\left(K\widetilde{\varepsilon }{\widetilde{C}}_{2}\right)$.

3.   Numerical experiments

In this section, three numerical experiments are proposed to verify the model.

3.1.   Single-line soliton.

First, we consider the single-line soliton solution of the KPI equation:

The following initial condition is employed:

Equation (50) describes the propagation process of a single-line soliton wave. The propagation speed is $4{\xi }^{2}-3{\eta }^{2}$, and $\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(-{\eta }^{-1}\right)$ is the angle between the propagation direction and the positive x-axis ^[48].

The initial condition and the splitting LB result when $t = 2$ are shown in Figure 1. The parameters are as follows: $\delta t = 0.0005$, lattice size $200\times 100$, $\tau = 1.3$, $\widetilde{\tau } = 1.0$, $c = 200$, $\widetilde{c} = 1$, $\gamma = -3$, $\xi = 1.0$, $\eta = -1.0/\sqrt{2}$, ${x}_{0} = 4.0$, $K = 200.0$, $\delta x = 0.1$ and $\delta y = 0.1$. This numerical simulation is carried out on the domain $20\times 10$, and the exact Dirichlet boundary conditions are applied.

Additionally, with the parameters given in Figure 1, the comparison between the theoretical solutions and splitting LB results at $y = 5.0$ when $t = 1.0$ is presented in Figure 2(a), from which it can be seen that the splitting LB results are in good agreement with the exact solutions.

The relative errors at line $y = 5.0$ when $t = 1.0$ are also plotted in Figure 2(b). Here, the relative errors are defined as $e = \left|{u}^{N}-{u}^{E}\right|/\left|{u}^{E}\right|$, in which ${u}^{N}$ is the splitting LB solution and ${u}^{E}$ is the exact solution. Figure 2(b) illustrates that the relative errors are lower than $1.5\times {10}^{-5}$. The result demonstrates that the splitting LB solutions are acceptable.

The comparison of the general relative errors $G$ between the splitting LB model and the traditional higher-order moment LB model is presented in Table 1. Here, the general relative errors $G$ are defined as:

From Table 1, it can be clearly seen that the general relative errors of the splitting LB scheme slightly fluctuate but are consistently lower than $1.0\times {10}^{-5}$. However, the general relative errors of the traditional higher-order moment LB model are approximately ${10}^{-1}$, which is significantly higher than the $G$ in the splitting LB scheme. Furthermore, in the higher-order moment LB model, the general relative errors increase over time. This is due to the error accumulation of the model.

With time step $\delta t = 0.01$, other parameters are the same as those in LB models, the general relative errors of the linearized local energy-preserving method ^[52] are also presented in Table 1, from which we can observe that the errors of the splitting LB scheme are also lower than those of the linearized local energy-preserving method.

In addition, the CPU times of these three numerical methods are considered. Here, the experiment is conducted on a computer with an Intel i5 (2.50 GHz) processor, and the computing program is written in Fortran and compiled in Gfortran. The CPU times of the splitting LB scheme, traditional higher-order moment LB model and linearized local energy-preserving method are 33.73s, 36.97s and 22.19s, respectively. Here, we can see that due to the large time step, the linearized local energy-preserving method has the lowest CPU time. The CPU time of the splitting LB scheme is lower than that of the higher-order moment LB model.

From the discussion above we can see that the splitting LB scheme has an advantage over the traditional high-order moment LB model in simulating KPI equations.

For the purpose of testing the convergence order of the splitting LB scheme, the relation between the general relative errors $G$ and spatial step $\delta x$ is shown in Figure 3(a). To compute the numerical convergence order, the linear fit is also adopted. The linear fit line is ${\mathrm{l}\mathrm{o}\mathrm{g}}_{10}G = 2.5547\times {\mathrm{l}\mathrm{o}\mathrm{g}}_{10}\delta x$ $+0.03765$, which means that the model's convergence order is 2.5547 in space.

To quantify the order of the truncation error, we introduce another parameter $\rho = e/{\left(\delta x\right)}^{2}$. Here, $\delta x$ is the spatial step in the x-direction, and $e$ is the absolute error. Parameter $\rho$ represents the scale factor of the truncation error $e$ and the square of the spatial step ${\left(\delta x\right)}^{2}$. In Figure 3(b), the curve of $\rho$ at times from $t = 0.0$ to $t = 1.0$ is plotted. From the graph, we can observe that $\rho$ is a finite small number in the range of $\left[0.12, \mathrm{ }0.16\right]$ in this example, which indicates that the splitting LB scheme in this paper has a second-order accuracy of the truncation error.

3.2.   Interaction of two line solitons

We consider the KPI equation with the following initial condition:

This problem has the following theoretical solution:

The computational domain is $\left[\mathrm{0, 40}\right]\times \left[\mathrm{0, 40}\right]$, and the exact Dirichlet boundary conditions are used in this experiment.

In Figure 4, the splitting LB solutions at different times $t = 0$, $t = 3$, $t = 5$ and $t = 8$ are carried out. The parameters are as follows: $\delta t = 0.0005$, lattice size $200\times 100$, $\tau = 1.3$, $\widetilde{\tau } = 1.0$, $c = 400$, $\widetilde{c} = 2$, $\gamma = -3$, ${\xi }_{1} = 1.0$, ${\xi }_{2} = 1.0/\sqrt{2}$, ${\eta }_{1} = -1.0/\sqrt{3}$, ${\eta }_{2} = -1.0$, ${x}_{\mathrm{0, 1}} = 4.0$, ${x}_{\mathrm{0, 2}} = 14.0$, $K = 200.0$, $\delta x = 0.2$ and $\delta y = 0.4$. Figure 4(a) shows the initial condition Eq (53) in which two line solitons are allowed a collision to take place. Figure 4(b) demonstrates that two solitons collide when $t = 3$. In Figure 4(c), when $t = 5$, we find that the solitons collide at $y > 10$. Finally, in Figure 4(d), the two solitons have separated completely when $t = 8$. The simulation results in Figure 4 are consistent with the classical results ^[48].

3.3.   Lump soliton

To simulate lump soliton, the KPI equation with the following initial condition is given:

The theoretical solution to this problem is as follows:

which shows that the lump soliton will propagate along the x-direction at a speed of $3{\mu }^{2}$. In this case, the exact Dirichlet boundary conditions are considered with the computation domain: $\left[\mathrm{0, 20}\right]\times \left[\mathrm{0, 20}\right]$.

In Figure 5(a), the initial condition is shown. The parameters are as follows: $\delta t = 0.0005$, lattice size $200\times 100$, $\tau = 1.3$, $\widetilde{\tau } = 1.0$, $c = 200$, $\widetilde{c} = 2$, $\gamma = -3$, $\mu = 1.0$, ${x}_{0} = 10.0$, ${y}_{0} = 10.0$, $K = 200.0$, $\delta x = 0.1$ and $\delta y = 0.2$. Figures 5(b) and 5(c) show the splitting LB results at $t = 0.5$ and $t = 1.0$, respectively. These numerical results are consistent with those in ^[48].

4.   Conclusions

In this paper, a splitting LB scheme for solving KPI equations was presented. In this algorithm, the variable substitution method was used. As a result, the KPI equation has been written into two macroscopic equations. Two sets of distribution functions are employed for the construction of these two macroscopic equations. This splitting scheme contains both time evolution and spatial evolution in which the second equation is based on a space variable rather than time. Three numerical examples were conducted to show that this splitting LB scheme was more suitable for simulating KPI equations than the traditional higher-order moment LB model.

It was laborious to apply the traditional LB model to higher order PDEs due to the complexity of the higher-order Chapman-Enskog expansion. However, the splitting LB method appeared promising in simulating higher order PDEs because the splitting technique can reduce the order of the macroscopic equation. How the splitting LB method can be applied to these equations still requires further research.

Use of AI tools declaration

The author declares he has not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The research was funded by the Doctoral Research Project of Northeast Electric Power University (BSJXM-2021235).

Conflict of interest

The author declares no conflicts of interest.

Appendix A

In this section, the multiscale expansion technique is applied to recover the KdV equation with a source term. First, we selected the following summations of the additional distribution functions:

By taking Eq (13) and summing over $\alpha$, according to Eqs (19) and (21), the following relation can be derived:

The macroscopic Eq (4) and the summation Eq (A.1) can be modified to obtain the following:

Equations (A.5) and (A.6) are called conservation laws in the ${t}_{0}$ scale.

Similarly, we take Eq (14) and sum over $\alpha$. Based on Eqs (19), (21), (22) and (A.2), the following can be derived:

Because of the macroscopic Eq (4), we can obtain the following:

Equations (A.7) and (A.8) are called conservation laws in the ${t}_{1}$ scale.

To recover the macroscopic equation, we take $\left(13\right)+\varepsilon \times \left(14\right)+{\varepsilon }^{2}\times \left(15\right)$ and sum over $\alpha$. Thus, we obtain the following:

Due to Functions (21)–(22), the third term on the left-hand side of Eq (9) is equal to zero. To recover the term ${\partial }^{3}u/\partial {x}^{3}$ in the macroscopic equation, the function $P\left(x, t\right)$ should be selected in Eq (23). Substituting Eqs (A.3) and (A.6) into Eq (A.9) yields the following:

This is the KdV equation with a source term if $\chi = 1/{\varepsilon }^{2}{C}_{3}$.

Furthermore, to improve the accuracy of the equation, we take Eq (16) and sum it over $\alpha$ on account of the conservation law Eq (A.5), thus getting:

If we select:

then

Moreover, several terms that are not equal to zero are:

We assume that the additional distribution functions ${\omega }_{\alpha }^{\left(i\right)}\left(i = \mathrm{1, 2}, \mathrm{3, 4}\right)$ are independent of $\alpha$. Hence, for the D1Q5 lattice, the following can be obtained:

Because of the summation Eq (A.4), we derive:

Taking $\left(13\right)+\varepsilon \times \left(14\right)+{\varepsilon }^{2}\times \left(15\right)+{\varepsilon }^{3}\times \left(16\right)$ and summing over $\alpha$, the KdV equation with a source term and fourth-order accuracy of truncation error can finally be achieved:

From the summations Eqs (A.1)–(A.4) that we select, the additional distribution functions can finally be obtained: