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.


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].
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.

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, f j (x, t) and f eq j (x, t) are the distribution function and local equilibrium distribution function, respectively, and F j (x, t) is the distribution function for the source term. {c j , j = 0, 1, · · · , 8} is the collection of discrete directions of the particle velocity, for D2Q9 model, {c j , j = 0, 1, · · · , 8} = {(0, 0), (±c, 0), (0, ±c), (±c, ±c)}, c = ∆x/∆t, ∆x is the spatial step, ∆t is the time step.
In contrast to the common Lattice BGK model, we define the first derivative of u(x, t) as the following conservation condition [53] ∂u ∂t To obtain the corresponding macroscopic evolution equation exactly, we should take f eq j as where the item I is the unit tensor, the item c s is referring to the sound speed, satisfy c 2 s = c 2 /3. w j are the weights coefficients and satisfy the following conditions: ∑ j w j = 1, ∑ j w j c j = 0, ∑ j w j c j c j = c 2 s I, then w 0 = 4/9, w 1 4 = 1/9, w 5,8 = 1/36.
Meanwhile, the corresponding source term F j is taken as where Then, f eq j and F j 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 F (1) j = w j F (1) . Employing the Taylor formula to discrete Boltzmann Equation (2) at point (x,t), we have where where Then we derive the first-and second-order equation in ε as Multiplying both sides of Equation (14) by the operator D 1j , 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) ×ε 2 is applied, we have To recover Equation (1) to order O(ε 2 ), we only need to let so, we have In the calculation process, to get u(x, t), using Equation (5), and applying the difference scheme then we get

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 f j (x, t) by setting to equal f eq j (x, t) for all grid points at t = t 0 . In addition, the macroscopic quantity u(x, t) in Equation (1) is also initialized by the given initial condition. The traditional explicit ]/∆t can be used for calculating ∂ t F j (x, t), 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 f j can be decomposed into its equilibrium and non-equilibrium parts where f Through the Chapman-Enskog multiscale analysis for the LBM, we can assume that f j . For better presentation, we assume x b is a boundary node, and x f is the nearest neighboring grid point of x b at a distance c j ∆t. Thus, the non-equilibrium part of the distribution at grid point x f can be given by Notice that the grid point x f is the nearest neighboring grid point x b at a distance ∆x = εc j , , and we obtain where f j (x f , t) can be obtained by Equation (4), and f j (x f , t) 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 F j and ∂ t F j 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 u(x, 0) and u t (x, 0) of Equation (2), initialize f eq j (x, 0) by Equation (6).
Step 2: Initialize f j (x, 0) by f eq j (x, 0) from Step 1 in all grid points.
Step 3: Calculate f j (x, t) of the inner points by the discrete Boltzmann Equation (4).
Step 4: Calculate u(x, t) by Equation (26) and u t (x, t) by Equation (5) of the inner points. If the specified termination time is reached, the program stops.
Step 5: Calculate u(x, t) and u t (x, t) of the boundary points by the given conditions (2).
Step 6: Calculate f eq j (x, t) of the boundary points by Equation (6).
Step 7: Calculate f j (x, t) of the boundary points by Equation (29).
Step 8: Calculate f j (x, t + ∆t) 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 L 2 , max error norm L ∞ , global relative error norm GRE and root mean square error norm RMS are generally defined as (1) The relative error norm (L 2 -error) (2) The max error norm (L ∞ -error) (3) The global relative error norm (GRE-error) (4) The root mean square error norm (RMS-error) Here, u(x i , y j , t), u * (x i , y j , t) 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. Example 1. Consider the following (2 + 1)-dimensional hyperbolic telegraph equation in the region 0 ≤ x ≤ 2, 0 ≤ y ≤ 2 as follows: the initial conditions are given below The exact solution for the current problem is given in Ref. [5] by u(x, y, t) = e x+y−t .
The boundary conditions are given from the exact solution.
In numerical simulation, we take F(x, t) = −u − 2∂u/∂t − 2e x+y−t , α = 2, β = 1, ∆x = ∆y = 0.025, c = 200. The computational domain is pinned to Ω = [0, 2] × [0, 2]. We present the surface graph of the numerical and exact solutions by the present lattice BGK model at t = 4.0, see Figure 1. For clarity of contrast, we also demonstrate the two-dimensional contrast diagrams at x = 1.0 for specific different times: t = 1.0, t = 2.0, t = 3.0 and t = 4.0, see Figure 2. The relative error norm L 2 , the max error norm L ∞ , 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.   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 t = 1.0 and t = 2.0 with different resolutions, range from ∆t = 1.25 × 10 −4 to ∆t = 10 −3 and c = 25 to 200, with N x = N y = 80. 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.

Example 2.
Consider the following (2 + 1)-dimensional hyperbolic equation in the area 0 ≤ x ≤ 2, 0 ≤ y ≤ 2 given by 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 'sin u' is a very good representation of nonlinearity. In the proceeding, we take F(x, t) = −2 sin u + 2 sin[e −t (1 − cos(πx))(1 − cos(πy))] − π 2 e −t [cos(πx) + cos(πy) − 2 cos(πx) cos(πy)], We present the spatio-temporal evolution of the numerical and exact solutions by the present model at t = 4.0, see Figure 4. For clarity of contrast, we also present the two-dimensional contrast diagrams at x = 1.0 for especial different times: t = 1.0, t = 2.0, t = 3.0 and t = 4.0, 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 L 2 , max error norm L ∞ , 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.   Example 3. Consider the following (2 + 1)-dimensional hyperbolic equation in the area 0 ≤ x ≤ 1, 0 ≤ y ≤ 1 given by The exact solution for this instance is given in Ref. [7] by u(x, y, t) = e −(x+y)t sin(πx) sin(πy).
The boundary conditions are given from the exact solution.
In the proceeding, we take F(x, t) = −u − 2(1 + π 2 )∂u/∂t, α = 2(1 + π 2 ), β = 1, ∆x = ∆y = 0.0125, c = 500. The computational domain is pinned to Ω = [0, 1] × [0, 1]. We present the spatio-temporal evolution of the numerical and exact solutions by the LBM at t = 4.0, see Figure 6. For clarity of contrast, we also present the two-dimensional contrast diagrams at x = 0.5 for some specific times: t = 1.0, t = 2.0, t = 3.0 and t = 4.0, see Figure 7. The relative error norm L 2 , max error norm L ∞ , 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. Example 4. Consider the following (2 + 1)-dimensional hyperbolic equation in the area 0 ≤ x ≤ 1, 0 ≤ y ≤ 1 given by the initial conditions are given below  The exact solution for this instance is given in Ref. [7] by u(x, y, t) = e −0.5t sin(πx) sin(πy).
The boundary conditions are given from the exact solution.
In the proceeding, we take F(x, t) = −4u − 0.5∂u/∂t + 2(π 2 + 2)e −0.5t sin(πx) sin(πy), α = 0.5, β = 1, ∆x = ∆y = 0.0125, c = 500. The computational domain is pinned to Ω = [0, 1] × [0, 1]. We present the spatio-temporal evolution of the numerical and exact solutions by the present model at t = 4.0, see Figure 8. For clarity of contrast, we also present the two-dimensional contrast diagrams at x = 0.5 for specific different times: t = 1.0, t = 2.0, t = 3.0 and t = 4.0, see Figure 9. The relative error norm L 2 , max error norm L ∞ , 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.   Example 5. Consider the following (2 + 1)-dimensional sine-Gordon equation [58] given by in the area −a < x < a, −b < y < b.
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 −14 ≤ x, y ≤ 14.
The second initial condition of two ring solitons is the following: The third initial condition of four ring solitons is as follows: (3,17), (17,3), (17,17)}. In our simulations, the zero gradient is used to deal with the boundary of the domain as u(x b , t) = u(x f , t), where x b is a boundary node, and x f is the nearest neighboring grid point of x b at a distance c j ∆t. The model parameters are set as α = 4, β = 4.13, γ = 1/0.436. The other parameters are given as α = 1, β = 0, ∆x = ∆y = 0.05, c = 50. 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 sin(u/2) at t = 0, t = 4.0, t = 8.0, t = 12.0, t = 13.0, t = 15.0, see Figure 10. The surface graph of numerical solutions of collision of two ring solitons in regard to sin(u/2) at t = 0, t = 4.0, t = 8.0, t = 12.0, t = 13.0, t = 15.0, see Figure 11. The surface graph of numerical solutions of collision of four ring solitons in regard to sin(u/2) at t = 0, t = 2.5, t = 5.0, t = 7.5, t = 10.0, t = 12.5, see Figure 12.
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].

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.