Higher-Order Numerical Solution of Two-Dimensional Coupled Burgers’ Equations ()
Received 10 October 2015; accepted 13 June 2016; published 16 June 2016

1. Introduction
The Burgers’ equation is an important non-linear parabolic partial differential equation widely used to model several physical flow phenomena in fluid dynamics teaching and in engineering such as turbulence, boundary layer behaviour, shock wave formation and mass transport [1] . Due to its wide range of applicability, several researchers, both scientists and engineers, have been interested in studying the properties of the two-dimensional coupled Burgers’ equation (TDCBE) using various numerical techniques.
There exist many different explicit and implicit numerical schemes with second-order approximation in the space variables, and a first or second-order approximation in the time variable. For example in [2] - [5] , the Crank-Nicolson scheme using the different fully/semi implicit finite-difference methods for the numerical solution of the TDCBE was applied. The implicit logarithmic and local discontinuous Galerkin finite-difference methods for the numerical solution of the TDCBE are proposed in [6] [7] . Also in [4] an explicit scheme using the finite-difference method was applied.
The implicit finite-difference methods with forth-order approximation in the space variables, and a second- order approximation in the time variable are proposed in [8] [9] . These methods based on the Crank-Nicolson scheme with Padé approximation of the finite-difference operator, and hybrid Crank-Nicolson Du Fort and Frankel scheme, respectively. However, the implicit methods on each time layer required to solve an algebraic system. In multidimensional case of the TDHE, it requires large calculation time for solving the algebraic systems till final time layer
, even taking into account the band structure of the matrices [10] .
The aim of the present paper is to construct a new stable and explicit finite-difference scheme to solve the two-dimensional heat equation (TDHE) with Robin boundary conditions. The proposed scheme has a fourth- order approximation in the space variables, and a second-order approximation in the time variable. We deve- loped the proposed scheme for solving a numerical solution of the TDCBE, which comes into the TDHE by the application of the Hopf-Cole transformation.
It is known that the time step of the explicit time-marching schemes must satisfy the so-called Courant- Friedrichs-Lewy condition, which usually enforces a limiting constraint on the time step. However, the main advantages of our explicit scheme considered are saving computing time and memory, and making paralle- lization easier compared to the other numerical methods applied to the TDHE.
The accuracy of the proposed numerical scheme is examined by comparing the numerical and exact solutions of the several TDCBE. The numerical results are found in good agreement with exact solutions for a wide rang of the Reynolds number and confirm the approximation orders of the proposed scheme. We also compared the efficiency of the proposed scheme and implicit fourth-order finite-difference method [8] . Both methods are comparable by the convergence of the solutions and total calculation times.
The structure of the paper is as follows. In Section 2, we present reductions of the TDCBE to a TDHE. The explicit fourth-order accurate finite-difference scheme for solving the TDHE and the fourth-order accurate finite-difference schemes for solving the TDCBE are given in Section 3. Numerical results are discussed in Section 4.
2. The Statement of the Problem
The TDCBE is given by
(1)
(2)
subject to the initial conditions
(3)
(4)
and Dirichlet boundary conditions
(5)
(6)
and the potential symmetry condition
(7)
Here
is the computational domain, and
is its boundary;
and
are the velocity components to be determined;
,
,
and
are known functions;
is the Reynolds number.
Using the Hopf-Cole transformations [7] [8]
(8)
(9)
Equations (1) (2) are reduced to the TDHE
(10)
where
is an arbitrary function depending on t only.
Theorem 1 [7] . Let
be the solution of Equation (10), the functions
and
are defined in Equations (8) and (9). Then
and
are independent of the function
.
By the above theorem, we can choose
, and Equation (10) is simplified to
(11)
The initial conditions (3), (4) and boundary conditions (5), (6) lead to
(12)
(13)
(14)
(15)
(16)
respectively. Here the initial-condition function
has the form [7] [8]
(17)
Thus, the TDCBE (1)-(6) are fully reduced to TDHE (11) with the initial and boundary conditions (12)-(16).
3. The Fourth-Order Accurate Explicit Finite-Difference Scheme
For the TDHE (11)-(16), we consider the following eleven-points explicit finite-difference scheme:
(18)
Here and throughout the work,
is the approximate solution of the
at the mesh point (
,
,
), where
and
are spatial steps by x and y,
is a time step, A, B, C and D are unknown coefficients. Let
be the error function. In this term the scheme (18) has the form
(19)
where
is an approximation error and
(20)
We suppose that the solution of Equations (11)-(16) is a sufficiently smooth function with respect to x, y and t. Using the Taylor expansions of
,
,
and
at the point
, and an identity
(21)
we have
(22)
Equating the coefficients of the partial derivatives to zero in (22), we obtain following system of equations
(23)
The above system has a unique solution if
:
(24)
where
. Using (24) and the higher-order Taylor expansions of the
at the point
we obtain
(25)
One can see that, the condition [11]
(26)
does not improve the order of the scheme (18), i.e., the truncation error of the scheme (18) is of the order of
for any
. If
, the scheme (18) is simplified to the canonical form:
(27)
To find the stability condition of the scheme (27), we seek the partial solution in the form:
(28)
From (27) we have
(29)
![]()
![]()
We have following theorem:
Theorem 2 [10] Let
, b and c are real numbers. Then roots of a quadratic equation
satisfy the condition
if and only if
(30)
Using the conditions (30), we obtain
(31)
The last inequality is true for any A under condition
or
(32)
The scheme (27) at
(or
) is a two-layer scheme in time, while at
(or
) is a three- layer one. Hence, if
in order to find
at level two, two values
and
are required. Using
the Taylor expansion of
at point
and Equation (11) we obtain
(33)
From the initial condition (12) and Taylor expansion (33), we find
with the accuracy ![]()
(34)
From the Robin boundary conditions (13)-(16) using the asymmetric fourth-order finite-difference approxi- mations of the first spatial derivative [12] , we find
,
,
and ![]()
(35)
Now we need to calculate values of the vertex points
,
,
and
. Each value of these points can be calculated using the boundary conditions (13)-(16) and a similar formula to (35) by direction x or y or a middle value of the values by the both directions. Below we presented formulas which used only the boun- dary conditions (13), (14):
(36)
Thus, we find
for
and
by Formulas (27), (35) and (36).
The higher-order finite-difference schemes presented in our previous papers [13] - [15] are applied for finding solutions
,
of Equations (1)-(6). We used the following fourth-order finite-difference scheme [13] :
(37)
with boundary conditions
(38)
Here
is an approximate solution of
. In a similar way we obtain
(39)
with boundary conditions
(40)
Here
is an approximate solution of
. The three-diagonal systems (37), (38) and (39), (40) are solved by the efficient elimination method [16] .
4. Numerical Results
Two exact solvable TDCBEs (1)-(6) are solved to show demonstrate the efficiency and robustness of the pro- posed schemes. To analyze the convergence of the proposed schemes, we used the maximum absolute errors of the solutions
,
and
:
![]()
(41)
![]()
The order (or Runge coefficient) of convergence of the proposed schemes is defined by the double-crowding spatial grids
(42)
The initial and boundary conditions (3), (4), (12) and (5), (6), (13)-(16) for the solutions
,
,
are taken from the analytical solutions. The computational domain is
.
All calculations were performed in double-precision arithmetic on a AMD Phenom II X6 processor using Intel FORTRAN Compiler.
Example 1. ( [7] [8] ). In this example, we solve the two-dimensional Burgers Equations (1), (2), for which the exact solutions are
(43)
The initial and boundary conditions are taken from the exact solutions. We solve the TDHE (11), for which the exact solution is
(44)
The initial (12) and boundary conditions (13)-(16) are taken from the exact solutions.
The convergence of the solutions
,
,
versus the inverse of the Reynolds number
and the numbers of grid N, M are presented in Table 1. To show that the method is fourth-order accurate in space, we fix the time step
as 0.0001 that in each the numbers of grid
holds the stability condition (32). The orders of convergence of the proposed schemes are consistent with the theoretical expectations
.
In Table 2 we compared the efficiency of the proposed scheme with the time step
and implicit the
fourth-order finite-difference method with the time step
[8] at
. From this Table the both methods are comparable by the convergence of the solution
and total calculation times.
Example 2 ( [7] [8] [17] ). In this example, we solve the two-dimensional Burgers Equations (1), (2), for which the exact solutions are
![]()
![]()
Table 1. The convergence of the solutions
,
,
and their corresponding orders of convergence for the Example 1 at T = 1,
versus the parameter
and the numbers of grid N, M. The first column shows the parameter
, the second ones displays the numbers of grid N, M. The third, fifth and seventh columns display the maximum absolute error
,
,
, while the second, forth, and sixth columns present their orders of convergence, respectively. The factor x in the brackets denotes 10x.
(45)
The initial and boundary conditions are taken from the exact solutions. We solve the TDHE (11), for which the exact solution is
(46)
The initial (12) and boundary conditions (13)-(16) are taken from the exact solutions.
The solutions (45) are so-called shock solutions of the TDCBE. It is well known that one of the difficulties in
![]()
Table 3. The same as in Table 1, but for for the Example 2. The factor x in the brackets denotes
.
solving Burgers’ equations is that shock of the solution may occur after some time, even if the initial functions are smooth. When the characteristic curves of Burgers’ equation cross, a shock of the solution occurs. A robust and accurate numerical algorithm should be able to capture the shock and the numerical solution should exhibit the correct physical behavior. From Table 3, we observe that for small values of
, one must consider a large numbers of N and M to obtain proper solutions. Here our proposed scheme works well, and the orders of convergence of the proposed schemes are consistent with the theoretical expectations
.
5. Conclusion
The proposed higher-order finite-difference schemes are easy for implementation and can be used for a numerical solution of two-dimensional coupled Burgers’ equation with higher accuracy. The numerical results show that the variation in the values of the Reynolds number does not adversely affect the numerical solutions. Since all numerical results obtained by the above methods show a reasonably good agreement with the exact one for modest values of
, and also exhibit the expected convergence as the mesh size is decreased, the proposed methods can be considered to be competitive and worth recommendation.
Acknowledgements
The work was supported partially by the Foundation of Science and Technology of Mongolia, and the JINR theme 05-6-1119-2014/2016 “Methods, Algorithms and Software for Modeling Physical Systems, Mathematical Processing and Analysis of Experimental Data”.