Burgers-Korteweg-de Vries复合方程的格子Boltzmann方法模拟

段雅丽, 陈先进, 孔令华（.中国科学技术大学数学科学学院,安徽合肥 3006；.江西师范大学数学与信息科学学院,江西南昌 3300）

Burgers-Korteweg-de Vries复合方程的格子Boltzmann方法模拟

段雅丽1, 陈先进1, 孔令华2
（1.中国科学技术大学数学科学学院,安徽合肥 230026；2.江西师范大学数学与信息科学学院,江西南昌 330022）

针对Burgers-Korteweg-de Vries（cBKdV）复合方程提出一种格子Boltzmann模型.通过恰当地处理色散项uxxx并运用Chapman-Enskog展开从格子Boltzmann方程推导出宏观方程,从而得到联系微观量与宏观量的局部平衡分布函数.对不同微分方程进行数值实验,数值解与解析解非常吻合,相比于其它数值结果,该格子Boltzmann模型的数值结果更精确,说明该数值模型的高效性.

格子Boltzmann模型；Burgers-KdV复合方程；Chapman-Enskog展开

0 Introduction

In the simulation of radiation hydrodynamic Lagrange problems,diffusion problem often appears as a crucial subproblem.It is important to study discrete schemes with high accuracy for these problems,especially for those with large deformation.

In recent years,to avoid difficulties in mesh generation for complex problems,meshless methods have received intensive attention,and some of them have been applied to solve diffusion problems.Shen et al[1].solved the Poisson equation as an example of application to awell-designed algorithm for selecting neighboring points.Reference[2]constructed minimal positive stencils in meshfree finite difference methods for the Poisson equation by a linear minimization approach.Reference[3]presented adaptive meshless discretization for the Poisson equation based on radial basisfunctionmethods（RBFs）.Reference[4]extended a numerically stable RBF-QR formulation of RBF approximation in the aspects of computing differentiation matrices and stencil weights,in which the Poisson equation was solved as an example of numerical implement.Reference[5]gave a generalized finite difference discretization to the Poisson equation with piecewise constant diffusion conductivity,which was based on an MLSapproach.Thismethod used a constrain condition for flux rather than governing equation at interface to derive a discrete scheme.More recently,Ref.[6] proposed a space-time diffuse approximation method,in which a weight function was introduced to remove some spurious oscillations in computation for problems with high temporally discontinuous heat sources.

In the present paper,we consider diffusion problems with nonlinear and/or discontinuous conductivity by a finite directional difference method（FDDM）.The FDDM is a meshlessmethoddefined on scattered point distributions,which is first proposed in Ref.[1]as a finite pointmethod （FPM）based on directional differences,and renamed as FDDM to distinguish from FPM proposed in Ref.[7].The FPM[7]is based on a weighted least square interpolation of point data and point collocation for evaluating approximation integrals,while the FDDM can be viewed as a generalization of classical finite differencemethod defined on uniform point distribution,which ismore difficult to perform due to disorders of scattered point distributions.In Refs.[1,8],explicit numerical formulae for approximations to directional differentialswere derived with expected accuracy by using information of proper scattered points.Above all,by virtue of explicit expressions,solvability conditions of numerical derivativeswere rigorously given,which gave a general guiding principle for selecting neighboring points avoiding singularity in computing derivatives.

In the FDDM regime,the present paper develops an approach to solving diffusion problemswith nonlinear conductivities,which has advantages as follows：It leads tominimal stencils,coefficients of the resulted scheme are given explicitly avoiding solving matrix equations,and well-designed method for selecting neighboring points guarantees that the issue of singularity never emerges.Reference[9]also investigated this problem.However,the present paper deals with the diffusion operator and the nonlinear term more rigorously than Ref.[9].

Furthermore,when the diffusion conductivity is discontinuous,we propose a scheme for discretizingmultimedia interface condition by the FDDM method.To discretize flux on interface, Ref.[5]employsmore neighbors of the master point,while we employ only five neighbors of the master point on each side of the interface,resulting flux with second-order accuracy.This idea is also explored to discretize energy flux of diffusion equation on unstructured meshes（for details see Ref.[10]）.

The rest of this paper is arranged as follows：Section 1 presents some preliminaries.Section 2 formulates numerical differentials on scattered point distributions.Section 3 constructsmethodology for solving nonlinear diffusion problems and discontinuous problems,and gives several numerical validations.Finally,concluding remarks aremade in Section 4.

1 Prelim inaries

To simplify presentation,we first introduce denotations and definitions as defined in Ref.[1].Let us denote

·i the index of point（xi,yi）and“O”a specific point（x0,y0）；

·ljthe j-th direction vector from point O and ejthe corresponding unit vector；·Δlithe distance between point“O”and i；

·ui＝u（xi,yi）the function value of u（x,y）at point i；and

·Δui＝ui－u0the difference of function u（x,y）.

We also have

·〈ij〉∶＝〈ijO〉,i.e.,k is a special point as“O”in the expression〈i jk〉.

Here,there is a little changemade from Ref.[1]in that the denotation（·,·）is changed into〈·,·〉0, because（·,·）is just a special case of〈·,·〉.

In this paper,since a large amount of operations for indices denoting direction vectors are required,we introduce an operation for indices defined as in Ref.[1]：

Definition 1 （Algorithm○k）Given i,j,k（k≥3）positive integers,an addition of i and j with period of k is defined by

where s is a nonnegative integer satisfying inequality sk＜i＋j≤（s＋1）k.The operation○k can be also expressed by

where（i＋j－1）（mod k）represents the remainder of（i＋j－1）modulo k.

Remark 1.1 The operation○k is defined only for indices i,j,…,which is irrelevant to the numbering and ordering of the directions denoted by these indices（See Fig.1）.

Fig.1 Illustration for k directions

Since the FDDM is based on directional differences,we need relations between directional derivatives,which will significantly help numerical discretization in the FDDM regime.

Firstwe state relevant resultswith constant coefficients.

Lemma 1 （see Ref.[1]）Given three unit vectors e1,e2,e3at point（x,y）,for a smooth function u（x,y）on domainΩ⊂R2,it follows that

Lemma 2 （see Ref.[1]）Given four unit vectors ei（i＝1,2,3,4）at point（x,y）,for secondorder directional derivatives of a smooth function u（x,y）on the domainΩ⊂R2,it follows that

To handle practical problems,more general results are mandatory.By simple deduction,wederive relations between directional derivatives with variable coefficients as follows.

Theorem 1 Given four arbitrary unit vectors ei（i＝1,2,3,4）,for smooth functions u（x,y）and κ（x,y）on domainΩ⊂R2,it follows that

2 Numerical differentials and solvability

Before constructing a discrete scheme for PDEs,we first deal with the basic issue,i.e.,the approximations to directional derivatives.

Suppose that for a given point O and its five neighboring points（xi,yi）（i＝1,…,5）whose indices i＝1,…,5 are numbered freely（see Fig.2）,the differencesΔui（i＝1,…,5）of the smooth function u（x,y）are available.

Fig.2 Point O and its five neighbors

Let us denote by

We also have

and

where

With detailed analysis on the five-point formulae presented in Ref.[1],we formulatenumerical approximations for the first-order and second-order directional derivatives of the smooth function u（x,y）at point O in thematrix forms

where

E is a unitmatrix,and

We call

as the solvability condition.

According to the result of Ref.[1],we note that Eq.（8）is second-order accurate as to the approximation to the first-order derivatives,and Eq.（9） is first-order accurate as to the approximation to the second-order derivatives.

We also notice that whether solvability condition（11）is satisfied or not is a key issue.Reference[1]dealtwith this issue,and presented an algorithm of selecting neighboring points for solving diffusion problems,which will be employed in the present paper.

From the above formulation we can learn that for the second-order directional derivatives,only information of themaster point and its five neighbors were used.One knows that five neighboring points are of the least number for approximating the second-order directional derivatives as consider the consistency.This is important to construct schemeswithminimal stencils.

3 Discretization methodology

In this section,we restrict our attention to constructing the discretization methodology for numerically solving diffusion equations,which have the form

with initial condition

and boundary condition

whereΩis an open bounded domain with smooth boundaryƏΩ,T is the final time,κis the nonlinear diffusion coefficient and maybe discontinuous,and f,g1and g2are given functions.

Hereafter,we always impose discretization toΩandƏΩby scattered point distribution,and denote the resulted discrete point set byΩhandƏΩh,respectively,and＝Ωh∪ƏΩh.

It is obvious that discretization of the diffusion operator∇·（κ（x,y,u）∇u）is a key issue.The first step to employ FDDM is to express the diffusion operator by directional differentials.

3.1 Expression of∇·（κ∇u）by directional differentials

In this section and Section 3.2,we consider the case thatκis a smooth function.

Given e1,e2,e3,three nonparallel unit vectors from point O（see Fig.3）,and eI,eJ,unitvectors in x,y axis directions,respectively,then by means of Eq.（5）,∇·（κ∇u）can be expressed by

Fig.3 Three nonparallel unit vectors

which can be simplified into

3.2 Discretization of∇·（κ∇u）

Fig.4 Point O and its five neighbors

We suppose that every point inΩhhas five steady distribution neighbors satisfying solvability condition（11）.

Given a point O（x0,y0）and its five neighbors（xi,yi）（i＝1, …,5）（see Fig.4）,denote“i′”as themiddle pointof the segment,κ0＝κ（x0,y0,u0）,κi＝κ（xi,yi,ui）,＝（κ0＋κi）/2, andκi′＝κ（ui′）,i＝1,2,…,5.It is obvious that

Motivated by the technique in constructing numerical formula （9）for the second-order directional differentials,we have

To simplify presentation,we denote by

Moreover,noticing that“i′”is themiddle point of,togetherwith Eq.（16）we have

By Eq.（17）,it consequently follows that

Therefore,by Eq.（15）,we have

Remark 3.1 （1）Note that scheme（21）reduces to a classical finite difference on uniform point distribution.

（2）It is obvious that Eq.（21）yields stencils ofminimal size.

3.3 Scheme for discretizing multimedia interface condition

In this section,we construct discrete schemes for discontinuous diffusion problems.

For simplicity,we consider the case of two subdomains,Ω1andΩ2,separated by an interface Γ,and suppose thatκis discontinuous throughΓwith respect to spatial variables x and y,but is continuous with respect to u.Denote the inward and outward unit normal vectors ofΓby n－and n＋, and the associated conductivities byκ－andκ＋,respectively.

Aswe have derived the scheme proposed in Section 3.2 for the case thatκis continuous,now we need only to focus on constructing discrete scheme for the interface condition.We first place appropriate points onΓ.For any given point O belonging toΓand corresponding functional values u0,the procedure of ourmethod can be outlined as follows：

1）Choose five neighbors of O on each side ofΓ,respectively,and denote them by1,2,3,4, 5 corresponding to n＋side,1′,2′,3′,4′,5′corresponding to n－side（see Fig.5,the selected neighbors are marked bold,and other points are indicated by white circles）,and their function values ui（i＝1,…,5；1′,…,5′）.In this step,we choose the nearest five points satisfying the solvability condition on each side.

Fig.5 Illustration of selecting neighboring points for a point on interface

The detailed procedure is as follows.

By the relation between the first-order directional derivatives（3）,we have

Here,Δl＋is an auxiliary quantity,which will be eliminated in the following deduction,M＋is defined following Eq.（10）,and a1j,a2j（j＝1,2,…,5）are as given in Eq.（18）.Therefore,wehave

Likewise,we have

Here,Δl－is also an auxiliary quantity similar toΔl＋,and M－,a′1jand a′2j（j＝1,2,…,5）are similar to M＋,a1jand a2j（j＝1,2,…,5）.

LetΔl＋＝Δl－,and denote by

Applying Eqs.（24）and（25）to the continuous condition of energy flux Eq.（22）,and eliminating Δl＋andΔl－,we have

which results in the discrete scheme for the interface condition as follows：

To summarize,for a point belonging toΩ1（orΩ2）,select its five neighbors from those belonging toΩ1∪Γ（orΩ2∪Γ）,and then discretize the control equation（12）by the scheme（21）；For a point belonging toΓ,the above scheme constructed for the interface condition is employed.After the discretization,we derive a large sparse system of algebraic equations,which can be solved by various iterativemethods.

3.4 Numerical Results

This section presents several numerical examples with different computational domains and different point distributions to investigate the accuracy and efficiency of the proposed approach.

In the subsequent computation,the resulted nonlinear systems are solved by a classical Picard iterative process.And the linear systems are solved by a biconjugate gradient stabilized algorithm.

To investigate convergence results of the proposed method,we define discrete norm errors as

The convergence rate is

where h1and h2are average distances corresponding to N1and N2,respectively.

Exam p le 1 Consider the problem

Fig.6 Point distribution on computational domain

where f（x,y）＝2（ey－x－ex－y）,Ωis a semi-circle with two semi-circles cut out as shown in Fig.6.

Note that discretizingΩby a good-quality mesh is a complex work as consider the corners （see e.g.the left bottom ofΩ）,meanwhile distributing scattered points on it is rather easy.We employ this example to investigate the proposed method on an irregular computational domain. Corresponding numerical results are presented in Table 1.

From Table 1, we can see that the approximation to the solution is almost secondorder accurate,and to the first-order directional derivatives of the solution is higher than first-order accurate.This indicates that the proposed method workswell on irregular computational domain.

Table 1 Errors of u and its first-order derivatives

Exam ple 2 Solve an nonlinear boundary problem

whereκ（u）＝u,and f（x,y）,g（x,y）are given by the exact solution u（x,y）＝2＋cos（πx）＋sin（πy）.Ωis a unit square,which is discretized by three types of scattered points as shown in Fig.7.

In this example,the tolerance of nonlinear iterative is‖Us＋1－Us‖≤10－8,where Usand Us＋1represent the numerical solution of two neighboring iterative steps,respectively.

Corresponding numerical results are presented in Table 2.From this table,we see that approximations to the solution in the-norm and the-norm are almost second-order accuracy, except that the convergence rate in the-norm on Z-type point distribution is slightly low,which is also satisfactory as taking account of the highly anisotropic point distribution in this case.

Fig.7 Three types of point distributions in square domain,left：uniform；middle：random；right：Z-type

Table 2 Convergence results for nonlinear diffusion problem on three types of point distributions

Fig.8 Random point distribution

Exam ple 3 Solve a parabolic problem

where T＝10－3,and f（x,y,t）,g1（x,y）and g2（x,y,t）are given by the exact solution u（x,y,t）＝e－π2 t（2＋cos（πx）＋sin（πy））.Ωis a unit square,which is discretized as shown in Fig.8.

Objective of this example is to compare FDDM with the classical least squaremethod（LSQ）.The time step is chosen asΔt＝10－5,and utis discretized by a backward difference formula.In LSQ,neighboring points are the nearest ones,the numbers of which are selected as 10,20,and 40,respectively.The numerical results are graphically depicted in Fig.9.

One can observe the following：

· Both methods have almost the same convergence rate.Errors of the LSQ are higher than those of the FDDM.

·For LSQ,the accuracy does not increase with increasing numbers of neighboring points.

By this test problem,FDDM is also compared with LSQmethod in terms of computational cost.For a sequence of point distributions,Fig.10 shows CPU times for setting up and solving the systemmatrices.One can observe that the LSQ method has lower computational efficiency due to more neighboring points resulting in large discrete stencils,while the FDDM greatly benefits from the sparsity of its stencils when the number of unknowns increases.

Fig.9 LSQ vs.FDDM,left：error,right：error

Fig.10 Computational cost for setting up and solving system matrices：LSQ vs.FDDM

Example 4 Solve a discontinuous coefficient problem（originally coming from Ref.[11]）in the form

whereΩ＝[0,1]×[0,1],the conductivityκis discontinuous and given by

and f,g are directly deduced from the exact solution

Fig.11 Random point distribution

It is obvious that this solution and its normal component of flux are continuous at x＝0.5.

The point distribution is almost the same as that in Example 3,but at x＝0.5 uniformly distributed points are placed to coincide with the multimedia interface（shown in Fig.11）.

Here,we takeκ＝10－3,10－6,and the corresponding results are shown in Fig.12 and Tables 3,4.

Figure 12 displays results on the interface for the caseκ ＝10－3,which indicates that the numerical solutions and energy fluxes are very close to the exact values；Table 3gives corresponding convergent results,where NIis the total number of interface points.It is obvious that both numerical solutions and energy fluxes on the interface are second-order accurate.Results for the caseκ＝10－6are similar to the above case,hence they are not displayed here.Table4 shows us that the solutions to the discontinuous coefficient problem with different coefficients are almost second-order accurate,which verifies good performance of the proposed method.

Fig.12 Results on interface asκ＝10－3,and N＝289,NI＝15,left：solutions：right：energy flux

Table 3 Convergence results of solutions and energy fluxes on interface of discontinuous coefficient problem asκ＝10－3

Table 4 Convergence results for discontinuous coefficient problem

4 Conclusions

We present an approach for numerically solving nonlinear diffusion equations in the FDDM regime.Taking advantage of a proper method for selecting steady neighboring point set in the procedure,the approach leads tominimal stencils,avoiding singularity in computing process.

Moreover,when the diffusion conductivity is discontinuous,discrete points are placed on the interface,and a scheme based on five-point formulae of the FDDM is proposed for discretizingmultimedia interface condition.In consequence,approximation to energy fluxes on interface is second-order accurate.

Finally,the approaches are demonstrated to have good accuracy and efficiency by numerical exampleswith different computational domains and different point distributions.

[1] Shen L J,Lv G X,Shen Z J.A finite pointmethod based on directional differences[J].SIAM JNumer Anal, 2009,47（3）：2224－2242.

[2] Seibold B.Minimal positive stencils inmeshfree finite differencemethods for the Poisson equation[J].Comput Methods Appl Mech Eng,2008,198（3－4）：592－601.

[3] Davydov O,Oanh D T.Adaptivemeshless centres and RBF stencils for Poisson equation[J].JComput Phys, 2011,230（2）：287－304.

[4] Larsson E,Lehto E,Heryudono A,Fornberg B.Stable computation of differentiation matrices and scattered node stencils based on Gaussian radial basis functions[J].SIAM JSci Comput,2013,35（4）：2096－2119.

[5] Iliev O,Tiwari S.A generalized（meshfree）finite differnce discretization for elliptic interface problems, numericalmethods and applications[C]∥Hutchison D.Lecture Notes in Computer Science.German：Springer, 2003：488－497.

[6] Sophya T,Silva A D,Kribèchea A.A space-time meshlessmethod that removes numerical oscillations when solving PDEswith high discontinuities[J].Numerical Heat Transfer,Part B,2012,62（1）：50－70.

[7] Oñate E,Idelsohn S,ZienkiewiczOC,Taylor R L.A finite pointmethod in computingmechanicsapplication to convective transport and fluid flow[J].Internat JNumer Methods Engrg,1996,39：3839－3866.

[8] Lv GX,Shen L J,Shen Z J.Study on finite pointmethod[J].Chinese Journal of Computational Physics,2008, 25（5）：505－524.

[9] Lv G X,Shen L J.A finite pointmethod based on directional derivatives for diffusion equation[C]∥Cemal A.World Academy of Science Engineering and Technology.Singapore：International Scientific Research and Experimental Development,2011：211－216.

[10] Lv G X,Shen L J,Shen Z J.Numerical methods for energy flux of temperature diffusion equation on unstructured meshes[J].Int JNumer Meth Biomed Engng,2010,26（5）：646－665.

[11] Shashkov M,Steinberg S.Solving diffusion equationswith rough coefficients in rough grids[J].JComput Phys, 1996,129：383－405.

A Finite Directional Difference M eshless M ethod for Diffusion Equations

LV Guixia, SUN Shunkai
（Laboratory ofComputational Physics,Institute of Applied Physics and Computational Mathematics,P.O.Box 8009-26,Beijing 100088,China）

1001-246X（2015）06-0649-13

An approach for numerically solving nonlinear diffusion equations on 2D scattered point distributions is developed with finite directional difference method.The approach yields stencils of minimal size using five neighboring points.And coefficients of discretization have explicitexpressions.A scheme employing five-point formulae is proposed to discretizemultimedia interface condition for discontinuous problems in which approximation to flux on interface is second-order accurate.The discretizationmethods show good performance in numerical exampleswith different computational domains and different point distributions.

meshless；finite directional differencemethod；nonlinear diffusion equations；multimedia interface；minimal stencil

O241.82 Document code：A

Received date：2014－12－17；Revised date：2015－02－05

Foundation items：Supported by National Natural Science Foundation of China（11371066,11372050）and Foundation of Laboratory of Computational Physics

Biography：Lv Guixia（1972－）,female,Dr.,professor,engaged in numerical solution of partial differential equations,E-mail：lvguixia＠126.com