Mesoscopic Simulation of the (2 + 1)-Dimensional Wave Equation with Nonlinear Damping and Source Terms Using the Lattice Boltzmann BGK Model.

In this work, we develop a mesoscopic lattice Boltzmann Bhatnagar-Gross-Krook (BGK) model to solve (2 + 1)-dimensional wave equation with the nonlinear damping and source terms. Through the Chapman-Enskog multiscale expansion, the macroscopic governing evolution equation can be obtained accurately by choosing appropriate local equilibrium distribution functions. We validate the present mesoscopic model by some related issues where the exact solution is known. It turned out that the numerical solution is in very good agreement with exact one, which shows that the present mesoscopic model is pretty valid, and can be used to solve more similar nonlinear wave equations with nonlinear damping and source terms, and predict and enrich the internal mechanism of nonlinearity and complexity in nonlinear dynamic phenomenon.

1. Introduction

Nonlinear dynamic phenomenon, which exists in many fields of science and engineering, such as hydrodynamic, nonlinear optics, biology, plasma physics, and so on, can be modeled by many systems of nonlinear partial differential equations (NPDEs) [1,2]. The dynamical processes of these nonlinear systems are very important for both production and scientific research, and they should be studied by a suitable method designed to treat the nonlinear problems. Many researchers use different analytical or numerical methods to investigate various nonlinear dynamic systems. Because of the complexity and particularity of the nonlinear evolution equations, there is no unity approach to find every solution to the nonlinear dynamic systems. Consequently, how to construct accurate and available methods to solve the nonlinear evolution equations has been an absorbing research career. In recent decades, with the vigorous development of computer science and technology, researchers have developed many different types of numerical methods to obtain numerical solutions, including the finite element, finite difference, finite volume, variational iteration, and spectral methods, etc. [3].

In this work, the generalized (2 + 1)-dimensional dynamical equation with nonlinear damping and source terms is considered to be as follows: where . The initial conditions associated with Equation (1) are given as follows: and the Dirichlet boundary conditions are given by where is the scalar variable; t is the time; is the Laplace operator; the parameter are supposed to be real number with . is the alleged dissipative term. When , Equation (1) is degraded into the undamped wave equation, while , to the damped one. The known functions and represent wave kinks or modes and velocity, respectively. and are known functions of their arguments. In recent years, many scholars used different types of methods to obtain the numerical solution, such as the implicit Lie-group iterative scheme [4], the meshless method [5], the space-time spectral method [6], the compact finite difference method [7], the nonconforming quadrilateral finite element method [8].

Recently, the mesoscopic lattice Boltzmann method (LBM) has made significant progress in the research nonlinear dynamical equations and evolving process of complexity micro-mesoscopic systems [9], especially in fluid mechanics [10,11,12]. Unlike traditional macroscopic numerical methods, which are independent of the discrete macroscopic evolution equations, the LBM is based on the mesoscopic kinetic evolution equations with distribution functions when the expression of the equilibrium distribution function is known. The fundamental idea is to take the place of the differential evolution equations of nonlinear system by the discrete kinetic Boltzmann equations. To obtain the macroscopic fluid behavior, we just need to calculate the discrete Boltzmann equations to obtain the evolution of the distribution function. From a computational resource perspective, the remarkable merits are brevity of programming, numerical potency, inherent parallelism, and ease treatment of intricate boundary conditions. This kind of method has comprehensive capacities in quite several fields, from phonon transport [13] to approximate incompressible flows [14,15,16,17,18,19,20,21,22,23,24,25], full compressible flows [26,27,28,29,30,31,32,33,34,35,36,37], dendrite growth [38,39] and thermal multiphase flows [40]. Recently, the mesoscopic kinetics method is also becoming increasingly popular in computational mathematics and engineering science for solving certain NPDEs, including Burgers’ equations [41,42], Korteweg-de Vries equation [43], Gross-Pitaevskii equation [44], convection-diffusion equation [45,46,47,48,49,50,51], Kuramoto-Sivashinsky equation [52], wave equation [53,54], Dirac equation [55], Poisson equation [56] etc.

Inspired by the successful promotion and application of the mesoscopic LBM in modeling nonlinear convection-diffusion system [45,46], the aim of this work is to further develop and apply the lattice Boltzmann Bhatnagar-Gross-Krook (BGK) method to solve (2 + 1)-dimensional wave equation with nonlinear damping and source terms. In the process of linking the mesoscopic Boltzmann equation to the nonlinear damped evolution system, we should choose suitable local equilibrium distribution functions to meet some constraints.

The content of this work is as follows: the next section presents the mesoscopic Boltzmann BGK model and deduces the wave equation with nonlinear damping and source terms through the multiscale expansion technique. Numerical verification of the model is presented in Section 3. Finally, a summary of the research is given in the last section.

2. Lattice Boltzmann BGK Model

In the present lattice Boltzmann BGK model, we use a single relaxation factor model for collision terms in this work. The discrete Boltzmann equation of the model with the BGK model takes the form [45] where is the dimensionless relaxation time which regulates the rate of access to equilibrium state, and are the distribution function and local equilibrium distribution function, respectively, and is the distribution function for the source term. is the collection of discrete directions of the particle velocity, for model, , , is the spatial step, is the time step.

In contrast to the common Lattice BGK model, we define the first derivative of as the following conservation condition [53]

To obtain the corresponding macroscopic evolution equation exactly, we should take as where the item is the unit tensor, the item is referring to the sound speed, satisfy . are the weights coefficients and satisfy the following conditions: , , , then .

Meanwhile, the corresponding source term is taken as where

Then, and should satisfy the following conservation conditions:

To obtain the macroscopic evolution Equation (1), we apply the Chapman-Enskog multiscale expansion to the distribution function, the first order time derivative, the spatial derivative and the source term as follows: where the item is as a small Knudsen number. Then, from Equation (7), and according to Equation (8), we obtain where .

Employing the Taylor formula to discrete Boltzmann Equation (2) at point (,t), we have where . Substitute Equation (10) into Equation (12), we have where .

Then we derive the first- and second-order equation in as

Multiplying both sides of Equation (14) by the operator , we obtain then substitute Equation (16) into Equation (15), we get

Summing Equations (14) and (17) over i and using Equations (5), (9) and (11), we get

Based on Equation (14), and using Equations (9) and (11), we have

Then, substituting Equation (20) into Equation (19), we obtain

Therefore, when Equation (18) + Equation (21) is applied, we have

To recover Equation (1) to order , we only need to let so, we have

In the calculation process, to get , using Equation (5), and applying the difference scheme to the item , we have then we get

3. Numerical Simulation

In this section, to show the efficiency of the present Lattice BGK model, we give some relevant numerical examples with and without damping terms. In addition, in order to compare with the exact solutions, the efficiency of present model is been tested. We set up the initial condition of distribution function by setting to equal for all grid points at . In addition, the macroscopic quantity in Equation (1) is also initialized by the given initial condition. The traditional explicit difference scheme can be used for calculating , here we use the analytic expression. The non-equilibrium extrapolation of distribution function proposed by Guo et al. [57] is applied to handle the boundary conditions. The instructions for the detailed process of boundary treatment are basic and detailed below.

Notice that the distribution function can be decomposed into its equilibrium and non-equilibrium parts where and are the equilibrium and the non-equilibrium parts of .

Through the Chapman-Enskog multiscale analysis for the LBM, we can assume that . For better presentation, we assume is a boundary node, and is the nearest neighboring grid point of at a distance . Thus, the non-equilibrium part of the distribution at grid point can be given by

Notice that the grid point is the nearest neighboring grid point at a distance , then we have , and we obtain where can be obtained by Equation (4), and and can be got by Equation (6). Therefore, as long as the macroscopic quantity of the boundary is given, the equilibrium distribution function of the boundary can be obtained, and the distribution function of the boundary can be obtained according to the above non-equilibrium extrapolation Formula (29).

Before simulation and calculation, we need to determine the expression of and according to the source term F of the given macroscopic Equation (1), then we can obtain the concrete discrete expression (4).

Next, we will introduce the calculation procedures of the present model as follows:

Step 1: Based on the given conditions and of Equation (2), initialize by Equation (6).

Step 2: Initialize by from Step 1 in all grid points.

Step 3: Calculate of the inner points by the discrete Boltzmann Equation (4).

Step 4: Calculate by Equation (26) and by Equation (5) of the inner points. If the specified termination time is reached, the program stops.

Step 5: Calculate and of the boundary points by the given conditions (2).

Step 6: Calculate of the boundary points by Equation (6).

Step 7: Calculate of the boundary points by Equation (29).

Step 8: Calculate of all grid points by Equation (4), then return to Step 4.

With the present mesoscopic model, we simulate several known exact solutions of the second-order (2 + 1)-dimensional hyperbolic telegraph equation and the (2 + 1)-dimensional damped, driven sine-Gordon equation, respectively. Furthermore, the (2 + 1)-dimensional undamped sine-Gordon equation with different initial condition of various ring solitons are studied to understand the nonlinear behavior characteristics of the system.

Meanwhile, we adopt four different kinds of error norms for measuring the present model’s precision. The root mean square error norm , max error norm , global relative error norm and root mean square error norm are generally defined as

The relative error norm (-error)

The max error norm (-error)

The global relative error norm (GRE-error)

The root mean square error norm (RMS-error)

Here, , are numerical solution and exact solution, respectively. The summation is added up from the information of all mesh points. Results show that the numerical solutions agree fairly well with the exact solutions over a considerable period of time.

Consider the following (2 + 1)-dimensional hyperbolic telegraph equation in the region , as follows: the initial conditions are given below

The exact solution for the current problem is given in Ref. [5] by

The boundary conditions are given from the exact solution.

In numerical simulation, we take , , , , . The computational domain is pinned to . We present the surface graph of the numerical and exact solutions by the present lattice BGK model at , see Figure 1. For clarity of contrast, we also demonstrate the two-dimensional contrast diagrams at for specific different times: , , and , see Figure 2. The relative error norm , the max error norm , the global relative error GRE norm and the root mean square error RMS norm for the solutions of the second-order hyperbolic telegraph equation at different instants of time can be found in Table 1.

Figure 1

Spatio-temporal evolution of the numerical (left) and exact (right) solutions at for Example 1.

Figure 2

Comparison of numerical solutions with exact ones at at different times for Example 1. The solid lines are drawn as the exact solutions.

Table 1

The maximum error , relative error norm, global relative error norm and root mean square error norm for at different times in Example 1.

t	1.0	2.0	3.0	4.0
L_2	1.3739×10^−4	1.0701×10^−4	3.9440×10^−5	1.5354×10^−4
L_∞	1.6459×10^−3	5.7798×10^−4	2.1405×10^−4	8.2911×10^−5
GRE	1.4147×10^−4	1.1327×10^−4	3.4975×10^−4	1.6388×10^−4
RMS	6.8635×10^−4	1.9667×10^−4	2.6666×10^−5	3.8189×10^−5

To measure the accuracy of the present mesoscopic model, the relative error, max error, global relative error and root mean square error are shown in Figure 3 at different time and with different resolutions, range from to and to 200, with . It is found that the present mesoscopic model is of second-order time accuracy. The order of the maximum error norm increases from 2.0720 to 2.2868, the order of the global relative error norm increases from 2.0789 to 2.2968, and the order of the root mean square error norm increases from 2.0720 to 2.2865.

Figure 3

Accuracy test at and with for Example 1.

Consider the following (2 + 1)-dimensional hyperbolic equation in the area , given by the initial conditions are given below

The exact solution for this instance is given in Ref. [6] by

The boundary conditions are given from the exact solution from Equation (39).

The term ’’ is a very good representation of nonlinearity. In the proceeding, we take , , , . The computational domain is pinned to . We present the spatio-temporal evolution of the numerical and exact solutions by the present model at , see Figure 4. For clarity of contrast, we also present the two-dimensional contrast diagrams at for especial different times: , , and , see Figure 5. The maximum value of the wave decays slowly over time due to the damping term and the source term. The relative error , max error norm , global relative error norm GRE and root mean square error norm RMS for the solutions of the second-order hyperbolic telegraph equation at specific times can be found in Table 2.

Figure 4

Spatio-temporal evolution of the numerical (left) and exact (right) solutions at for Example 2.

Figure 5

Comparison of numerical solutions with exact ones at at different instants of time for Example 2. The solid lines are drawn as the exact solutions.

Table 2

The maximum error norm , relative error norm , global relative error norm and root mean square error norm for at specific instants of time in Example 2.

t	1.0	2.0	3.0	4.0
L_2	6.6512×10^−4	1.1301×10^−3	1.7607×10^−3	2.3686×10^−3
L_∞	1.2614×10^−3	4.5286×10^−4	4.2965×10^−4	1.3513×10^−4
GRE	8.3720×10^−4	1.4311×10^−3	2.0109×10^−3	2.7853×10^−3
RMS	3.6249×10^−4	2.2658×10^−4	1.2987×10^−4	6.4271×10^−5

Consider the following (2 + 1)-dimensional hyperbolic equation in the area , given by the initial conditions are given below

The exact solution for this instance is given in Ref. [7] by

The boundary conditions are given from the exact solution.

In the proceeding, we take , , , , . The computational domain is pinned to . We present the spatio-temporal evolution of the numerical and exact solutions by the LBM at , see Figure 6. For clarity of contrast, we also present the two-dimensional contrast diagrams at for some specific times: , , and , see Figure 7. The relative error norm , max error norm , global relative error norm GRE and root mean square error norm RMS for the solutions of the second-order hyperbolic telegraph equation at specific times can be found in Table 3.

Figure 6

Spatio-temporal evolution of the numerical (left) and exact (right) solutions at for Example 3.

Figure 7

Comparison of numerical solutions with exact ones at at specific times for Example 3. The solid lines are drawn as the exact solutions.

Table 3

The maximum error norm , relative error norm , global relative error norm and root mean square error norm for at specific times in Example 3.

t	1.0	2.0	3.0	4.0
L_2	1.0885×10^−4	1.0395×10^−4	1.6103×10^−4	3.3706×10^−4
L_∞	4.8807×10^−5	2.7630×10^−5	1.8712×10^−5	2.0280×10^−5
GRE	1.0063×10^−4	9.2696×10^−5	1.7823×10^−4	3.8410×10^−4
RMS	2.1102×10^−5	8.9647×10^−6	6.9148×10^−6	7.9351×10^−6

Consider the following (2 + 1)-dimensional hyperbolic equation in the area , given by the initial conditions are given below

The exact solution for this instance is given in Ref. [7] by

The boundary conditions are given from the exact solution.

In the proceeding, we take , , , , . The computational domain is pinned to . We present the spatio-temporal evolution of the numerical and exact solutions by the present model at , see Figure 8. For clarity of contrast, we also present the two-dimensional contrast diagrams at for specific different times: , , and , see Figure 9. The relative error norm , max error norm , global relative error norm GRE and root mean square error norm RMS for the solutions of the second-order hyperbolic telegraph equation at specific times can be found in Table 4.

Figure 8

Spatio-temporal evolution of the numerical (left) and exact (right) solutions at for Example 4.

Figure 9

Comparison of numerical solutions with exact ones at at specific times for Example 4. The solid lines are drawn as the exact solutions.

Table 4

The maximum error norm , relative error norm , global relative error norm and root mean square error norm for at specific times in Example 4.

t	1.0	2.0	3.0	4.0
L_2	9.7057×10^−5	3.4711×10^−4	2.8500×10^−4	1.5080×10^−4
L_∞	5.9367×10^−5	1.2789×10^−4	6.3485×10^−5	2.0128×10^−5
GRE	9.6015×10^−5	3.4693×10^−4	2.8454×10^−4	1.5180×10^−4
RMS	2.9071×10^−5	6.3060×10^−5	3.1404×10^−5	1.0078×10^−5

Consider the following (2 + 1)-dimensional sine-Gordon equation [58] given by in the area , .

We simulate some particular cases of specific initial conditions with various numbers of circular ring solutions to study the nonlinear behaviors of the system. Numerical examples are carried out for three cases:

The first initial condition of one ring solitons is as follows: where .

The second initial condition of two ring solitons is the following: where , , .

The third initial condition of four ring solitons is as follows: where , , .

In our simulations, the zero gradient is used to deal with the boundary of the domain as , where is a boundary node, and is the nearest neighboring grid point of at a distance . The model parameters are set as , , . The other parameters are given as , , , . To better display the nonlinear propagation process of the wave, we present the surface graph of numerical solutions of collision of one ring solitons in regard to at , , , , , , see Figure 10. The surface graph of numerical solutions of collision of two ring solitons in regard to at , , , , , , see Figure 11. The surface graph of numerical solutions of collision of four ring solitons in regard to at , , , , , , see Figure 12.

Figure 10

Collision of one ring solitons: the function at , , , , , , for Example 5 with the initial condition (47).

Figure 11

Collision of two ring solitons: the function at , , , , , , for Example 5 with the initial condition (48).

Figure 12

Collision of four ring solitons: the function at , , , , , with the initial condition (49).

It can be found that the solitons reveal potent nonlinear evolution characteristics as time passed. One of the distinct features of solitons is that they can evolve without phanic changes in their identity after interplay. These numerical simulation results are in accordance with those results in Ref. [58].

4. Conclusions

Based on the mesoscopic lattice BGK method, we have investigated the numerical solution of (2 + 1)-dimensional wave equation with nonlinear damping and source terms, such as the hyperbolic telegraph equation, damped or undamped sine-Gordon equation, and so on. With the help of the Chapman-Enskog multiscale expansion, the macroscopic dynamical evolution equation can be precisely obtained from the present mesoscopic scheme in the continuity system without appending any amending term. Through observation, we can find that for the sine-Gordon system without damping terms and other source terms, the crest of the wave will oscillate up and down, and at the same time, the waveform will deform in the form of two or four crests that have evolved into one over time. All these phenomena reflect the evolution characteristics of nonlinear systems. Numerical examples for some test issues have been held to check the present mesoscopic model. The numerical solutions are in well coincident with the exact ones. From the convergence research of the Example 1, it can be found that the present mesoscopic model has the second-order accuracy in time. It is believed that with this model, we can predict and enrich the characterization and description of the nonlinear behavior characteristics in complex nonlinear dynamic systems.