# Finite Element Method of BBM-Burgers Equation with Dissipative Term Based on Adaptive Moving Mesh

Authors: Changna Lu; Qianqian Gao; Chen Fu; Hongwei Yang
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

A finite element model is proposed for the Benjamin-Bona-Mahony-Burgers (BBM-Burgers) equation with a high-order dissipative term; the scheme is based on adaptive moving meshes. The model can be applied to the equations with spatial-time mixed derivatives and high-order derivative terms. In this scheme, new variables are needed to make the equation become a coupled system, and then the linear finite element method is used to discretize the spatial derivative and the fifth-order Radau IIA method is used to discretize the time derivative. The simulations of 1D and 2D BBM-Burgers equations with high-order dissipative terms are presented in numerical examples. The numerical results show that the method keeps a second-order convergence in space and provides a smaller error than that based on the fixed mesh, which demonstrates the effectiveness and feasibility of the finite element method based on the moving mesh. We also study the effect of the dissipative terms with different coefficients in the equation; by numerical simulations, we find that the dissipative termuxx plays a more important role than uxxxx in dissipation.


---

## Body

## 1. Introduction

In mathematics, physics, and engineering, nonlinear partial differential equations play an important role due to their rich mathematical structures and features [1–4] as well as important applications in fluid dynamics, plasma physics, fiber optics, condensed matter physics, and chemistry [5–11]. It is important to search for solutions. These nonlinear partially differential equations can be solved exactly or numerically to understand the properties of solutions.In this paper, we consider the nonlinear multidimensional Benjamin-Bona-Mahony-Burgers (BBM-Burgers) equation with a dissipative term and use the finite element method based on the moving mesh to solve its numerical solution. The equation has the following form:(1)ut+∑i=1dαiuxi+∑i=1dβiuuxi-γΔu-μΔut+∑i=1d∑n=1pδn∂2n+2u∂xi2n+2=F,X∈Ω,t∈0,TuX,t=gX,t,X∈∂Ω,t∈0,TuX,0=u0X,X∈Ω,where X=(x1,x2,…,xd), d stands for the dimension, xi stands for the corresponding spatial vector, the coefficients αi,βi∈R, γ,μ,δn>0, Ω⊆Rn, Δ is a Laplace operator, and p is one constant to define the order of the dissipative term. In the equation, F is a sufficiently smooth function and stands for the source term, and g(X) and u0(X) are given functions.Equation (1) is related to the well-known Benjamin-Bona-Mahony (BBM) equation, which was advocated by Benjamin et al. [12] in 1972 as a refinement of the Korteweg-de-Vries (KdV) equation. There are several different transformations of this equation; when the coefficients μ=δn=0, F=0, the equation becomes ut+∑i=1dαiuxi+∑i=1dβiuuxi=0, which is the Burgers equation with a nonlinear term. When parameters γ,δn=0, (1) becomes ut+∑i=1dαiuxi+∑i=1dβiuuxi-μΔut=F, which is called Regularized Long Wave (RLW) equation, also called BBM equation. The RLW equation is usually used to analyze small amplitude long waves on the surface of water in fluid dynamics. When δn=0, (1) becomes a generalized BBM-Burgers equation, which is a mathematical model of propagation of small-amplitude long wave in nonlinear dispersive media and better than Korteweg-de-Vries-Burgers (KdV-B) equation in fluid mechanics; when F=0, d=1, and g(x)=0, the BBM-Burgers equation satisfies two conservation laws:(2)MassfunctionE1t=∫Ωudx,EnergyfunctionE2t=∫Ω12u2+μ2ux2dx+α∫0t∫Ωux2ςdxdς.If d=2, p=1, (1) becomes(3)ut+α1ux+α2uy+β1uux+β2uuy-γΔu-μΔut+δuxxxx+uyyyy=F,x,y∈Ω,t∈0,Tux,y,t=gx,y,t,x,y∈∂Ω,t∈0,Tux,y,0=u0x,y,x,y∈Ω. In physical phenomena, the effect of dispersion in the BBM-Burgers equation with a high-order dissipative term is similar to the RLW equation; the effect of dissipation is similar to the Burgers equation.To get the approximate solution of the BBM-Burgers equation with a high-order dissipative term, some numerical and analytical methods have been proposed in recent years. In the aspect of mathematical theory, Zhao et al. consider the existence of the generalized BBM-Burgers equation with a dissipative term and verify the corresponding convergence in H-measures [13, 14]. Kondo and Webler apply kinetic decomposition for the generalized BBM-Burgers equations with a dissipative term to obtain the approximate transport equation, and then take advantage of the averaging lemma to get the convergence [15]. Seok states the existence and conservation laws of the generalized BBM-Burgers equation with a dissipative term [16]. For more theoretical analyses, refer to the references in [13–16]. In the aspect of numerical solutions, many methods have been used, like Painlev’e test, Darboux Transformation, bilinear method, symmetry method [17], and so on [18–21]. Hong et al. apply the homotopy perturbation method and He’s variational iteration method to solve the BBM-Burgers equation. Ganji et al. put forward the notion that the method they used in [22] is better than the homotopy analysis method proposed in [23]. Finite element Galerkin methods have been discussed by Kadri et al. [24], Lee [25], Yin and Piao [26], and Dehghan et al. [27] to solve the BBM-Burgers equations with a dissipative term. Finite difference methods have also been widely used to solve the equations in [28–30]; as for other numerical methods, refer to [31–34].The references introducing the numerical simulations above are mostly based on fixed meshes; in this paper, we consider the finite element method based on the adaptive moving mesh to solve the BBM-Burgers equation with a high-order dissipative term. The moving mesh method is a kind of adaptive mesh method which adds elements near the place where the numerical solution changes rapidly and decreases elements in the solution changing slowly; the total element number still remains unchanged. When we compare the finite element method based on adaptive moving mesh with methods based on fixed mesh or proposed in [26], we find that the moving mesh method takes fewer elements to get the same distinguishability and simulates steep waves and the transition of gap distinctly. In recent years, the moving mesh method has witnessed further development and extensive application [35]. This method was applied to approximate the anisotropic diffusion-convection-reaction problems to satisfy a discrete maximum principle, and it was adopted to study the porous medium equation and approximate the solution [36]. For more information about the moving mesh method, please refer to [37–41].The remainder of this paper is divided into three sections. In Section2, we introduce the method we used, the discretization of the BBM-Burgers equation with a fourth-order dissipative term in one dimension, and the implementation of the moving mesh method. Section 3 presents 1D and 2D numerical examples of BBM-Burgers equations with a dissipative term. Finally, Section 4 shows the conclusion we draw from the numerical examples.


## 2. Finite Element Method Based on Moving Mesh

In this section, we mainly introduce the linear finite element method based on moving mesh and make use of the method to solve the BBM-Burgers equation with a high-order dissipative term. For convenience, in Section2.2, we introduce the finite element method of the one-dimensional equation with a fourth-order dissipative term uxxxx; the other parts of this section will focus on two dimensions.


### 2.1. Adaptive Moving Mesh Method

In this subsection, we mainly want to show the generation of physical meshThn+1 at a new time level by an adaptive moving mesh strategy based on the Moving Mesh Partial Differential Equation (MMPDE) [42]; at the current time level, the physical mesh Thn=Xjn and the numerical solution (uh)n are known (the mesh and value will be introduced clearly in the following parts). During the strategy, one of the crucially important rules is to view the nonuniform mesh (i.e., the physical mesh in this paper) as a uniform mesh under the metric tensor M; the tensor M is a symmetric and positive definite matrix. Under the metric tensor M, the physical mesh Th and the computational mesh Tc can be transformed into each other if and only if it satisfies (the following notations and concepts about the mesh are the same as those in [36])(4)Kdet⁡MK1/2=KcσhΩc,∀K∈Th,(5)12trFK′-1MK-1FK′-T=det⁡FK′-1MK-1FK′-T1/2,∀K∈Th,where K is an element of the physical mesh Th, Kc is an element of the computational mesh Tc corresponding to K, and |K| and |Kc|, respectively, show the volume of the corresponding element. Meanwhile, Ωc=∑Kc∈Tc|Kc|, det⁡(·) denotes the determinant of a matrix, and tr(·) stands for the trace of a matrix. σh=∑K∈ThKdet⁡(MK)1/2 represents the total error of the numerical solution, FK′ is the Jacobian matrix of the affine mapping FK: KC→K, and MK is the average of the metric tensor over the elements K. In this paper, we take M=det⁡αh+H((uh)n)-1/6αhI+H((uh)n), where I is an identity matrix and H(·) stands for a Hessian matrix. Equations (4) and (5) also satisfy the equidistribution principle and alignment principle on moving mesh, which means that, under the metric tensor M, each element has the same numerical error and is similar to a reference element K^; the reference element K^ is uniform. From (4) and (5), we find that the metric tensor M plays an important role in the adaptive moving mesh and controls the adaptation.If we want to get the mesh which satisfies the two principles, we need to minimize the energy function that contains the metric tensorM:(6)IhTh,Tc=θ∑KKdet⁡MK1/2trFK′-1MK-1FK′-Tdp/2+1-2θddp/2∑KKdet⁡MK1/2det⁡FK′-1det⁡MK-1/2p, where θ and p are dimensionless parameters. When we consider the equation to be just one-dimensional, (6) will have a relationship with (4) and θ=0; when the equation is two-dimensional, parameters θ,p are usually defined as 1/3 and 2, and d stands for the dimension in space. The energy function Ih is also a function about the physical mesh coordinate Xj=(xj,yj) and the computational mesh coordinate ξj=(ξj,ηj)(j=0,…,J), where J+1 stands for the number of mesh vertices. The coordinate ξ=ξ(X,t) is an inverse transformation of X=X(ξ,t); when j=0 or j=J, the coordinate ξj is the same as Xj. The energy function is too long to write; we write it in a simple form:(7)Ih=∑KKGFK′-1,det⁡FK′-1,MXK,XK.Assume that the current physical mesh Thn and (uh)n are known; by minimizing the energy function, we can get the corresponding computational mesh of the next time level, which can be expressed by (8)Tcn+1=arg⁡minTc⁡IhThn;Tc.In order to solve the minimum of energy function in (6), we need to take the derivative of the coordinate with respect to the energy function; that is,(9)∂Ih∂Xj=∂Ih∂ξj,j=0,…,J.It is not difficult to find that the equations are nonlinear and difficult to solve directly; in order to solve this problem, many scholars put forward a large number of methods to improve it. In this paper, we use the MMPDE [42] method to get the minimum value of the energy function (6). For convenience and efficiency, ξ-formulation MMPDE is made use of rather than x-formulation MMPDE; during the calculation, ξ-formulation only involves the differentiation of G with respect to its first two arguments in the energy function, and there is no need for the metric tensor M to update. The minimization process can be replaced as the derivative of the computational mesh coordinate with respect to t, which is an ordinary differential equation and has the following form:(10)dξjdt=Pjτ∂Ih∂ξjT,where the derivative ∂Ih/∂ξj is seen as a row vector, τ>0 is a relaxation factor which is used for adjusting the time scale on moving mesh, and P={P0,…,PJ} is a positive function used to make (10) have some desired properties above; in this paper, τ=10-4, Pj=-det⁡(M(Xj))1/2.At the current physical meshThn(n=0,…,K) and computational solution (uh)n, we define a reference computational mesh T^c=ξ^0,…,ξ^J, which has uniform elements; the mesh T^c also has the same relationship, grid connection, and number of mesh vertices as the physical mesh Th and computational mesh Tc. The three kinds of meshes just differ in the location of the mesh coordinate. During the mesh generation, by integrating (10) from tn to tn+1, we get the new computational mesh Tcn+1={ξjn+1}; there is no need for the metric tensor M to update. The relationship between the physical mesh and the computational mesh can be expressed by an affine mapping: Thn=Φh(Tcn+1); the physical mesh in the next time level can take the reference computational mesh by using a form of linear interpolation to get the approximate result: Thn+1=Φh(T^c).


### 2.2. Linear Finite Element Discretization of the 1D BBM-Burgers Equation with a Dissipative Term

In this subsection, we mainly introduce the linear finite element discretization of the 1D BBM-Burgers equation with a dissipative term. When(uh)n, xn, and the discretization of the equation are known, we will get the value (uh)n+1 at the new time level. The 1D BBM-Burgers equation with a fourth-order dissipative term has the formulation(11)ut+αux+βuux-γuxx-μuxxt+δuxxxx=F,x∈Ω,t∈0,T.Before solving (11) by the finite element method based on moving mesh, we have to introduce a new variable w=uxx rather than solving the equation directly; because a mixed derivative of the finite element approximation based on the moving mesh method is unable to be defined, then (11) becomes a coupled system. Before the discretization of the 2D equation with a fourth-order dissipative term, we need to introduce two new different variables to replace uxx and uyy; the equation then becomes three coupled equations; the other steps in the 2D case are similar to the 1D case. The new formulation of (11) is as follows: (12)ut-μwt+αux+βuux-γw+δwxx=F,x∈Ω,t∈0,T,w=uxx,x∈Ω,t∈0,T,ux,t=0,x∈∂Ω,t∈0,T.We divide the domain Ω into N uniform elements and J+1 mesh vertices; the initial physical mesh is Th0={xj0}(j=0,1,…,J). We define the linear basis function of the finite element method as φ. On a fixed mesh, the basis function φ is just related to x; on a moving mesh, it is also related to x and t. The basis function spans the space V0h={φ1(x,t),…,φJ-1(x,t)}, ∀φ∈V0h. We need to find whether uh(x,t)∈V0h and wh(x,t)∈V0h with the Dirichlet boundary condition equal to zero can satisfy the following system:(13)∫Ωut-μwt+αux+βuux-γw-Fφdx-∫Ωδwxφxdx=0,∀φ∈V0hΩ,t∈0,T,∫Ωwφ+μuxφxdx=0,∀φ∈V0hΩ,t∈0,T.Take uh(x,t)=∑j=1J-1uj(t)φj(x,t), wh(x,t)=∑j=1J-1wj(t)φj(x,t), where the derivative of the basis function [43] φj with respect to t is ∂φj/∂t=-∂φj/∂xX˙(x,t), and X˙ has the following formulation: X˙(x,t)=∑j=1J-1x˙jφj(x,t) (in the two-dimensional equation, X˙(x,y,t)=∑j=1J-1x˙jy˙jφjx,y,t), so (wh)t can be expressed as(14)wht=∑j=1J-1dwjdtφj-∑j=1J-1wj∂φj∂xX˙=∑j=1J-1dwjdtφj-∂wh∂xX˙.Take the basis function φ=φk(k=1,2,…,J-1) and substitute (14) into (13); we can get(15)∑j=1J-1∫Ωdujdt-μdwjdtφjφk-δwj∂φj∂x∂φk∂xdx+∫Ωα∂uh∂x+βuh∂uh∂x-γwh-F-∂uh∂x∑j=1J-1x˙jφj-μ∂wh∂x∑j=1J-1x˙jφjφkdx=0,∑j=1J-1∫Ωwjφjφk+uj∂φj∂x∂φk∂xdx=0.We rewrite (15) into a matrix form:(16)Mddtu-μw+f-δAw=0,Mw+Au=0,where M and A stand for Mkj=∫Ωφjφkdx and Akj=∫Ω(∂φj/∂x)(∂φk/∂x)dx and the vector f can be expressed as fk = ∫Ω(α∂uh/∂x+βuh∂uh/∂x-γwh-F-((∂uh/∂x)∑j=1J-1x˙jφj-μ(∂wh/∂x)∑j=1J-1x˙jφj))φkdx. On a moving mesh, with the boundary condition, the solution of (11) can be obtained by solving an ordinary differential equation firstly:(17)MddtA-1Mw+μMdwdt-f-A-1Mw,w,F+δAw=0.Only if the ODE (17) is solved can the solution u of the BBM-Burgers equation with a dissipative term be expressed, where u=-A-1Mw.All the process of discretization is based on the homogeneous boundary conditionu0=g; when the boundary condition is nonhomogeneous, we just need to subtract the integration at the boundary.


### 2.3. The Implementation of the Method Combined with Finite Element and Moving Mesh

In this subsection, we mainly introduce the implementation of the two-dimensional method. The ways to solve the equation usually can be divided into two kinds: a simultaneous solution procedure and an alternate solution procedure. The simultaneous solution procedure combines the finite element discretization with the mesh equation (10) into a large system and solves them together; the alternate solution procedure solves one by one. In this paper, we introduce a formulation of the second method which is called MP (the method solves the moving equation firstly and solves the physical equation later), which is more flexible and efficient. The method contains four important parts: initialize the mesh, move the mesh, discretize the equation, and integrate time from t0 to tK.Step 1 (initialize the mesh).
In the whole computational domain, we first divide the domain intoN uniform elements; the corresponding initial physical mesh is Th0, and the initial coordinate, initial value of the BBM- Burgers equation with a dissipative term, initial time step, and initial time are known; they are Xj0=(xj0,yj0), (uh)0, dt0, and t0.Step 2 (move the mesh).
(a) Suppose the current timet=tn(n=0,1,…,K), the numerical solution (uh)n of the equation, the coordinate Xjn, Thn of the mesh, and the time step dtn are known; (b) we make use of the energy function (9) and ODE function (10) to obtain the coordinate ξjn+1=(ξjn+1,ηjn+1), and then we derive the relationship Φh(·) between the new computational mesh Tcn+1={ξjn+1} and the physical mesh Thn, where Thn can be expressed by Φh(Tcn+1). When the reference computational mesh T^c and Φh(·) are known, we get the new physical mesh Thn+1=Φh(T^c) approximately; the new coordinate is X~jn+1; (c) when t∈[tn,tn+1], using linear interpolation Xj(t) = ((tn+1-t)/dtn)Xjn+(t-tn/dtn)X~jn+1, the velocity of the corresponding grid node is X˙j(t) = (X~jn+1-Xjn)/dtn. The fifth-order Radau IIA method is used in time discretization; if you are interested in more details, please refer to [44]. From the method, we can determine the time step dtn~; if the real time step dtn~<dtn, the next time level, the mesh coordinate, and the time step of next time level update: tn+1←tn+dtn~, Xjn+1←Xjn+dtn~X˙j, and dtn+1←dtn~.Step 3.
On the new physical meshThn+1, use the discrete version (16) based on the linear finite element method to solve (12) from time tn to tn+1; we get the new numerical solution at the new time level; the solution is (uh)n+1.Step 4.
Whentn+1<T, n←n+1, go to Step 2; else, stop the computation and save the value of solution.


## 2.1. Adaptive Moving Mesh Method

In this subsection, we mainly want to show the generation of physical meshThn+1 at a new time level by an adaptive moving mesh strategy based on the Moving Mesh Partial Differential Equation (MMPDE) [42]; at the current time level, the physical mesh Thn=Xjn and the numerical solution (uh)n are known (the mesh and value will be introduced clearly in the following parts). During the strategy, one of the crucially important rules is to view the nonuniform mesh (i.e., the physical mesh in this paper) as a uniform mesh under the metric tensor M; the tensor M is a symmetric and positive definite matrix. Under the metric tensor M, the physical mesh Th and the computational mesh Tc can be transformed into each other if and only if it satisfies (the following notations and concepts about the mesh are the same as those in [36])(4)Kdet⁡MK1/2=KcσhΩc,∀K∈Th,(5)12trFK′-1MK-1FK′-T=det⁡FK′-1MK-1FK′-T1/2,∀K∈Th,where K is an element of the physical mesh Th, Kc is an element of the computational mesh Tc corresponding to K, and |K| and |Kc|, respectively, show the volume of the corresponding element. Meanwhile, Ωc=∑Kc∈Tc|Kc|, det⁡(·) denotes the determinant of a matrix, and tr(·) stands for the trace of a matrix. σh=∑K∈ThKdet⁡(MK)1/2 represents the total error of the numerical solution, FK′ is the Jacobian matrix of the affine mapping FK: KC→K, and MK is the average of the metric tensor over the elements K. In this paper, we take M=det⁡αh+H((uh)n)-1/6αhI+H((uh)n), where I is an identity matrix and H(·) stands for a Hessian matrix. Equations (4) and (5) also satisfy the equidistribution principle and alignment principle on moving mesh, which means that, under the metric tensor M, each element has the same numerical error and is similar to a reference element K^; the reference element K^ is uniform. From (4) and (5), we find that the metric tensor M plays an important role in the adaptive moving mesh and controls the adaptation.If we want to get the mesh which satisfies the two principles, we need to minimize the energy function that contains the metric tensorM:(6)IhTh,Tc=θ∑KKdet⁡MK1/2trFK′-1MK-1FK′-Tdp/2+1-2θddp/2∑KKdet⁡MK1/2det⁡FK′-1det⁡MK-1/2p, where θ and p are dimensionless parameters. When we consider the equation to be just one-dimensional, (6) will have a relationship with (4) and θ=0; when the equation is two-dimensional, parameters θ,p are usually defined as 1/3 and 2, and d stands for the dimension in space. The energy function Ih is also a function about the physical mesh coordinate Xj=(xj,yj) and the computational mesh coordinate ξj=(ξj,ηj)(j=0,…,J), where J+1 stands for the number of mesh vertices. The coordinate ξ=ξ(X,t) is an inverse transformation of X=X(ξ,t); when j=0 or j=J, the coordinate ξj is the same as Xj. The energy function is too long to write; we write it in a simple form:(7)Ih=∑KKGFK′-1,det⁡FK′-1,MXK,XK.Assume that the current physical mesh Thn and (uh)n are known; by minimizing the energy function, we can get the corresponding computational mesh of the next time level, which can be expressed by (8)Tcn+1=arg⁡minTc⁡IhThn;Tc.In order to solve the minimum of energy function in (6), we need to take the derivative of the coordinate with respect to the energy function; that is,(9)∂Ih∂Xj=∂Ih∂ξj,j=0,…,J.It is not difficult to find that the equations are nonlinear and difficult to solve directly; in order to solve this problem, many scholars put forward a large number of methods to improve it. In this paper, we use the MMPDE [42] method to get the minimum value of the energy function (6). For convenience and efficiency, ξ-formulation MMPDE is made use of rather than x-formulation MMPDE; during the calculation, ξ-formulation only involves the differentiation of G with respect to its first two arguments in the energy function, and there is no need for the metric tensor M to update. The minimization process can be replaced as the derivative of the computational mesh coordinate with respect to t, which is an ordinary differential equation and has the following form:(10)dξjdt=Pjτ∂Ih∂ξjT,where the derivative ∂Ih/∂ξj is seen as a row vector, τ>0 is a relaxation factor which is used for adjusting the time scale on moving mesh, and P={P0,…,PJ} is a positive function used to make (10) have some desired properties above; in this paper, τ=10-4, Pj=-det⁡(M(Xj))1/2.At the current physical meshThn(n=0,…,K) and computational solution (uh)n, we define a reference computational mesh T^c=ξ^0,…,ξ^J, which has uniform elements; the mesh T^c also has the same relationship, grid connection, and number of mesh vertices as the physical mesh Th and computational mesh Tc. The three kinds of meshes just differ in the location of the mesh coordinate. During the mesh generation, by integrating (10) from tn to tn+1, we get the new computational mesh Tcn+1={ξjn+1}; there is no need for the metric tensor M to update. The relationship between the physical mesh and the computational mesh can be expressed by an affine mapping: Thn=Φh(Tcn+1); the physical mesh in the next time level can take the reference computational mesh by using a form of linear interpolation to get the approximate result: Thn+1=Φh(T^c).


## 2.2. Linear Finite Element Discretization of the 1D BBM-Burgers Equation with a Dissipative Term

In this subsection, we mainly introduce the linear finite element discretization of the 1D BBM-Burgers equation with a dissipative term. When(uh)n, xn, and the discretization of the equation are known, we will get the value (uh)n+1 at the new time level. The 1D BBM-Burgers equation with a fourth-order dissipative term has the formulation(11)ut+αux+βuux-γuxx-μuxxt+δuxxxx=F,x∈Ω,t∈0,T.Before solving (11) by the finite element method based on moving mesh, we have to introduce a new variable w=uxx rather than solving the equation directly; because a mixed derivative of the finite element approximation based on the moving mesh method is unable to be defined, then (11) becomes a coupled system. Before the discretization of the 2D equation with a fourth-order dissipative term, we need to introduce two new different variables to replace uxx and uyy; the equation then becomes three coupled equations; the other steps in the 2D case are similar to the 1D case. The new formulation of (11) is as follows: (12)ut-μwt+αux+βuux-γw+δwxx=F,x∈Ω,t∈0,T,w=uxx,x∈Ω,t∈0,T,ux,t=0,x∈∂Ω,t∈0,T.We divide the domain Ω into N uniform elements and J+1 mesh vertices; the initial physical mesh is Th0={xj0}(j=0,1,…,J). We define the linear basis function of the finite element method as φ. On a fixed mesh, the basis function φ is just related to x; on a moving mesh, it is also related to x and t. The basis function spans the space V0h={φ1(x,t),…,φJ-1(x,t)}, ∀φ∈V0h. We need to find whether uh(x,t)∈V0h and wh(x,t)∈V0h with the Dirichlet boundary condition equal to zero can satisfy the following system:(13)∫Ωut-μwt+αux+βuux-γw-Fφdx-∫Ωδwxφxdx=0,∀φ∈V0hΩ,t∈0,T,∫Ωwφ+μuxφxdx=0,∀φ∈V0hΩ,t∈0,T.Take uh(x,t)=∑j=1J-1uj(t)φj(x,t), wh(x,t)=∑j=1J-1wj(t)φj(x,t), where the derivative of the basis function [43] φj with respect to t is ∂φj/∂t=-∂φj/∂xX˙(x,t), and X˙ has the following formulation: X˙(x,t)=∑j=1J-1x˙jφj(x,t) (in the two-dimensional equation, X˙(x,y,t)=∑j=1J-1x˙jy˙jφjx,y,t), so (wh)t can be expressed as(14)wht=∑j=1J-1dwjdtφj-∑j=1J-1wj∂φj∂xX˙=∑j=1J-1dwjdtφj-∂wh∂xX˙.Take the basis function φ=φk(k=1,2,…,J-1) and substitute (14) into (13); we can get(15)∑j=1J-1∫Ωdujdt-μdwjdtφjφk-δwj∂φj∂x∂φk∂xdx+∫Ωα∂uh∂x+βuh∂uh∂x-γwh-F-∂uh∂x∑j=1J-1x˙jφj-μ∂wh∂x∑j=1J-1x˙jφjφkdx=0,∑j=1J-1∫Ωwjφjφk+uj∂φj∂x∂φk∂xdx=0.We rewrite (15) into a matrix form:(16)Mddtu-μw+f-δAw=0,Mw+Au=0,where M and A stand for Mkj=∫Ωφjφkdx and Akj=∫Ω(∂φj/∂x)(∂φk/∂x)dx and the vector f can be expressed as fk = ∫Ω(α∂uh/∂x+βuh∂uh/∂x-γwh-F-((∂uh/∂x)∑j=1J-1x˙jφj-μ(∂wh/∂x)∑j=1J-1x˙jφj))φkdx. On a moving mesh, with the boundary condition, the solution of (11) can be obtained by solving an ordinary differential equation firstly:(17)MddtA-1Mw+μMdwdt-f-A-1Mw,w,F+δAw=0.Only if the ODE (17) is solved can the solution u of the BBM-Burgers equation with a dissipative term be expressed, where u=-A-1Mw.All the process of discretization is based on the homogeneous boundary conditionu0=g; when the boundary condition is nonhomogeneous, we just need to subtract the integration at the boundary.


## 2.3. The Implementation of the Method Combined with Finite Element and Moving Mesh

In this subsection, we mainly introduce the implementation of the two-dimensional method. The ways to solve the equation usually can be divided into two kinds: a simultaneous solution procedure and an alternate solution procedure. The simultaneous solution procedure combines the finite element discretization with the mesh equation (10) into a large system and solves them together; the alternate solution procedure solves one by one. In this paper, we introduce a formulation of the second method which is called MP (the method solves the moving equation firstly and solves the physical equation later), which is more flexible and efficient. The method contains four important parts: initialize the mesh, move the mesh, discretize the equation, and integrate time from t0 to tK.Step 1 (initialize the mesh).
In the whole computational domain, we first divide the domain intoN uniform elements; the corresponding initial physical mesh is Th0, and the initial coordinate, initial value of the BBM- Burgers equation with a dissipative term, initial time step, and initial time are known; they are Xj0=(xj0,yj0), (uh)0, dt0, and t0.Step 2 (move the mesh).
(a) Suppose the current timet=tn(n=0,1,…,K), the numerical solution (uh)n of the equation, the coordinate Xjn, Thn of the mesh, and the time step dtn are known; (b) we make use of the energy function (9) and ODE function (10) to obtain the coordinate ξjn+1=(ξjn+1,ηjn+1), and then we derive the relationship Φh(·) between the new computational mesh Tcn+1={ξjn+1} and the physical mesh Thn, where Thn can be expressed by Φh(Tcn+1). When the reference computational mesh T^c and Φh(·) are known, we get the new physical mesh Thn+1=Φh(T^c) approximately; the new coordinate is X~jn+1; (c) when t∈[tn,tn+1], using linear interpolation Xj(t) = ((tn+1-t)/dtn)Xjn+(t-tn/dtn)X~jn+1, the velocity of the corresponding grid node is X˙j(t) = (X~jn+1-Xjn)/dtn. The fifth-order Radau IIA method is used in time discretization; if you are interested in more details, please refer to [44]. From the method, we can determine the time step dtn~; if the real time step dtn~<dtn, the next time level, the mesh coordinate, and the time step of next time level update: tn+1←tn+dtn~, Xjn+1←Xjn+dtn~X˙j, and dtn+1←dtn~.Step 3.
On the new physical meshThn+1, use the discrete version (16) based on the linear finite element method to solve (12) from time tn to tn+1; we get the new numerical solution at the new time level; the solution is (uh)n+1.Step 4.
Whentn+1<T, n←n+1, go to Step 2; else, stop the computation and save the value of solution.


## 3. Numerical Examples

Example 1 (a single solitary wave in 1D BBM-Burgers equation).
In the one-dimensional model,(18)ut+αux+βuux-γuxx-μuxxt=F,x∈-20,30,t∈0,T,u-20,t=u30,t=0,t∈0,T, where the coefficients α=γ=μ=1, β=-1, the final time T=10, the exact solution of the equation is uex(x,t)=sech⁡(x-t), and the source term F=F(x,t)=(1-6tanh3⁡x-t-2tanh2⁡x-t+sech⁡x-t+5tanh⁡(x-t))sech⁡(x-t).We solve the numerical solution with the method we introduced; with the different element numbers, the corresponding error and the order of convergence for both moving mesh and fixed mesh are reported in Table1. When the element number is the same, the error based on moving mesh is always smaller than that on fixed mesh, and the convergence order of the two kinds of meshes can lead to the second order in the spatial direction. Comparison between the exact solution and the numerical solution based on moving mesh at the final time T=10 is shown in Figure 1. From the figure, we can see that the grid points mainly focus on the domain, where the solution fluctuates greatly, and match the exact solution well. With time going on, the moving mesh trajectory is plotted in Figure 1; the feature of the moving process indicates that most of the points focus on nearby solitary waves and move with the propagation of the solitary wave, which also demonstrates the superiority of the moving mesh method.Table 1
Example1. The L2 and L∞ error and convergence order on moving and fixed meshes for the single solitary wave.
N
Moving mesh
Fixed mesh
L
2 error
Order
L
∞ error
Order
L
2 error
Order
L
∞ error
Order
20
3.65
E
-
2
1.70
E
-
1
8.35
E
-
1
3.20
E
-
0
40
9.12
E
-
2
2.00
3.99
E
-
2
2.09
1.54
E
-
1
2.44
1.12
E
-
0
1.51
80
2.49
E
-
3
1.87
1.05
E
-
2
1.93
4.09
E
-
2
1.91
3.74
E
-
1
1.58
160
6.53
E
-
4
1.93
2.74
E
-
3
1.93
1.04
E
-
2
1.98
1.00
E
-
1
1.90
320
1.84
E
-
4
1.83
7.65
E
-
4
1.84
2.60
E
-
3
1.99
2.55
E
-
2
1.97
640
5.01
E
-
4
1.88
2.06
E
-
4
1.89
6.51
E
-
4
2.00
6.41
E
-
3
1.99Figure 1
Example1. (a) The exact solution and the numerical solution based on adaptive moving mesh with element number N=200 at t=10; (b) the corresponding mesh trajectories.
(a)
Computed solution and exact solution
(b)
Mesh trajectoriesExample 2 (the elastic collision of three solitary waves in 1D BBM-Burgers equation).
In this example, we study the elastic collision of three solitary waves and simulate numerical solutions based on moving mesh and fixed mesh; the BBM-Burgers equation has the following form:(19)ut+ux+uux-0.01uxx-uxxt=0,x∈-50,300,t∈0,50,u-50,t=u300,t=0,t∈0,50,ux,0=u0x,x∈-50,300, with the initial condition u0(x)=∑j=133djsech2⁡(kj(x-xj))2, where k1=0.39, k2=0.3, k3=0.3, x1=10, x2=28, x3=52, dj=4kj2/(1-4kj2), and the final time T=50.
In this example, the exact solution of the equation is unknown; we mainly want to show the superiority of the moving mesh. In Figure2, the conservation laws (2) are plotted; we can see that ΔE1(T) is close to zero and ΔE2(T) is a little bigger than zero. The two conservation laws are not conserved by finite element approximation on the moving mesh, because, in this case, X˙≡0 and both matrices A and M are time-dependent (if you want to know more, read [45]). Since the exact solution is unknown, we take the element number N=8000 on the fixed mesh as a reference solution and plot the numerical solutions of moving mesh and fixed mesh with element number N=400 at final time T=50 in Figure 3. It can be seen that the numerical solution of fixed mesh is on the left of the reference solution and the wave peaks are shorter than those in the reference solution; there are several small waves in the left of the first wave on the fixed mesh. The result also shows that the numerical solution of moving mesh can match the reference better than that based on the fixed mesh. If we consider the computer memory and computational efficiency, the moving mesh is better for that it moves its mesh points quickly to the place where the numerical solution changes rapidly to get higher distinguishability with fewer elements. The propagation process of the three solitary waves is plotted in Figure 4.Figure 2
Example2. The variation of mass function (E1) and energy function (E2) for the elastic collision of three solitary waves in the BBM-Burgers equation.Figure 3
Example2. The comparison between the numerical solution and the reference solution on moving and fixed meshes for the interaction of three solitary waves with N=400.Figure 4
Example2. The propagation process for the elastic collision of three solitary waves at different times with N=400.Example 3 (different coefficients in 1D BBM-Burgers equation with a high-order dissipative term).
In this case, we mainly consider the following BBM-Burgers equation with a fourth-order dissipative term; the equation is(20)ut+ux+2uux-uxxt-γuxx+δuxxxx=0,x∈-100,150,t∈0,20,u-100,t=u150,t=0,t∈0,20,ux,0=u0x,x∈-100,150. The initial condition u0(x)=0.15sech⁡kx-402, where k=0.5×1/11; when the coefficients γ,δ change, the maximum value Umax of the numerical solution changes as well.
In this example, we mainly want to talk about the dissipation effect of the dissipative termsuxx and uxxxx in the BBM equation. Numerical solutions obtained with different coefficients and one of the coefficients fixed at T=20 are shown in Figure 5. When two BBM equations with one dissipative term uxx or uxxxx have the same dissipation effect or similar wave peak at the final time, the corresponding coefficient of dissipative term is plotted in Figure 6. From these figures, we find that when the coefficient γ of the second-order dissipative term increases and the coefficient δ of the fourth-order dissipative term is fixed, the numerical solution decreases a lot; however, when the coefficient δ increases and the coefficient γ is fixed, the numerical solution changes a little. The numerical simulation provides three sets of data to show the result in Figure 5: δ=1; when γ=0.1,1,10, the corresponding maximum values of the numerical solutions are Umax=0.139,0.1324,0.1091, γ=1; when δ=0.1,1,10, the corresponding values are Umax=0.0947,0.0936,0.0866. A lot of numerical tests have been carried out to adjust and confirm the corresponding coefficient of different dissipative terms to keep the equation with dissipation effect or similar wave peak at T=20. The value of the wave peak is Umax=0.0764 at T=20, and the corresponding data are shown in Figure 6. When δ=0, γ=1.92, with time going on, Umax drops steadily; however, when γ=0, δ=100, the value drops quickly before time t=2, and then the value drops slowly and steadily. We simulate a series of numerical examples with different coefficients; all of these data show us that the dissipative term uxx plays a more important role in dissipation than the dissipative term uxxxx.Figure 5
Example3. The numerical solutions at the final time for different coefficients. (a) The numerical solutions of the equation with δ=1, γ=0.1,1,10; (b) the numerical solutions of the equation with γ=1, δ=0.1,1,10.
(a)
δ
=
1, γ changes
(b)
γ
=
1, δ changesFigure 6
Example3. The corresponding dissipative term coefficient when Umax has the same value at the beginning and final time; (a) Umax of the equation with just one dissipative uxx and the coefficient γ=1.92; (b) Umax of the equation with just one dissipative uxxxx and the coefficient δ=100.
(a)
γ
=
1.92, δ=0
(b)
γ
=
0, δ=100Example 4 (two solitary waves in 2D BBM-Burgers equation).
This example has two solitary waves and we want to verify the convergence order for a 2D BBM-Burgers equation; the equation has the following formulation:(21)ut+α+βcos⁡uux+uy-γuxx+uyy-μuxxt+uyyt=F,x,y∈Ω,t∈0,T, where α,γ,μ=1, β=-1, and the exact solution of the equation is uex(x,y,t)=et(sech2⁡x1+sech2⁡(x2)), where x1=x+y-1, x2=x+y+1, Ω=(-2,2)×(-2,2), T=1. To maintain the balance of the equation, F = F(x,y,t) = et((9(sech2⁡x1+sech2⁡x2) − 4(sech2⁡x1tanh⁡x1+sech2⁡x2tanh⁡x2) − 24(sech2⁡x1tanh2⁡x1+sech2⁡x2tanh2⁡x2) + 4cos⁡(et(sech2⁡x1+sech2⁡x2)))(sech2⁡x1tanh⁡x1+sech2⁡x2tanh⁡(x2))).Numerical solutions are shown in Table2 and Figure 7. From Table 2, we can find that the finite element method leads to the same second order of convergence for both fixed and moving meshes, and the convergence result is similar to Example 1. The error of the numerical solution on the moving mesh is smaller than that on the fixed mesh with the same element number. What is more, with fewer mesh elements, the moving mesh reaches second-order convergence more quickly than the fixed mesh. The trajectory and the numerical solution of both the moving mesh and the fixed mesh at final time can be seen in Figure 7. From the figure, we can see that most points focus on the place where the numerical solution changes a lot.Table 2
Example4. The numerical error and convergence order for the 2D BBM-Burgers equation.
N
Moving mesh
Fixed mesh
L
2 error
Order
L
∞ error
Order
L
2 error
Order
L
∞ error
Order
100
3.08
E
-
2
1.41
E
-
1
4.10
E
-
2
2.28
E
-
1
400
8.04
E
-
3
1.94
4.44
E
-
2
1.67
1.20
E
-
2
1.77
7.52
E
-
2
1.60
1600
2.09
E
-
3
1.94
1.12
E
-
2
1.99
3.11
E
-
3
1.95
2.02
E
-
2
1.90
6400
5.48
E
-
4
1.93
2.87
E
-
3
1.96
7.83
E
-
4
1.99
5.19
E
-
3
1.96
25600
1.42
E
-
4
1.95
7.41
E
-
4
1.95
1.96
E
-
4
2.00
1.30
E
-
3
1.99Figure 7
Example4. Meshes (a, c) and the corresponding numerical solutions (b, d) with N=6400.
(a)
Moving mesh
(b)
The corresponding numerical solution
(c)
Fixed mesh
(d)
The corresponding numerical solution


## 4. Conclusion

In this paper, we study the numerical model based on adaptive moving mesh finite element method to simulate 1D and 2D BBM-Burgers equations with a high-order dissipative term. The method can be used in equations with mixed derivatives and high-order derivative terms; on the moving mesh, new variables need to be introduced to transform the equation into coupled equations. A series of numerical examples demonstrate that the finite element method based on moving mesh keeps a second-order convergence in space as that based on fixed mesh. With the same element number, the finite element moving mesh method improves the distinguishability and provides a smaller error than the fixed mesh. Meanwhile, numerical tests indicate that the dissipative termuxx plays a more important role than uxxxx in dissipation. We compare two BBM equations with different dissipative terms uxx and uxxxx; to keep the maximum value Umax of the numerical solutions the same at the final time, the coefficient of uxxxx is around 50 times that of uxx in the same BBM equation with different dissipative terms in the dissipation aspect.



---

*Source: 3427376-2017-11-29.xml*