Journal of Applied Mathematics

Journal of Applied Mathematics / 2018 / Article

Research Article | Open Access

Volume 2018 |Article ID 3416860 | 18 pages |

Numerical Solution to Coupled Burgers’ Equations by Gaussian-Based Hermite Collocation Scheme

Academic Editor: Mustafa Inc
Received08 May 2018
Revised21 Aug 2018
Accepted04 Sep 2018
Published27 Sep 2018


One of the most challenging PDE forms in fluid dynamics namely Burgers equations is solved numerically in this work. Its transient, nonlinear, and coupling structure are carefully treated. The Hermite type of collocation mesh-free method is applied to the spatial terms and the 4th-order Runge Kutta is adopted to discretize the governing equations in time. The method is applied in conjunction with the Gaussian radial basis function. The effect of viscous force at high Reynolds number up to 1,300 is investigated using the method. For the purpose of validation, a conventional global collocation scheme (also known as “Kansa” method) is applied parallelly. Solutions obtained are validated against the exact solution and also with some other numerical works available in literature when possible.

1. Introduction

Fluid is known to dominate most part of the planet and can be modelled mathematically by the well-known Navier-Stoke equation. To understand its behavior under a certain set of conditions, one seeks solutions (if any or if possible) of the equation and this is not an easy task. The combination of challenging factors, nonlinearity, coupling, and transient nature, of the equation makes it one of the most complex systems both analytically and numerically.

Acknowledged as a simple case of the Navier-Stoke equation, the famous “Burgers” equations, a system of equations that describes the interaction between two crucial physical manners of nature, convection, and diffusion, are used to model variety of applications. This includes flows through a shock wave traveling in viscous fluid, the phenomena of turbulence, sedimentation of two kinds of particles in fluid suspensions under the effect of gravity, and also dispersion of pollutants in rivers (see [13]).

The standard form of this type of equation is expressed asThe ratio of inertial forces to viscous forces is represented by the Reynolds number, Re, and the x- and y-velocity components are noted by and , respectively.

For a domain with its boundary , the above system is subject to the initial conditionsand the boundary conditionsHere , and are known functions.

In 1983, Fletcher [4] applied the Hopf-Cole transformation and successfully solved the equation analytically. Ever since, the coupled two-dimensional Burgers’ Equations have been attracting a number of investigations both numerically and analytically. Some of the recent computational researches include the application of a fully implicit finite-difference (see [5]), by Adomain-Pade technique (see [6]), by adopting variational iteration method (see [7]), and by applying a mesh-free technique (see [8, 9]). Even more recently, more numerical works have been successfully carried out and nicely documented in literature: the application of the dual reciprocity boundary element method (see [10, 11]), POD/discrete empirical interpolation method- (DEIM-) reduced-order model (ROM) (see [12]), and a combination between the homotopy analysis method and finite differences (see [13]).

Amongst the well-known numerical scheme, finite volume, finite difference, and finite element method that have been invented, developed, and applied in a wide range of science and engineering problems, a rather young idea was discovered and has been categorized as “meshless/mesh-free” methods. The methods under this category have recently become promising alternative tools for numerically solving variety of science and engineering problems. Generally, these meshless schemes can be grouped into two main classes:(i)Strong forms: the finite point method (see [14]), The hp-meshless cloud method (see [15]), the collocation method (see [16]), and references therein.(ii)Weak forms: the diffuse element method (see [17]), the element-free Galerkin method (EFGM; see [18]), the point interpolation method (see [19]), and references therein.

This depends on the requirement for integration. In this work, particularly, the main attention is paid to the development of the so-called “Collocation Method” characterized in the strong form family. The method is truly mesh-free meaning that no mesh is at all required at any point of the whole computing process. Amongst several versions of its further improvements, this work focuses on an application of the Hermite form of the conventional collocation where the involved differential operator is applied twice, explained in more detail in Section 2.

The paper is organized as follows. Section 2 provides the necessary mathematical background concerning the Gaussian radial basis function and also the fundamental idea of collocation mesh-free schemes in general. The implementation of the schemes to the couple-Burgers equations is then detailed in Section 3. Section 4 demonstrates the numerical experiments and discusses general aspect of the application before some interesting findings are listed in Section 5.

2. The Gaussian-Based Collocation Meshless Method

Letting be a point in a -dimensional Euclidean space , for any point that depends only on the distance, , from the fixed , then the function that is written as is called a Radial Basis Function (RBF), with being called the center and some nonnegative parameter .

Among different types of RBFs invented, improved, and applied nowadays, in this work, one of the most widely used forms of radial functions, known as “Gaussian,” is focused on. The standard form of this type of RBF is expressed aswhere and is called “shape parameter,” that is related to the variance of the normal distribution function by . Under the context of collocation schemes, this relation is linked to the Euclidean distance function, , via the following multivariate function:for a fixed center and . It can then be seen that the connection between and is as follows:

The collocation scheme starts with considering the following elliptical partial differential equation defined on a bounded and connected domain :where is the domain boundary containing two nonoverlap sections: and , with . These differential operators, and are applied on the domain and the two boundary sections, respectively. Three known functions can well be dependent of space and/or time. Let be a set of randomly selected points, known as “collocation” or “centers”, on the domain, where are those contained within, where and are those on the boundaries and , respectively.

2.1. The Conventional Kansa Collocation Method

Initially, a collocation scheme writes the approximate solution, , as a linear combination of the basis function , shown in the following form:where are coefficients and is the Euclidean norm. The basis function used now is the radial type as defined previously. Applying and to both domain and boundary sections, satisfying the governing system of equations, allows the system to arrive atwhere , and the known vector expressed is as follows:and by setting to be a matrix with entries for we have Once are obtained by (15), the approximate solution is straightforward yielded. This method is known as “Kansa,” in honor of a great mathematician E. J. Kansa [22] who discovered the idea in 1990. The method has been applied to a wide range of problems ever since [23, 24]. It, nevertheless, has shortcomings as it is known to suffer the problem of unsymmetric interpolation matrix, , and the rigorous mathematical proof of its solvability is still not available [25], and very often produces low quality results particularly in boundary-adjacent region [26, 27].

2.2. The Hermite Collocation Method

In order to alleviate the difficulty mentioned earlier, in 1997, Fasshauer [28] proposed a new way of interpolating by applying the self-adjoint operators and to the governing system of equations and rewriting the approximate solution (14) asThis leads to a new interpolation matrix , shown as follows:

In this study, nevertheless, the problem domain has only one continuous boundary with no differential operator. An application of this Hermite type of collocation method to elastostatic problem was successfully investigated by Leitao [29]. Additionally, some interesting implementations of the scheme to the transient and nonlinear plate problems can be found in [3032].

3. Implementation to Burgers’ Equations

As previously mentioned, the coupled Burgers’ equations can be seen as a simple model that mimics the coupling of fluid flow to thermal dynamics, as well as some other scientific fields. High and low Reynolds numbers play important roles in both theoretically modeling and numerical simulation. In this section, the Hermite collocation is now applied based on the Gaussian type of radial basis function.

Since the principle implementation of the conventional collocation method to nonlinear and time dependent problem has already been documented in some of our previous studies [33, 34], this section is therefore dedicated to the application of the Hermite version of collocation as explained in details in the following subsections.

3.1. Space Discretization with the Hermite Scheme

With Hermite interpolation technique, it begins with writing the approximation of solution for and , respectively, with the same radial basis function , as follows: and Note that it is now assumed that there is only one section of boundary . Therefore, at derivative, this leads toand Therefore, the governing equations are now rewritten, at a center node, aswhere leading their matrix forms expresses as below:Similarly,where Consequently, the systematics of matrices involved in this coupling process is

3.2. Time Discretization

In order to carry out the time-dependent terms, the -order Runge Kutta is applied to both terms appearing in the above equations, that is, by settingThat is,The time progressing of variables involved can be estimated as follows:wherewhere the similar manner is applied for the other case, , and the systems are treated and advanced in time with an ODE method. Then the approximate solution obtained by this final process is expressed in terms of the radial basis function at collocation nodes as follows:

3.3. Gaussian-Based Hermite Collocation

Letting be a radial basis function depending on a distance function, , where , it can be seen that

This leads to the Laplacian expressed as follows:Furthermore, when the Hermite concept of collocation is employed, it is necessary that the Lapacian be applied twice resulting in the following so-called fourth-order biharmonic form; For this work, Gaussian type of radial basis function, is used and its first four orders of derivatives can, be respectively, expressed asfor each pair of center nodes: .

What follows is the brief demonstration of how the self-adjoint operator works.

Hermite starts with the following solution approximation:From (24),(25), (27), and (28), when considering for , we have . That is, For this particular example we have and there is no derivative operation taking place on the boundary. Therefore, the above equation becomesFor , we have =Now, we have and with the absence of the derivative operation on the boundary, , the remaining term of the above equation can be written asSimilarly, By substituting (49), (51), and (52) back into the governing equationsTherefore, for each center , the following is reached:

The whole process explained above is applied to the other equation of the coupled governing equations. Note that all the terms with shall be cancelling each other out so there is no zero-division problem in computing process.

3.4. Computing Setup and Algorithm

Each computation under this work was completed following the computing steps detailed below.

Step 1. Choose collocation or center nodes , on the domain .

Step 2. Specify the desired values of(2.1)The Reynolds number, Re.(2.2)The final time .(2.3)The time step .(2.4)The Gaussian shape parameter .

Step 3. Compute the initial collocation matrices and using the initial Hermite approximation equations (20)-(21).

Step 4. Apply the initial conditions, , to get and from (30), (31), (33), (34), (35), (36), (37), (38), (40), (41), (42), and (43).

Step 5. Compute the solutions for the next time step via the time-discretization explained above.

Step 6. Construct the collocation matrices and in the forms expressed in (30), (31), (33), (34), (35), (36), (37), (38), (40), (41), (42), and (43) using the solutions previously obtained in Step 5.

Step 7. Carry on Steps 5 and 6 until reaching the final time step .

All computing experiments carried out in this study were executed on a computer notebook: Intel(R) Core(TM) i7-5500U CPU @ 2.40GHz with 8.00 GB of RAM and 64-bit operating system.

4. Numerical Experiments and General Discussion

The test case chosen for numerical investigation carried out in this study is the most popular form of Burgers’ equations where the corresponding exact solutions for validation are provided by using a Hopf-Cole transformation nicely carried out by Fletcher [4] in 1983 and are expressed as follows:where both the initial and the boundary conditions to be imposed to the equation system are generated directly from the above exact forms over the domain . With their rich of challenging features, the equations have received interest from several researchers and been treated with variety of numerical techniques as mentioned in Section 1. Nevertheless, it can be clearly seen that, with the reasons previously mentioned, most studies are concerned only with only cases of relatively low Reynolds numbers, i.e., . In this work, the investigation covers a wide range of Reynolds number from low to moderate and, moreover, goes above , where no numerical works have ever reached before.

Solutions obtained from applying the Hermite collocation method are validated by comparing to both those documented in literatures and those produced by the exact formula. Error measurements are carried out using the following error norms:(1)Maximum error ();(2)Root-mean-square error ();

In addition to these norms, the results presented below are sometimes referred to as “Average-Error” and it is to be understood aswhere and are of U- and V-velocity component, respectively.

The whole study is split into two main parts, one concerned with low- to moderate-Reynolds number cases, , and that with relatively high value of Reynolds number, . The so-called Reynolds number plays crucial rules in determining the ratio of forces exerted in the system and this is remained as one of the most challenging tasks to simulate and mimic both numerically and experimentally.

In order to demonstrate the overall effectiveness of the Hermite type of collocation scheme, investigations using the conventional global collocation known as “Kansa” method, (14), are also studied parallelly. The results obtained from using the Hermite type and Kansa are referred to hereafter as “Hermite” and “Conventional Kansa/Kansa,” respectively.

4.1. Low- to Moderate-Reynolds Number Cases

As it is known to be one of the most challenging tasks for any collocation method, finding a suitable shape parameter, , is to be firstly presented. In the past, some attempts to pinpoint the optimal value of involve the classic work of Hardy [35] where it was shown that, by fixing the shape at , where and is the distance from the node to its nearest neighbor, good results should be anticipated. Also, in the work of Franke and Schaback [36] where the choice of a fixed shape of the form with is the diameter of the smallest circle containing all data nodes, can also be a good alternative. In 2000, Zhang et al. [37] demonstrated and concluded that the optimal shape parameter is highly problem dependent. In 2002, Wang and Lui [38] pointed out that, by analyzing the condition number of the collocation matrix, a suitable range of derivable values of can be found. Later in 2003, Lee et al. [39] suggested that the final numerical solutions obtained are found to be less affected by the method when the approximation is applied locally rather than globally. A rather recent work is the selection of an interval for variable shape parameter by Biazar and Hosami in 2016 [40], where a novel algorithm for determining and interval was proposed.

It is to be noted, however, that even though this topic has long been widely studied (see also in [4144]), the majority of studies are done on “Multiquadric RBF” or “Inverse-Multiquadric RBF,” whereas works on Gaussian type are rather rare. A confirmation on problem-dependent nature of Gaussian RBF shape value was confirmed in the work of Davydov and Oanh [45]. A more recent work is the application of the traditional Kansa collocation to linear PDE based on a variable shape parameter called symmetric variable shape parameter (SVSP) carried out by Ranjbar [46].

Even though Gaussian RBF was focused on in these two later works, the scheme does not contain any self-adjoint operator. Moreover, the complications introduced by the structure of the governing equation at hand are known to greatly affect the final choice of the shape. Taking these facts into consideration, with the transient, nonlinear, and coupling nature of Burgers’ equations, this work chose to directly measure the overall errors generated in the system. To complete this task, a large number of numerical experiments were carried out in this work aiming to reach the optimal value of . The impact of shape parameter on the quality of the computed solution is demonstrated in Figures 14. Figures 1 and 2 show that, within the range of , both methods reveal slight differences in the trends of , with some strong fluctuations being found for the case of Kansa. For this particular set of conditions, at Re = 10 with , it can be noted that a more suitable shape parameter may be found at for the Hermite scheme. This observation can also be useful for Kansa case as indicated in the Figure. Figures 3 and 4 indicate clearly that both of the methods become less sensitive to the shape value when it goes beyond . Above this point, it should be observed that remains under when using Hermite scheme whereas stays constantly above .

Information provided previously on the effect of shape parameter leads to decision made for the case with Re = 50. At randomly chosen points, Table 1 provides the solutions of both U- and V-velocity components compared against the corresponding exact values. It can be seen that the solutions are in good agreement for both choices of shape parameter; and . A bigger picture of results obtained in this work is shown in Table 2 and Table 3, with cases of Reynolds number starting from 1 up to 400.