International Journal of Scientific & Engineering Research, Volume 6, Issue 4, April-2015 215

ISSN 2229-5518

Numerical Solutions of Two-Dimensional Burgers’ Equations

Vildan Gülkaç-Department of Mathematics, Faculty of Science and Arts, Kocaeli University, Kocaeli/Turkey vgulkac@kocaeli.edu.tr

Abstract— Two-dimensional Burgers’ equations are reported various kinds of phenomena such as turbulence and viscous fluid. In this paper, we illustrate the LOD method for solving the two-dimensional coupled Burgers’ equations. W e extend our earlier work [1] and a stability analysis by Fourier method of the LOD method is also investigated. The computational results obtained by present method are in excellent agreement with earlier results. Present method can be easily implemented for solving nonlinear problems evolving in several branches of engineering and science.

Keywords— Burgers’ equations, Finite-difference, LOD method, Reynolds number

——————————  ——————————

1. Introduction

Burgers’ equation is a fundamental partial differential equation from fluid mechanics. It has been found to describe various kinds of phenomena such as modeling of dynamics, heat conduction, acoustic waves, a mathematical model of turbulence, and the approximate theory of flow through a shock wave traveling in a viscous fluid.
Consider two-dimensional coupled nonlinear
Burgers’ equations taken from [2]
Fletcher [5] has discussed the comparison of a number of different numerical approaches. Goyon [6] used several multilevel schemes with ADI. A.R Bahadır [2] has applied a fully implicit method. V.K. Srivastava et al. [7] has applied a Crank-Nicolson scheme, El-Sayed and Kaya has applied a decomposition method [8], Zhu et all. [9] developed numerical solutions by discreate Adomian decomposition method.
In this paper, Locally One Dimensional (LOD)
method is used to solve two-dimensional Burgers’ equations. Computed results are compared with analytical and other numerical results.

𝜕𝑢 + 𝑢 𝜕𝑢 + 𝑣 𝜕𝑢 = 1 �𝜕2 𝑢

𝜕2 𝑢� (1)

𝜕𝑡

𝜕𝑥

𝜕𝑦

𝑅 𝜕𝑥 2 + 𝜕𝑦2

2. LOD Method and Adaptation of Solution Methodology

𝜕𝑣 + 𝑢 𝜕𝑣 + 𝑣 𝜕𝑣 = 1 �𝜕2 𝑣

𝜕2 𝑣� (2)

𝜕𝑡

𝜕𝑥

𝜕𝑦

𝑅 𝜕𝑥 2 + 𝜕𝑦2

Adaptation of solution methodology numerical computations is always active areas for solutions of
subject to the conditions
𝑢(𝑥, 𝑦, 0) = 𝑓(𝑥, 𝑦), 𝑥, 𝑦 ∈ 𝐷
𝑣(𝑥, 𝑦, 0) = 𝑔(𝑥, 𝑦), 𝑥, 𝑦 ∈ 𝐷
and boundary conditions
𝑢(𝑥, 𝑦, 𝑡) = 𝑓1 (𝑥, 𝑦, 𝑡), 𝑥, 𝑦 ∈ 𝜕𝐷, 𝑡 > 0
𝑣(𝑥, 𝑦, 𝑡) = 𝑔1 (𝑥, 𝑦, 𝑡), 𝑥, 𝑦 ∈ 𝜕𝐷, 𝑡 > 0
(3)
(4)
differential equations. The finite-difference methods are
easy to use for numerical solutions, this methods are still
used extensively in practical computations. Recently, there have been a renewed interests in the worked and the application of finite-difference methods for the solutions of the multi-dimensional partial differential equations [2, 7].
Today, new difference methods have been
constantly presented and for multi dimensional problems LOD and ADI scheme get much attention for their unconditional stability and high efficiency. Gülkaç [1, 10] suggested a LOD method for the solution of multi
Here 𝐷 = {(𝑥, 𝑦): 𝑎 ≤ 𝑥 ≤ 𝑏, 𝑎 ≤ 𝑦 ≤ 𝑏} , 𝜕𝐷 denotes the
boundary of D, u(x,y,t) and v(x,y,t) are the velocity
components to be determinant f,g and 𝑓1 , 𝑓2 are known
functions and R is Reynolds number.
The analytic solution of equations (1) and (2) were
proposed by Fletcher using the Hope-Cole transformation
[3]. The numerical solutions of this equation system have been studied by several authors. Jain and Holla [4] developed two algorithms based on cubic spline technique.
dimensional phase change problems.
The two dimensional Burgers’ equations (1) and (2) can be written by splitting it into two one- dimensional equations, respectively eqns. (5), (6) and (7), (8), as seen Gülkaç [1].
The domain of definition was separated into sets of
sub domains defined along the x and y variable mesh such as equations (5) and (6). Each of the equations was then solved over half of the time step used for the complete two-

IJSER © 2015 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 6, Issue 4, April-2015 216

ISSN 2229-5518

dimensional equation using techniques for the one-
dimensional problems.

2 2

Equations (7) and (8) can be written similarly as equation
(9) and (10).
We consider the use of equal mesh spacing

1 𝜕𝑢 = 1 ��𝜕 𝑢�

+ �𝜕 𝑢

1 𝑢 �𝜕𝑢

1 𝑣 �𝜕𝑢

∆𝑥 = ∆𝑦 over each sub domain for the problem. It should be

2 𝜕𝑡

and

2𝑅

𝜕𝑥 2

2

𝑖+1,𝑗

2 � � −

𝑖,𝑗

2

� −

𝜕𝑥 𝑖,𝑗 2

�

𝜕𝑦 𝑖,𝑗

(5)
noted that the solution algorithm possesses high flexibility for using unequal mesh spacing provided that the stability of equations are valid for each spatial mesh spacing separately.

1 𝜕𝑢 = 1 ��𝜕 𝑢

𝜕 𝑢�

� − 1

𝜕𝑢� − 1

𝜕𝑢�

(6)

2 𝜕𝑡

2𝑅

2 �

𝑖,𝑗+1

+ �

𝜕𝑦2

𝑖,𝑗

𝑢 �

2 𝜕𝑥

𝑖,𝑗

𝑣 �

2 𝜕𝑦

𝑖,𝑗

3. Stability of Equations

similarly,

2 2

The basic idea defining von Neumann stability [11]
is that this numerical algorithm used exactly, should limit

1 𝜕𝑣 = 1 ��𝜕 𝑣 �

+ �𝜕 𝑣

1 𝑢 �𝜕𝑣

1 𝑣 �𝜕𝑣

2 𝜕𝑡

and

2𝑅

𝜕𝑥 2 𝑖+1,𝑗

2

� � −

𝑖,𝑗 2

2

�

𝜕𝑥

−

𝑖,𝑗 2

�

𝜕𝑦 𝑖,𝑗

(7)
the amplification of all elements of initial conditions.
The Burgers’ equation express the initial values at
the mesh points along t=0 in terms of a finite Fourier series,
then regards the growth a function that reduces to this series for t=0 by a variable separable method

1 𝜕𝑣 = 1 ��𝜕 𝑣

𝜕 𝑣�

� − 1

𝜕𝑣 � − 1

𝜕𝑣 �

(8)

2 𝜕𝑡

2𝑅

𝜕𝑦2 �

𝑖+1,𝑗

+ �

𝜕𝑦2

𝑖,𝑗

𝑢 �

2 𝜕𝑥

𝑖,𝑗

𝑣 �

2 𝜕𝑦

𝑖,𝑗

indistinguishable to that commonly used for solving partial
differential equations. The Fourier series can be formulated
In this method, we replace all spatial derivatives with the
in complex exponential form [11].
In order to show the von Neumann stability of the
average of their values at the 𝑛 and 𝑛 + 1�2
time levels and
present method, we replace
then substitute the central finite form all derivatives. Eqns.
(5), (6) and (7), (8) can be written as

𝑛+1�2 𝑛

𝑖,𝑗

𝑢𝑖,𝑗,𝑛 = 𝑢𝑝,𝑞,𝑟 = 𝑒𝑖𝛽𝑝ℎ 𝑒𝑖𝛽𝑞𝑘 𝑒𝑖𝛼𝑡 = 𝑒𝑖𝛽𝑝ℎ 𝑒𝑖𝛽𝑞𝑘 𝜉𝑟

and 𝑟𝑥 = 𝑅𝑥 , 𝑟𝑦 = 𝑅𝑦 .

1 � 𝑖,𝑗 1

2 − 1 𝛿 4 + 1

6 ∓ ⋯ � �𝑢𝑛 + 𝑢𝑛 �

2 ∆𝑡

2𝑅∆𝑥 2 �𝛿𝑥

12 𝑥 90 𝛿

𝑖+1,𝑗

𝑖,𝑗

Where 𝜉 = 𝑒𝛼𝑡 , 𝑖 = √−1 and 𝛼 in general is a complex


− 1 𝑢 �𝜕𝑢�


− 1 𝑣 �𝜕𝑢�
constant, 𝜉 is often called amplification factor [11]. The

2 𝜕𝑥 𝑖,𝑗 2

or

𝜕𝑦 𝑖,𝑗

finite difference equations will be stable by von Neumann
definition if |𝜉| ≤ 1 [ 11].
Equation (9) can be written as
�1 + � 1 − 1 𝑟 � 𝛿 2� 𝑢

𝑛+1�2 = �1 + � 1 + 1 𝑟 � 𝛿 2� 𝑢𝑛

( ) ⁄ ⁄

12 2 𝑥

𝑥 𝑖,𝑗

12 2 𝑥

𝑥 𝑖 ,𝑗

� 1 − 1 𝑅 � 𝑒𝑖𝛽 𝑝+1 ℎ 𝑒𝑖𝛽𝑞𝑘 𝜉𝑟+1 2 + (𝑅
− 5)𝑒𝑖𝛽𝑝ℎ 𝑒𝑖𝛽𝑞𝑘 𝜉𝑟+1 2 +

12 2 𝑥 𝑥

� 1 − 1 𝑅 � 𝑒𝑖𝛽(𝑝−1)ℎ 𝑒𝑖𝛽𝑞𝑘 𝜉𝑟+1⁄2 = 1

𝑖𝛽𝑝ℎ

𝑖𝛽𝑞𝑘 𝑟

12 2 𝑥

�(−5 − 𝑅 )𝑒

𝑅

𝑒 𝜉 +

𝑛 −𝑢𝑛

𝑛 −𝑢𝑛

� 1 + 1 𝑅 � 𝑒𝑖𝛽(𝑝+1)ℎ 𝑒𝑖𝛽𝑞𝑘 𝜉𝑟 + � 1 + 1

𝑖𝛽(𝑝−1)ℎ

𝑖𝛽𝑞𝑘 𝑟

− 1 𝑢𝑛 𝑢𝑖,𝑗

𝑖−1,𝑗

1 𝑣 𝑛 𝑢𝑖,𝑗

𝑖−1,𝑗

12 2 𝑥

𝑅 � 𝑒

12

𝑒 𝜉 � −

2 𝑖 ,𝑗

−

∆𝑥 2

𝑖,𝑗

(9)

∆𝑦

1

2∆𝑥

𝑢0 �𝑒𝑖𝛽𝑝ℎ 𝑒𝑖𝛽𝑞𝑘 𝜉𝑟 − 𝑒𝑖𝛽(𝑝−1)ℎ 𝑒𝑖𝛽𝑞𝑘 𝜉𝑟 � −
for ∀𝑖 = 1, … , 𝑁, ∀𝑗 = 1, … 𝑁, 𝑟𝑥 = ∆𝑡⁄∆𝑥2 , 𝑟𝑦 = ∆𝑡⁄∆𝑦2
𝑡𝑛+1⁄2 → 𝑡𝑛+1 and equation (6) can be written as,
then

1 𝑣 �𝑒

2∆𝑦

𝑖𝛽𝑝ℎ

𝑒𝑖𝛽𝑞𝑘
𝜉𝑟
− 𝑒

𝑖𝛽𝑝ℎ

𝑒𝑖𝛽(𝑞−1)𝑘
𝜉𝑟 � (11)

1 1

2 𝑛+1

𝑢0 𝑎𝑛𝑑 𝑣0 initial conditions then we let 𝑢0 = 𝑣0 = 0 and
∆𝑥 = ∆𝑦, and division eqn.(11) by 𝑒𝑖𝛽𝑝ℎ 𝑒𝑖𝛽𝑞𝑘 𝜉𝑟 equation (11)
can be written as,
�1 + �12 − 2 𝑟𝑦 � 𝛿𝑦 � 𝑢𝑖,𝑗

1 1

2 𝑛+1⁄2

1� 1 1

−𝑖𝛽ℎ

1 1

𝑖𝛽𝑝ℎ

= �1 + �12 + 2 𝑟𝑦 � 𝛿𝑦 � 𝑢𝑖,𝑗

𝜉 2 �� − 𝑅 � 𝑒

12 2

+ (𝑅𝑥 − 5) + � − 𝑅𝑥 � 𝑒 � =

𝑛+1⁄2 − 𝑢𝑛+1⁄2

1 1 1 𝑅 � 𝑒−𝑖𝛽𝑝ℎ − (𝑅

+ 5) + � 1

1 𝑅 � 𝑒𝑖𝛽𝑝ℎ �

− 1 𝑢𝑛 𝑢𝑖,𝑗
2 𝑖 ,𝑗
∆𝑥

𝑖−1,𝑗 −

�� +

𝑅 12 2

and then
+

12 2

𝑛+1⁄2 −𝑢𝑛+1⁄2

1 𝑣 𝑛 𝑢𝑖,𝑗

𝑖,𝑗−1

1 � 1 +1𝑅𝑥 �2 cos 𝛽ℎ−(𝑅𝑥 +5)

2 𝑖,𝑗

(10)

∆𝑦

𝜉 �2 =

12 2

1

𝑅�� + 𝑅 � 2cos 𝛽ℎ�+(𝑅 −5)

12 2

IJSER © 2015 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 6, Issue 4, April-2015 217

ISSN 2229-5518

𝑢(𝑥, 0, 𝑡) = 1 + 𝑠𝑖𝑛𝜋𝑥, 𝑢(𝑥, 0.5, 𝑡) = 𝑠𝑖𝑛𝜋𝑥

1

� 0 ≤ 𝑥 ≤ 0.5, 𝑡 ≥ 0
let 𝜉 �2 = 𝐾 the condition is |𝐾| ≤ 1.This required. Clearly,
0 < 𝐾 ≤ 1 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑅𝑥 > 0 𝑎𝑛𝑑 𝑎𝑙𝑙 𝛽 𝑎𝑛𝑑 𝑅 =
𝑅𝑒𝑦𝑛𝑜𝑙𝑠 𝑛𝑢𝑚𝑏𝑒𝑟 ≥ 1, |𝐾| ≤ 1.
Therefore, equations unconditionally stable.
Similarly, it is easily shown by same method that equations (10) and (7), (8) are unconditionally stable.

4. Numerical Examples and Conclusions

4.1 Problem I

The exact solutions of Burgers’ equations (1) and (2)can be generated by using the Hope-Cole transformation [3] which are:
𝑣(𝑥, 0, 𝑡) = 𝑥, 𝑣 = 𝑥, 0.5, 𝑡) = 𝑥 + 0.5
The numerical methodology used is similar to that of
Gülkaç [1] and Gülkaç and Öziş [10]. We presented the
numerical method for solving two-dimensional Burgers’ equations using the LOD method and then substituted Douglas-like equation form for all derivatives [1, 11].
Equations (1), (2) are discredited using the LOD method. The stability analysis of the scheme is also investigated and the scheme is therefore unconditionally stable. The accuracy of the numerical solutions indicates that the method is well suited for the solution of two- dimensional non-linear Burgers’ equations.

Table 1

Comparison of numerical values of u for R=500 at t=0.625


𝑢(𝑥, 𝑦, 𝑡) = 3 − 1 ,

4 4[1+exp((−4𝑥+4𝑦−𝑡)𝑅⁄32]

𝑣(𝑥, 𝑦, 𝑡) = 3 + 1 ,

4 4[1+exp((−4𝑥+4𝑦−𝑡)𝑅⁄32]

Here the computational domain is taken as a square
domain 𝐷 = {(𝑥, 𝑦): 0 ≤ 𝑥 ≤ 1, 0 ≤ 𝑦 ≤ 1} . The initial and
boundary conditions are taken from the exact solutions.
The numerical computations are performed using uniform
grid, with a mesh width ∆𝑥 = ∆𝑦 = 0.05.
From Tables 1-4, it is clear that the results from the
present study are in good agreement with the exact solution
for different values of Reynolds number (R) and some typical mesh points demonstrate that the present scheme achieves similar results as those of Jain and Holla [4], Bahadır [2], Srivastava et all [7].

4.2. Problem II

Here the computational domain is taken as
𝐷 = {(𝑥, 𝑦): 0 ≤ 𝑥 ≤ 0.5, 0 ≤ 𝑦 ≤ 0.5}
and Burgers’ equation (1) and (2) are taken with the initial
conditions,
𝑢(𝑥, 𝑦, 0) = 𝑠𝑖𝑛𝜋𝑥 + 𝑐𝑜𝑠𝜋𝑦
𝑣(𝑥, 𝑦, 0) = 𝑥 + 𝑦 � 0 ≤ 𝑥 ≤ 0.5, 0 ≤ 𝑦 ≤ 0.5
and boundary conditions,

𝑢(0, 𝑦, 𝑡) = 𝑐𝑜𝑠𝜋𝑦, 𝑢(0.5, 𝑦, 𝑡) = 1 + 𝑐𝑜𝑠𝜋𝑦

𝑣(0, 𝑦, 𝑡) = 𝑦, 𝑣(0.5, 𝑦, 𝑡) = 0.5 + 𝑦 � 0 ≤ 𝑦 ≤ 0.5, 𝑡 ≥ 0,

Table 2

Comparison of numerical values of v for R=500 at t=0.625

IJSER © 2015 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 6, Issue 4, April-2015 218

ISSN 2229-5518

Table 3

Comparison of numerical values of u for R=50 at t=0.625

Table 4

Comparison of numerical values of v for R=50 at t=0.625

5. References

[1] V. Gülkaç, A numerical solution of two-dimensional fusion problem with convective boundary conditions, International Journal for Computational Methods in Engineering Science and Mechanics, 11 (2010) 2-26.
[2] A.R. Bahadır, A fully implicit finite-difference scheme for two-dimensional Burgers’ equations, Appl. Math. Comput. 206 (2008) 131-137.
[3] C.A.J. Fletcher, Generating exact solutions of two- dimensional Burgers’ equation, Int. J. Numer. Meth. Fluids
3 (1983) 213-216.
[4] P.C. Jain, D.N. Holla, Numerical solutionof coupled
Burgers’ equations, Int. J. Numer. Meth. Eng. 12 (1978) 213-
222.
[5] C.A.J. Fletcher, A comparison of finite element and finite
difference solution of the one-and two-dimensional
Burgers’ equations, J. Comput. Phys. 51 (1983) 159-188.
[6] O. Goyon, Multilevel schemes for solving unsteady
equations, Int. J. Numer. Meth. Fluids 22 (1996) 937-959.
[7] M.K. Srivastava, M. Tamsir, U. Bhardwaj, Y. Sanyasiraju,
Crank-Nicolson scheme for numerical solutions of two dimensional coupled Burgers’ equations, 2 (2011) 1-7.
[8] S.M. El-Sayed, D. Kaya, On the numerical solution of the
system of two-dimensional Burgers’ equations by decomposition method, 158 (2004) 101-109.
[9] H. Zhu, H. Shu, M. Ding, Numerical solutions of two- dimensional Burgers’ equations by discreate Adomian decomposition method, Computer and Mathematics with Applications, 60 (2010) 840-848.
[10]V. Gülkaç, T. Öziş, On a LOD method for solution of
two-dimensional Fusion problem with convective boundary conditions, International Communications in Heat and Mass Transfer, 31 (2004) 597-606.
[11] G.D. Smith, Numerical Solutions of Partial Differential Equations: Finite Difference Methods, Clarendon Press, Oxford, 1985.

IJSER © 2015 http://www.ijser.org