Element-free characteristic Galerkin method for Burgers’ equation

ARTICLE IN PRESS Engineering Analysis with Boundary Elements 33 (2009) 356–362

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Element-free characteristic Galerkin method for Burgers’ equation Xiao Hua Zhang, Jie Ouyang , Lin Zhang Department of Applied Mathematics, Northwestern Polytechnical University, Xi’an 710072, China

a r t i c l e in f o

a b s t r a c t

Article history: Received 7 December 2007 Accepted 9 July 2008 Available online 28 August 2008

A new meshfree method named element-free characteristic Galerkin method (EFCGM) is proposed for solving Burgers’ equation with various values of viscosity. Based on the characteristic method, the convection terms of Burgers’ equation disappear and this process makes Burgers’ equation self-adjoint, which ensures that the spatial discretization by the Galerkin method can be optimal. After the temporal discretization, element-free Galerkin method is then applied to solve the semi-discrete equation in space. Moreover, the process is fully explicit at each time step. In order to show the efficiency of the presented method, one-dimensional and two-dimensional Burgers’ equations are considered. The numerical solutions obtained with different values of viscosity are compared with the analytical solutions as well as the results by other numerical schemes. It can be easily seen that the proposed method is efficient, robust and reliable for solving Burgers’ equation, even involving high Reynolds number for which the analytical solution fails. & 2008 Elsevier Ltd. All rights reserved.

Keywords: Burger’s equation Meshfree method EFCGM Characteristic

1. Introduction Burgers’ equation, which has a convection term, a viscosity term and a time-dependent term, is a nonlinear equation and very similar to the Navier–Stokes equations except pressure gradient term. This equation usually exhibits a ‘‘mixed’’ property, that is, it is parabolic when diffusion term or viscous term is dominative, or else it is hyperbolic. In latter case, shock waves or sharp discontinuities may appear during finite time because of the nonlinearity of Burgers’ equation, even if the initial condition is smooth enough. Such discontinuities are the cause of difficulties that occur in obtaining a solution, thus these properties make Burgers’ equation a proper model for testing numerical algorithm in flows. Many researchers have proposed various kinds of numerical methods for Burgers’ equation. In general, those methods usually fall into the following classes: finite difference method (FDM) [1–4], finite element method (FEM) [5–7], boundary element method (BEM) [8], spectral methods [9], and so on. Up to the present, the development of an innovative, robust and efficient numerical method for seeking accurate numerical solutions of Burgers’ equation with large value of Reynolds numbers (Re) remains as a challenging task [10]. Among the methods mentioned above, the spatial domain where the partial differential governing equations are defined is often discretized into meshes. Generally speaking, the creation of suitable meshes is very essential for acquiring accurate results.

 Corresponding author. Tel.: +86 029 8849 5234; fax: +86 029 8849 1000.

E-mail address: [email protected] (J. Ouyang). 0955-7997/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2008.07.001

However, mesh generation process consumes a lot of time and labor for some problems, especially for discontinuous and high gradient problems, for which these methods will become complicate. The root of these difficulties is the use of mesh in the formulation stage. Recently, a kind of so-called meshfree or meshless methods has developed quickly. Meshfree methods only use a set of nodes scattered within the problem domain as well as a set of nodes scattered on the boundary. To date, there exist many meshfree methods such as element-free Galerkin (EFG) method [11–13], meshless local Petrov–Galerkin (MLPG) method [13], reproducing kernel particle method (RKPM) [12], and so on. The more details of these meshfree methods can refer to Ref. [13]. As a matter of fact, Ouyang et al. [14] developed a series of nonstandard element-free Galerkin methods to solve the onedimensional and two-dimensional Burgers’ equation. But among these methods, the choice of a proper stabilization parameter is difficult. Hashemian et al. [10] applied the gradient reproducing kernel particle method (GRKPM) to solve one-dimensional Burgers’ equation, while Xie et al. [15] used RKPM to solve onedimensional Burgers’ equation. Young et al. [16] applied the Eulerian–Lagrangian method of fundamental solutions (ELMFS) to solve two-dimensional Burgers’ equation. It is well known that the standard Galerkin FEM is not ideally suited to deal with the spatial discretization of Burgers’ equation with large value of Re, nor is the standard EFG method [14]. Many different ideas and approaches have been proposed to overcome the deficiencies of the Galerkin approach in highly convective situations, such as streamline upwind Petrov–Galerkin (SUPG) method [17,18], Galerkin least square (GLS) method [17,18], Taylor–Galerkin (TG) method [18], characteristic Galerkin (CG)

ARTICLE IN PRESS X.H. Zhang et al. / Engineering Analysis with Boundary Elements 33 (2009) 356–362

method [18–20], and so on. Among these methods, the CG method is optimal for transient problems and gives identical stabilizing terms to that derived by the use of SUPG, GLS and other procedures when the time step used is near the stability limit [18]. So in this paper we develop a new meshfree method, which is called element-free characteristic Galerkin method (EFCGM), to solve Burgers’ equation. The paper is organized as follows. In Section 2, the details of the EFCGM in two-dimensional are presented. In Section 3, a onedimensional and two-dimensional numerical experiments including comparisons with the available analytical solution and the numerical results in the literatures are provided. Conclusions are drawn in Section 4.

2. The element-free characteristic Galerkin method

!

qv qv qv 1 q2 v q2 v þu þv ¼ þ qt qx qy Re qx2 qy2

(2)

where Re is the Reynolds number. For a lager value of Re, Burgers’ equation behaves merely as hyperbolic partial differential equation, and the problem becomes very difficult to solve as steep shock-like wave fronts develop. We can write (1) and (2) along the characteristic as [18]     qu 0 q 1 qu q 1 qu ðx ðtÞ; y0 ðtÞ; tÞ  0  0 ¼0 (3) 0 0 qt qx Re qx qy Re qy 

In the EFCGM, the field variable u(x,t) is approximated by moving least squares (MLS) approximation which was initially introduced for data fitting and surface construction in 1981 [11,12]. The MLS approximation uh(x,t) is given by [11–13] uh ðx; tÞ ¼

m X

pj ðxÞaj ðx; tÞ ¼ PT ðxÞaðx; tÞ

(7)

j¼1

where P(x) is a complete polynomial basis of arbitrary order and a(x, t) is coefficient (to be determined), which is the function of the space coordinate x and time t. For the sake of simplicity, linear basis is chosen in the paper. The unknown coefficient a(x, t) in (7) is obtained at any point x by minimizing the following weighted, discrete error norm: n X

wðx  xI Þ½uI ðtÞ  PT ðxI Þaðx; tÞ2

(8)

i¼1

Consider the following unsteady two-dimensional Burgers’ equation ! qu qu qu 1 q2 u q2 u þu þv ¼ þ (1) qt qx qy Re qx2 qy2



2.2. Spatial discretization using meshfree method

J¼

2.1. The temporal discretization

357



qv 0 q 1 qv q 1 qv ðx ðtÞ; y0 ðtÞ; tÞ  0  0 qt qx Re qx0 qy Re qy0



¼0

(4)

As we can see, in the moving coordinates x0 and y0, the convective acceleration terms disappear and diffusion terms are averaged quantities along the characteristic. For such a selfadjoint system, it is known that the standard Galerkin approximation in space is optimal, but the inconvenience of a moving coordinate system is introduced. In order to overcome the difficulty, we can apply local Taylor expansion and coordinate transformation within the time step. This process is fully described in Refs. [18–20] and the explicit form can be written in fixed coordinates as !n unþ1  un qun qun 1 q2 u q2 u ¼ u v þ þ Dt qx qy Re qx2 qy2     Dt q qu qu n Dt q qu qu n u þv u þv þu þv (5) 2 qx 2 qy qx qy qx qy !n vnþ1  vn qvn qvn 1 q2 v q2 v ¼ u v þ þ Dt qx qy Re qx2 qy2  n   Dt q qv qv Dt q qv qv n u þv u þv þu þv 2 qx qx 2 qy qy qx qy

(6)

where un+1 ¼ u(tn+1), un ¼ u(tn), vn+1 ¼ v(tn+1), vn ¼ v(tn), Dt ¼ tn+1tn. The higher-order terms appear due to time discretization using the characteristic concept, which act as a smoothing or stabilizing operator that reduces the oscillations arising from the spatial discretization of the convection terms.

where w(xxI) is a weight function of compact support (often called the domain of influence of node I) and n is the number of nodes whose support includes point x. uI(t) is the parameter associated with node I of the approximation field. The choice of the weight function is more or less arbitrary, the exponential function and spline functions are often used in practice, and in the paper the following cubic spline function is chosen: 8 2 2 3 0pzp0:5 > < 3  4z þ 4z ; wðzÞ ¼

4

3 > :

 4z þ 4z2  43z3 ; 0;

0:5ozp1

(9)

otherwise

Minimization of (8) with respect to a(x, t) then yields to the following system of linear equations for the coefficient a(x, t): AðxÞaðx; tÞ ¼ BðxÞuðtÞ

(10)

where A(x) and B(x) can be easily obtained from (7), and u(t) is the vector of nodal unknowns. If A is invertible, the coefficient a(x, t) can be expressed as aðx; tÞ ¼ A1 ðxÞBðxÞuðtÞ

(11)

Substituting the above equation back into (7) leads to h

u ðx; tÞ ¼ PT ðxÞA1 ðxÞBðxÞuðtÞ ¼ Nðx; tÞuðtÞ

(12)

where N(x,t) is the matrix of MLS shape functions. The MLS approximation is obtained by a special least-squares method, thus the function obtained by the MLS approximation is a smooth curve and it does not pass through the nodal values. Therefore, the MLS shape functions do not, in general, satisfy the Kronecker delta condition at each node, i.e. NI ðxJ ÞadIJ

(13)

Consequently, the imposition of essential boundary conditions is more complicated than that for the standard FEM. Several methods have been proposed, including Lagrange multipliers [11], penalty methods [13], coupled FEM method [12,21], and so on. The MLS approximation has two major features that make it popular: (a) the approximation field function is continuous and smooth in the entire problem domain and (b) it is capable of producing an approximation with the desired order of consistency. As mentioned above, the Galerkin spatial approximation is optimal when the characteristic form is used. Thus using the MLS shape function, the field variables can be easily written as follows: u ¼ Nu;

v ¼ Nv

And using the weighting NT in the integrated residual expression, we can obtain the following equations.

ARTICLE IN PRESS 358

X.H. Zhang et al. / Engineering Analysis with Boundary Elements 33 (2009) 356–362

For x Burgers’ equation

3. Results and discussion

Mðunþ1  un Þ ¼ Dt½Cun þ Kun  DtK s un þ f u 

(14)

3.1. Burgers’ equation in 1D

(15)

Burgers’ equation is now solved in [0,1] with homogeneous Dirichlet boundary conditions. The problem is defined by [1,2,5,6,10,23]

For y Burgers’ equation Mðvnþ1  vn Þ ¼ Dt½Cvn þ Kvn  DtK s vn þ f v 

qu qu 1 q2 u þu ¼ qt qx Re qx2

In the above equations, Z

N T N dO

M¼

(16)

O

  qN qN þv dO C¼ NT u qx qy O Z

1 K¼ Re

(17)

!

Z O

qN T qN qN T qN þ dO qx qx qy qy

(18)

uð0; tÞ ¼ uð1; tÞ ¼ 0:0 and with the following initial condition u(x,0) ¼ sin(px), 0oxo1. The analytical solution of this problem given by Cole [2,22] is P1 1 a expðn2 p2 t=ReÞn sinðnpxÞ n¼1 P n (21) uðx; tÞ ¼ 2p 2 2 Re a0 þ 1 n¼1 an expðn p t=ReÞ cosðnpxÞ where Z a0 ¼

1

0

Z





qN T qN qN u þv dO u qx qx qy O   Z 1 qN T qN qN u þv dO  v 2 O qy qx qy

1 Ks ¼  2

an ¼

Z 0

(19)

and the forcing terms fu and fv are the results of the application of Green’s lemma to the second-order derivatives of the differential equations. It can be noted that no stabilization parameters are involved in the EFCGM, which implies that the EFCGM can be widely applied. To evaluate the integrals in (16)–(19), it is necessary to define integral cells over the domain of the problem, and then the Gauss quadrature scheme is employed to perform the integrations numerically over these cells. Though the shape functions are time dependent, if the nodes are unchanged during the analysis, we can calculate shape functions and its derivatives only one time in order to save computational cost. Therefore, the EFCGM for Burgers’ equation can be summarized as in Table 1.

(20)

1

  Re exp  ½1  cosðpxÞ dx 2p   Re exp  ½1  cosðpxÞ cosðnpxÞ dx 2p

ðn ¼ 1; 2; 3; . . .Þ.

The EFCGM results for the example are presented for different values of Re ranging from 10 to 10,000. It must be noted that the results obtained from (21) may not be reliable for large values of Re.

Table 2 Comparison analytic and EFCGM results for Re ¼ 10 x

T

Numerical

Exact

0.25

0.4 0.6 0.8 1.0

0.308892 0.240747 0.195679 0.162560

0.30889 0.24074 0.19568 0.16256

0.5

0.4 0.6 0.8 1.0

0.569629 0.447214 0.359241 0.291919

0.56963 0.44721 0.35924 0.29192

Table 1 Flowchart of the EFCGM program for Burgers’ equation

0.75

1 Pre-processing phase  Define EFCGM nodes and their initial position and value  Compute shape function and their derivatives: – Loop over cells of the domain * Loop over quadrature points xQ in cell C * (a). If quadrature point is outside the physical domain, go to (d) (b). Check all nodes in cell and surrounding cells to determine n nodes xI, I ¼ 1 to n such that xQ is in their domain of influence (c). For each of the n neighbor nodes, compute NI(xQ) and its spatial derivatives (d). Continue; * End quadrature point loop

0.4 0.6 0.8 1.0

0.625446 0.487214 0.373923 0.287472

0.62544 0.48721 0.37392 0.28747

Table 3 Comparison analytic and EFCGM results for Re ¼ 100 x

T

Numerical

Exact

0.25

0.4 0.6 0.8 1.0

0.341914 0.268958 0.221491 0.188193

0.34191 0.26896 0.22148 0.18819

0.5

0.4 0.6 0.8 1.0

0.660712 0.529426 0.439142 0.374423

0.66071 0.52942 0.43914 0.37442

0.75

0.4 0.6 0.8 1.0

0.910256 0.767247 0.647438 0.556051

0.91026 0.76724 0.64740 0.55605

– End cell loop 2 Computational phase  Loop over time steps – Compute matrix M,C,K and Ks at each quadrature point – Compute forces vector fu and fv along the Neumann boundary – Use penalty method imposed essential boundary – Update u and v

 End loop over until the given time steps

ARTICLE IN PRESS X.H. Zhang et al. / Engineering Analysis with Boundary Elements 33 (2009) 356–362

As a matter of fact, the analytical solution of the equation involves Fourier series solutions, and the convergence is very slow due to series solutions for the larger value of Re [2]. Therefore, the numerical solutions are compared with the results obtained from (21) only for the cases of Re ¼ 10 and 100. For the case of Re ¼ 10,000, we compare the EFCGM results with some results available in the literatures [2,10,23].

Table 4 Comparison between the numerical solutions for Re ¼ 10,000 x

Ref. [23]

Ref. [2]

Ref. [10]

EFCGM

0.05 0.11 0.16 0.22 0.27 0.33 0.38 0.44 0.50 0.55 0.61 0.66 0.72 0.77 0.83 0.88 0.94 0.96 0.98 0.99

0.0424 0.0843 0.1263 0.1684 0.2103 0.2522 0.2939 0.3355 0.3769 0.4182 0.4592 0.4999 0.5404 0.5805 0.62.0 0.6600 0.6957 – – –

0.0379 0.0834 0.1213 0.1667 0.2044 0.2469 0.2872 0.322 0.3769 0.4140 0.4584 0.4951 0.5388 0.5749 0.6179 0.6533 0.6952 – – –

0.0379 0.0834 0.1213 0.1667 0.2044 0.2497 0.2872 0.3322 0.3769 0.4141 0.4584 0.4951 0.5388 0.5749 0.6179 0.6533 0.6952 0.7090 0.7228 0.7296

0.037923 0.083419 0.121288 0.166682 0.204445 0.249661 0.287239 0.332186 0.376943 0.414070 0.458385 0.495088 0.538817 0.574947 0.617909 0.653303 0.695216 0.709033 0.722767 0.729600

359

Tables 2 and 3 present the comparison between the EFCGM solutions and the analytical values for Re ¼ 10 and 100, respectively. In these cases, 20 and 50 nodes are distributed uniformly in [0,1] for Re ¼ 10 and 100, respectively. The dimension of the influence domain for each node is assumed to be 1.3Dx, where Dx is the distance between two adjacent nodes. Meanwhile, in each integral cell, single Gauss quadrature point is used and Dt ¼ 0.001. It can be noted that our numerical results are identical with the analytical solutions that are obtained from (21). For the case of Re ¼ 10,000, 500 nodes are distributed uniformly in the domain and Dt ¼ 0.0001, the other computational parameters are the same as that in the cases of Re ¼ 10 and 100. Table 4 presents the EFCGM results and those available in the literatures [2,10,23]. It can be found that our results are in agreement with the results presented in Refs. [2,10]. From these results, we conclude that the EFCGM gives remarkable accuracy for moderate values of Re and large values of Re. Fig. 1 illustrates the numerical results at different time for Re ¼ 10, 100, 1000 and 10,000, respectively. In order to demonstrate the stability of the EFCGM, the profiles of u at times 0.2, 0.4, 0.6, 0.8 and 1.0 are shown. It can be observed that for tX0.4 with large value of Re (e.g. ReX1000), the profiles of u develop extremely sharp gradients near the right-hand boundary, but the numerical results of the proposed method are stable and no spurious oscillations appear.

3.2. 2D Burgers’ equation A two-dimensional Burgers’ problem is now considered over the square domain O ¼ [0,1]  [0,1]. The problem is defined

Fig. 1. The profiles of u at different times for Re ¼ 10, 100, 1000 and 10,000. (a) Re ¼ 10, (b) Re ¼ 100, (c) Re ¼ 1000 and (d) Re ¼ 10,000.

ARTICLE IN PRESS 360

X.H. Zhang et al. / Engineering Analysis with Boundary Elements 33 (2009) 356–362

Fig. 2. The profiles of u velocity at different times for Re ¼ 100.

by coupled equations [17]

qu qu qu 1 q2 u q2 u þu þv ¼ þ qt qx qy Re qx2 qy2

!

!

qv qv qv 1 q2 v q2 v þu þv ¼ þ qt qx qy Re qx2 qy2

and by the following boundary and initial conditions: u(0,y,t) ¼ u(1,y,t) ¼ 0.0 v(x,0,t) ¼ v(x,1,t) ¼ 0.0 qu/qn(x,0,t) ¼ qu/qn(x,1,t) ¼ 0.0 qv/qn(0,y,t) ¼ qv/qn(1,y,t) ¼ 0.0 u(x,y,0) ¼ sin(px)cos(py) v(x,y,0) ¼ cos(px)sin(py) The above problem exhibits various symmetries. For example, there exist u(x,y,t) ¼ v(y,x,t) and u(x,y,t) ¼ u(1x,1y,t) over the

global domain, and u(x,x,t) ¼ v(x,x,t) and u(x,x,t) ¼ u(1x,1x,t) along the transverse section of the domain. When convection dominates the nonlinear transport, the solutions include the formation of a discontinuity along the diagonal of the domain passing through the points (0,1) and (1,0). The EFCGM results obtained on a uniform distribution of 41  41 nodes and 3  3 Gaussian quadrature points in each integral cell for Re ¼ 100 and 1000 are shown in Figs. 2 and 3, respectively. It can be found that for tX0.4, the solutions appear discontinuous lines along the domain diagonal in both Figs. 2 and 3. However, in Fig. 3 a few insignificant overshooting and undershooting are observed near the discontinuous lines for tX0.6. But compared with the results via the nonstandard FEM [17], the results by the proposed method are fairly accurate. The results of u along the domain diagonal (y ¼ x) are presented in Fig. 4 for Re ¼ 100 and 1000, respectively. It can be seen that the solutions of u(x,x,t) appear as a sharp front near the point (0.5,0.5) at about t ¼ 0.4 and afterwards the amplitude of the sharp front starts to decay. Also, the results show that the peak of shock waves for Re ¼ 1000 is sharper than that for Re ¼ 100.

ARTICLE IN PRESS X.H. Zhang et al. / Engineering Analysis with Boundary Elements 33 (2009) 356–362

Fig. 3. The profiles of u velocity at different times for Re ¼ 1000.

Fig. 4. Solutions u(x,x,t) of the EFCGM at different times for Re ¼ 100 and 1000, respectively. (a) Re ¼ 100 and (b) Re ¼ 1000.

361

ARTICLE IN PRESS 362

X.H. Zhang et al. / Engineering Analysis with Boundary Elements 33 (2009) 356–362

4. Conclusions A numerical method, which is based on characteristic and element-free Galerkin method, is developed for solving Burgers’ equation. The above procedure shows that: (1) the presented EFCGM provides a very efficient and accurate method to solve Burgers’ equation for wide range of Re, (2) based on CG method, the EFCGM is optimal for transient Burgers’ equation and no stabilization parameters need to be chosen, (3) the proposed method is fully explicit, thus it is very simple to implement and parallelize and (4) the EFCGM is a meshfree technique in which no mesh generation and mesh recreation are involved. Since there is no explicit mesh, one can easily add or delete nodes in the desired regions. Therefore, the EFCGM is an attractive approach to deal with nonlinear and high gradient problems. It is believed that this simple approach will also be useful for solving more general problems in fluid dynamics.

Acknowledgments The support from the National Natural Science Foundation of China (NSFC) (No: 10590353) and Natural Science Foundation in Shaanxi province (No: 2005A16) are fully acknowledged. References [1] Gu¨lsu M. A finite difference approach for solution of Burgers’ equation. Appl Math Comput 2006;175:1245–55. [2] Hassanien IA, Salama AA, Hosham HA. Fourth-order finite difference method for solving Burgers’ equation. Appl Math Comput 2005;170:781–800. [3] Radwan SF. Comparison of higher-order accurate schemes for solving the two-dimensional unsteady Burgers’ equation. J Comput Appl Math 2005;174: 383–97. [4] Kadalbajoo MK, Awasthi A. A numerical method based on Crank–Nicolson scheme for Burgers’ equation. Appl Math Comput 2006;182:1430–42.

¨ zis- T, Aksan EN, O ¨ zdes- A. A finite element approach for solution of Burgers’ [5] O equation. Appl Math Comput 2003;139:417–28. [6] Dogan A. A Galerkin finite element approach to Burgers’ equation. Appl Math Comput 2004;157:331–46. [7] Aksan EN. A numerical solution of Burgers’ equation by finite element method constructed on the method of discretization in time. Appl Math Comput 2005;170:895–904. [8] Chino E, Tosaka N. Dual reciprocity boundary element analysis of timeindependent Burgers’ equation. Eng Anal Boundary Elem 1998;21:261–70. [9] EI-Hawary HM, Abdel-Rahman EO. Numerical solution of the generalized Burger’s equation via spectral/spline methods. Appl Math Comput 2005; 170:267–79. [10] Hashemian A, Shodja HM. A meshless approach for solution of Burgers’ equation. J Comput Appl Math 2007. [11] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Eng 1994;37(2):229–56. [12] Belytschko T, Krongauz Y, Organ D, et al. Meshless method: an overview and recent developments. Comput Methods Appl Mech Eng 1996;139(1–4): 3–47. [13] Liu GR, Gu YT. An introduction to meshfree methods and their programming. Berlin and New York: Springer; 2005. [14] Ouyang J, Zhang L, Zhang XH. Nonstandard element-free Galerkin method for solving unsteady convection dominated problems. Acta Aerodyn Sin 2007;25:287–93. [15] Xie SS, Heo S, Kim S, et al. Numerical solution of one-dimensional Burgers’ equation using reproducing kernel function. J Comput Appl Math 2008; 214(2):417–34. [16] Young DL, Fan CM, Hu SP, Atluri SN. The Eulerian–Lagrangian method of fundamental solutions for two-dimensional unsteady Burgers’ equations. Eng Anal Boundary Elem 2008;32(5):395–412. [17] Donea J, Huerta A. Finite element methods for flow problems. Chichester: Wiley; 2003. [18] Zienkiewicz OC, Taylor RL. The finite element method, 5th edn, vol. 3. Butterworth-Heinemann: London, 2000. [19] Lewis RW, Nithiarasu P, Seetharamu KN. Fundamentals of the finite element method for heat and fluid flow. Chichester: Wiley; 2004. [20] Zienkiewicz OC, Nithiarasu P, Codina R, et al. The characteristic-based-split procedure: an efficient and accurate algorithm for fluid problems. Int J Numer Methods Fluids 1999;31:359–92. [21] Zhang XH, Ouyang J. Meshless analysis of heat transfer due to viscous dissipation in polymer flow. Eng Anal Boundary Elem 2008;32:41–51. [22] Cole JD. On a quasi linear parabolic equation occurring in aerodynamics. Q Appl Math 1951;9:225–36. [23] Hon YC, Mao XZ. An efficient numerical scheme for Burgers’ equation. Appl Math Comput 1998;95:37–50.