Applied and Computational Mathematics
Volume 4, Issue 4, August 2015, Pages: 275-285

An Efficient Scheme of Differential Quadrature Based on Upwind Difference for Solving Two-dimensional Heat Transfer Problems

Abdul-Sattar Jaber Ali Al-Saif

Department of Mathematics, College of Education for Pure Science, University of Basrah, Basrah, Iraq

Email address:

To cite this article:

Abdul-Sattar Jaber Ali Al-Saif. An Efficient Scheme of Differential Quadrature Based on Upwind Difference for Solving Two-dimensional Heat Transfer Problems. Applied and Computational Mathematics. Vol. 4, No. 4, 2015, pp. 275-285. doi: 10.11648/j.acm.20150404.16

Abstract: In this paper, a new technique of differential quadature method called the upwind difference - differential quadature method (UDDQM) for solving two-dimensional heat transfer (convection-diffusion) problems is proposed. Also, investigated the effects of physical quantities on behavior of flow problems, and combined effects of upwind difference mechanism together with differential quadrature method to modified the numerical solutions of heat transfer problems are presented. To validate our proposed UDDQM, two convection-diffusion problems ((i) Steady-state incompressible flow problem has exact solution and (ii) Natural convection motion of the incompressible fluid flow problem hasn't exact solution) are solving numerically. Graphical results on the effects of parameter variation on velocity, temperature, Peclet number, Grashof number, and Prandtl number are presented and discussed. Numerical experiments are conducted to test its accuracy and convergence and compare it with the standard DQM and other numerical methods that are available in literature. The numerical results show the efficiency of the proposed method to handle the problems, and it is more accurate and convergent than other methods.

Keywords: Upwind Difference Scheme, Differential Quadrature Method, Two-dimensional Convection-diffusion, Accuracy, Convergence

1. Introduction

The heat transfer (convection and diffusion) problems attract many researchers interest and they have wide applications, for example, in energy system, which includes the solar collector, nuclear reactor, heat transfer, natural convective motion of fluid flow…etc. Most numerical simulations for these problems can be currently carried out by using finite differences method(FDM) and finite elements method(FEM) [1,3,6,7,10,11,14,16,17,18,22]. These techniques for solving numerically these problems may be very complex and require a large number of grid points to obtain accurate solutions. The differential quadrature method which is introduced by Bellman et al.[4], is able to overcome these difficulties with few grid points and reasonable computational workloads. This fact is mentioned in many articles [2,4,10,11,23,24]. Recently, using differential quadrature procedures to solve one- and two- dimensional differential equations are given accurate results using less discretization points frequently, for example; Izadian et al.[8] they studied Burger-Huxley equation, which is one of the type of convections-diffusion equation, by two important numerical methods, spectral method and DQM. The result show that the spectral method is relatively better than DQM.

Jiwari et al.[9] have applied a weighted average of DQM on one dimensional Burgers’ equation. The authors reduced the equations into a system of linear equations by DQM and then solved by Gauss-elimination method. Marji and his associate’s [12] they improved the application of DQM for the solution of Navier–Stokes equations based on the point pressure–velocity iteration method. In [13] Meral was used DQM to discretize one- and two- dimensional nonlinear heat- and mass- transfer equations and using the Runge–Kutta method to solved the resulting nonlinear system of ordinary differential equations numerically. Öğüt [17], investigate the effects of the direction of magnetic force, heater length, and heater location on heat and fluid flow in the enclosure, and he showed that isotherms and streamlines are employed generally to depict the mechanism of heat and fluid flow, and the solution of the governing equations by using the differential quadrature method is an effective technique. Verma et al.[25] have used DQM to obtain the numerical solutions for nonlinear diffusion equations. In addition to the above mentioned facts, the upwind mechanism is important for the computation of fluid flow problems, and that mechanism is absent(or used little) in DQM. According to these reasons, the solution of flow problems by DQ motivates the present work.

In this paper, the system of convection and diffusion problem can be described by the entire system of non-dimensional governing equation as


where, is interpreted as column vector of unknown functions, is the column vector of the velocity of flow,  is the column vector of the variable diffusivity of diffusion terms, and is the column vector of potential functions, with appropriate initial and boundary conditions.

In recent years, many models have been used to describe flow problems that involve two mechanism of the flow (convection and diffusion)[8,17,20]. In this paper, we will develop the differential quadrature method to solve heat transfer problems. Upwind difference- differential quadrature method is presented by using the traditional upwind difference scheme for the convective terms and the traditional DQ formulas for diffusive terms in the flow models that under consideration. Using the upwind difference-differential quadrature method for solving two-dimensional convection and diffusion problems excellent numerical results are obtained. Compare with other methods; the new method with a few grid points appears that it has better convergence and accuracy than the other methods.

This paper is organized as follows: the ideas of differential quadrature method are illustrated in section 2. Upwind difference-differential quadrature method is illustrated in section 3. To demonstrated the efficiency of new method two test problems are presented in section 4. In section 5 analysis of errors and stability are introduced. Finally, the conclusions are reported.

2. Differential Quadrature Formulations

The essence of the DQM is that the partial derivatives of a function with respect to a variable (ordirection) in a governing equation are approximated by a weighted linear sum of function values at all discrete points. We consider a function defined on a rectangular domain , with and  fixed. Suppose that, the grid is obtained by taking  and points in the - and- directions, respectively. Then, the –order -partial derivative () of the function at a point, and  –order -partial derivative () of the function at a point, respectively may be approximately written as


where ,  ,  , , and  are the respective weighting coefficients. The determinant of the weighting coefficients and the choice of sampling points are very important factors for the accuracy of the DQ method. The weighting coefficients for the derivatives may be obtained directly, and most accurately, irrespective of the number and positions of the sampling points. From Shu’s [19] the weighting coefficients of DQ method are as follows:

The weighting coefficients for the first-order derivative are given as

 ,  but       (3)

where  ,

The weighting coefficients for the second-order and high-order derivatives are given as

 For  but,                                          (4)


 For,        (5)

The weighting coefficient of derivatives with respect tocan also be obtained in similar form like those in formulae (3-5) .

We choose type of sampling points from extremely useful formulas, taken from Refs. [10,19]. This kind of sampling points is the Chebyshev-Gauss-Lobatto points; forspace given as,

 , for  .

3. Upwind Difference- Differential Quadrature Method (UDDQM)

Here, a UDDQM method will be suggested to solve the two-dimensional incompressible flow problems. This method includes two points: First, the upwind difference scheme is used to approximate the convective terms subject to mechanisms of the flow directions, and DQM for other terms in space. The second two dependent stages are adopted to calculate the results at every step, similar to the operator-splitting algorithm employed by [15] to solve heat transport problem. In order to demonstrate these ideas: we replace a convective terms by the upwind difference scheme, after the average axial evaluated at half a grid forward and backward from (,) in the - and - directions respectively. The ideas of upwind difference scheme are used to approximate the convective terms extensively in many texts (for example see[5]). It can easily be verified that the upwind difference form is automatically preserved when the following numerical formulas are used;


where, and  and  are forward and backward difference respectively.

The remaining terms are approximated by DQM in spatial and forward difference scheme in temporal. By these methods mentioned above the Equation (1) is discretized into following equation


On the first stage at every step,  is obtained by (7). In the second stage at every step, the discrete system is obtained by traditional DQM and the  is taken as predictions of. The formulae for the second computation stage at every step are written as;


where ,s should be all internal grid points

The can be obtained by solving approximation (8). Generally, if we assume that

Equations (7) and (8) can be written in matrix form, respectively as:



From Equations (9) and (10) we can computed  by the flowing form




and are  unsymmetrical matrices,

The subscription  and are defined by  in the following form and ,for and .

 are the Kroncker deltas.

To show the effectiveness of the proposed method and to give a clear overview of the methodology, two examples of the convection and diffusion flow problem are going to be discuss in the following section.

4. Applications

We shall apply the upwind difference- differential quadrature method on two flow problems. Nonetheless, such an approach is needed to evaluate the accuracy and convergence of the new numerical method. All the results are calculated by using PC of kind Pentium 4, and the programs are written by Fortran language. These applications are given as in the illustrative examples.

Example 1

For a two –dimensional, steady state, incompressible flow, the heat transfer equation in Cartesian co-ordinates is consider as;


where is the temperature of fluid, and  are variable of mass velocity of convection terms and related with stream function by  and , and  variable diffusivity of diffusion terms. Equation (12) reduced from Equation (1) by setting, and.The recurrence relations (9),(10),and (11) are obtained by applying the new method into equation (12) with appropriate boundary conditions.

Figure 1. Control volume of the flow.

The problem selected to test the implication of new numerical method is that of the fluid state of solid-body rotation[18] (see Figure1), the non-dimensional stream function at any point in the region of solution, given by;


where  and are dimensionless coordinates. The analytic solution of Equation (12) is


where  and are dimensionless quantities. The numerical solution of Equation (12) can be obtained by two different distribution of. In the first case: , where is constant. The logarithmic corresponding exact distribution of temperature is


For the second case,  was assumed to vary as  where is constant. the parabolic corresponding exact distribution of temperature is


where,  are constants.

For the simple test and analysis of the method, were noted that none of convection terms and diffusion terms are identically zero for all grid nodes. Since the numerical method doses not know of existence of an analytical solution, it goes on to solve it in its own way as a two-dimensional problem. Moreover, the influence of the convective terms is absent in the analytical solution (14), it inters into the differential quadrature calculations. The behaviour of the numerical solution is determined by the relative value of convective terms and diffusive terms and can be characterized in terms of a non-dimensional parameter, which is called Peclet number () and written as;


the numerical results were obtained for a square-shaped control volume(Figure 1), such that, side of square equal to units, the corner nearest to the origin situated at  , and a square's covered by grid with, and initial condition at all interior nodes was put zero. The boundary condition is useful to calculate the constants in temperature Equations (15) and (16), where the first value of equal to a minimum value of zero at outermost corner and the other equal to maximum value of unity at the innermost corner. The numerical solution in both cases results from combining the forces of diffusion and convection, which the values of  will be in it, in which iteration Gauss-Seidel iterative method is used to find.

The numerical solution for the equations system, which results from using UDDQM to approximate Equation (12) with initial, and boundary conditions is obtained. The points that are used to divide the solution region were  in ordirections. The results which we are going to discuss later will be related withand, which its effect will be clear through the numerical results and the conclusions (i.e. andrelated with Peclet number (17)). To illustrate that, the definition of for logarithm distributing of temperature, Equation (15), by using Equation(13) and  in the form:


Largest value to  finds at the largest value of ,and the largest value in any internal point  will be, for example; when result, in the same analysis yield for the parabolic distributing of temperature, Equation(16), there are many schemes that are unstable numerically when [18,22]. Figures 2-5 are shown that the curves of iteration number (IN) and the relative errors (Re) of the numerical results obtained by used DQM and a new technique UDDQM. From these figures, we notice that for the two distributions of temperature the IN increases with increases N, and the numerical results are convergence and stable at the lowerest value of and fail for  and approximately, i.e. Also, Re decreases with decrease in , but Re increases with increase in , for parabolic distribution. From these results, we can say that UDDQM gives high accuracy with few points. UDDQM always convergence and gives accurate results than DQM for all sample of grid points and for all values of and. Here, to prove the efficiency of UDDQM in accuracy and speed of convergence, we compare it with other numerical methods at , such as central difference method(CDM),upwind difference method(UDM),and Spalding difference method(SDM),for more details look at [18,22]. From the results are representing in the Figures 2-7, we see that UDDQM for all values of and is better in accuracy and convergence than the other methods that are used for comparison.

Figure 2. Convergence of DQM –logarithmic and parabolic distribution temperature.

Figure 3. Convergence of UDDQM –logarithmic and parabolic distribution temperature.

Figure 4. Relative errors of DQM –logarithmic and parabolic distribution temperature.

Figure 5. Relative errors of UDDQM –logarithmic and parabolic distribution temperature.

Figure 6. Comparison iteration numbers between UDDQM and others numerical methods-logarithmic and parabolic temperature distribution.

Figure 7. Comparison relative errors between UDDQM and others numerical methods-logarithmic and parabolic temperature distribution.

Example 2

We consider liquid enclosed by a rectangular box of highest, length, shown in Figure 8. The fluid is heated from below and cooled at the upper wall. At, the upper wall is of a constant temperature  and  on the lower plate, as time increasers fluid temperature changes but the values at solid walls are always kept at the initial condition, the velocities are satisfying no-slip condition on boundaries.

Figure 8. Schematic diagram of the test region with boundary conditions.

In order to facilitate the analysis, the following assumptions and dimensionless variables are considered:

1.   The fluid is Newtonian and the flow is two-dimensional laminar.

2.   The fluid is assumed to be incompressible[7]

3.   The Boussinesq approximation is valid.



where,  are the references of length, velocity and temperature difference respectively. Using above assumptions and dimensionless variables on the Navier –Stokes equations [5], the non-dimensional governing equations for vorticity (), stream function, and energy can be written as




 and                          (19d)

In which  the velocity components are inanddirections respectively,  is the Grashof number, and is the Prandtl number. Equations(19 a,c) can be reduced from(1) by setting , ,  ,and  ,and using the above assumptions and dimensionless variables.

Boundary conditions:

Because the symmetry of the temperature field , the boundary condition along the axis is, on the three solid walls comes from the fact that there is no net flow across those boundaries. Also, because the symmetry about the axis, it is necessary to seek a solution for the left of the flow. The boundary conditions along the axes are:

, , , , , ,

, , , , , ,

,,, , , ,

, ,, , , ,

The computationally efficiency of UDDQM with Chebyshev-Gauss-Lobatto points [2,10,19] on the numerical accurate results have been well demonstrated here. In the present computations, we adopted Gauss seidel method to solve vorticity and energy equations (Equations(19a,c)), and SOR with damping factor, to solve stream function(Equation(19b)). The sufficient condition for convergence of numerical solution of UDDQM is Max,(where,is the column vector of unknown variables). If it not satisfied, the maximum iterations can be determined by trail to stop the iterative procedure and the steady state solutions of the incompressible flow obtained. The streamlines in three dimensions are shown in Figure 9, at prime time the secondary flow-vortex presented at a corner where the cold and hot walls intersect. In the remaining region flow, the fluid is practically motionless. In this moment, the stream function value is positive () everywhere, that is; the fluid motion is in the counterclockwise direction. Conversely, the motion is clockwise along closed streamlines of negative () values for the stream function. Thus, the fluid descends along the cold vertical wall and then rises after flowing over the hot surface. At this time, the convective motion appears weakly. This motion becomes stronger gradually with the increase in time; also, the isotherm profiles have a slow variation with respect to time. Although a steady state has not been reached yet, the flow does not seem to have any further significant changes. Thus, additional computations with maximum step greater than 2000 have not been attempted. From Figure 10, we see that the thermo-gravitational convection is weak due to the low Grashof number. As a result, the values of  at are positive () and atare () . The isotherms gather densely near the horizontal walls and the fluid flows mainly in the clockwise direction like a rotating flow around the core region. The regions of the isotherms gathering densely are in a long region of the horizontal walls. It is observed that the fluid motion at is not steady but oscillating in nature, this is the same phenomena with[2].

Figure 9. Surface with contour of streamlines for , .

Figure 10. Streamlines and isotherm for different values of Grashof number with .

We shall discuss the effect of temperature on the behaviors of fluid motion in aspect of secondary flow-vortex (which is representing by sequential plots of the streamlines). Here, an important matter is to specify the critical values, which are indicated to the changing of the streamlines or the total number of the convection cells that appear in the channel. Figures 11, shows the temperature effect on the behaviors of motion of the fluid for, in the channel. From these figures, we notice that, there is one cell at, the stream function value is positive (+ve) everywhere, that is; the fluid motion is in the counterclockwise direction. The numbers of cells are changing into two cells on. Moreover, when the two cells are different in size and motion direction (the left cell is counterclockwise direction and the right cells is clockwise direction). The size of the right cell is changing onclearly. This case goes on for. In the remaining region flow, the fluid is practically motionless. Conversely, the motion is clockwise along closed streamlines of negative () values for the stream function. Thus, the fluid descends along the cold vertical wall and then rises after flowing over the hot surface. In this case, the convective motion appears weakly. This motion becomes stronger gradually with the increase in temperature; also, the isotherm profiles have a slow variation (identical with those in Figure 10). We conclude that the total number of the cells (secondary flow vortex) is changing with respect to temperature. That is, total number of cells is increased with increase in temperature. The same phenomena happen for the streamlines of cells. We see that the thermo-gravitational convection is weak due to the low Grashof number and low temperature, and become stronger gradually with increases in temperature.

Figure 11. Effect of the temperature of the shape on the streamlines.

The change in the number of cells of the secondary flow continues with increase of temperature values, Another change occurs in the number of cells at , where the total number of cells become three in the channel. These cells are different in size and direction (the left +ve , the middle is –ve , the right is +ve). The size of the right cell starts to change until it becomes identical with other cells approximately.

Figure 12. Surface of in 3D with contour lines for UFDM & UDQM at .

From the results are shown in Figure 12, we see that the results are obtained by five-point UDDQM agreement with existing results, three-point FDM and DQM. It pointed out the secondary –flow vortex appears starting the wall of channel toward the whole domain consequently.

Table-1 shows that comparison of the solver method and numerical results of the present work with these are presented by DQM and QUICK scheme, which stands for Quadratic Upstream Interpolation for Convective Kinematics, is a higher-order differencing scheme that considers a three

point upstream weighted quadratic interpolation for the cell phase values[21]. It is most noteworthy that the present computation is very close to DQM and QUICK. The present results are in agreement with those given by [14] in the qualitative analysis for the distribution of streamlines and isotherms (see pages 2152-2157).

Table 1. Comparison of the present study with [2,21], for , .

Method Grids



0.4579 0.4552 0.6882 0.5917 0 .7951 0.6904

0.4573 0.6776 0.7943

5. Error Analysis

Analyzing the errors resulting from approximation of a function and derivatives is useful work. Depending on the DQ is identical to Lagrange polynomial interpolation of order , author [19] have given a thorough error analysis for the first-order derivative () and the second-order derivative (). These errors written as;

These residual estimates show that very high accuracy can be obtained if the number of grid points  is large. Accuracy is proportional to  or its powers. By "its powers", we mean here that accuracy may also be proportional to squared or cubic, or even higher order term of. However, too large may lead to instability this shown in [24]. For dynamical problems, small error may accumulate with each time step. A simple estimate for the first-order derivative is then that the accumulated error at each time is proportional to , where is time step. Similar arguments lead us to the conclusions that the accumulated error at each time step for the second order derivative is proportional to. If largeis used, time step must be small to keep the errors within controllable range. High order differential quadrature discretization becomes unstable faster than low-order DQ discretization. This simple analysis and the numerical results lead us to the conclusion that DQ is of high accuracy, but of poor stability. The more grid points are used, the high accuracy we obtained, but the poorer the stability is. Stability of the function values and its derivative Lagrange polynomial interpolation is a complicated problem. From this rude estimate, however, we conclude that accuracy and stability are conflicting. Accuracy requires large number of grid points, but stability requires the opposite.

6. Conclusions

The new method UDDQM is applied successfully to solve the steady state two-dimensional convection and diffusion equation. In addition, the effects of a heated on natural convection of fluid flow in solid-body rotating (Figure 1) and rectangular box (Figure 8) are examined. The results of a new version of differential quadrature technique are compared with finite difference techniques and other numerical techniques [2],[18] and [21]. The numerical results obtained show that the UDDQM owns advantages including the higher accuracy and convergence by using few grid points, compared with conventional low-order finite difference method, in which a large number of grid points must usually be used.


  1. Al-Saif A.S.J. and Al-HumediM.O " Aweightedfinitedifferencemethodinvolvingnine-pointformula for two-dimensional convection-diffusion equation", Int. Math. Forum, 5, 2010, 3379 – 3398.
  2. Al-Saif A.S.J." Numerical study for convection motion stability of the incompressible two- dimensional fluid flow by DQM", Journal Basrsh Res.Sci.,33(1),2007,103-110.
  3. Anderson J.D. " Computational Fluid Dynamics" The basic with application, McGraw-Hill Companies, Inc.(1995).
  4. Bellman R., Kashef G., and Casti J." Differential quadrature: A technique for the rapid solution of nonlinear partial differential equations" J. Comput. Phys , 10 (1),1972, 40- 52.
  5. Biringen S. and Chow C.Y. "An Introduction to Computational Fluid Mechanics by example" John Wiley &Sons, Inc.,2011
  6. Dehghan M. and Mohebbi A., High-order compact boundary value method for the solution of unsteady convection–diffusion problems, mathematics and computer in simulation,79,2008,683-699.
  7. Fu W.S. and Shieh W.J.,"A study of thermal convection in an enclosure induced simultaneously by gravity and vibration",Int. J. of Heat and Mass Transfers, 35(7) , 1992,1695-1710.
  8. Izadian J., Nateghi F., and Jalili M."Comparison of spectral and differential quadrature methods for solving the Burger-Huxley equation" Communications in Numerical Analysis, 2013, 2013, 1-7
  9. Jiwari R.,Mittal R.C.,and Sharma K.K." A numerical scheme based on weighted average differential quadrature method for the numerical solution of Burgers’ equation" Appl. Math . Comput. , 219, 2013, 6680-6691.
  10. Kajal V." Numerical solutions of some reaction-diffusion equation by differential quadrature method" Int. J. Appl. Math. And Mech., 6(14), 2010,68-80.
  11. Kaya B." Solution of the advection-diffusion equation using the differential quadratur method" KSCE J.Civil Engng.,14(1), 2010,69-75.
  12. Meraji S. H.,Ghaheri A, Malekzadeh P. "An efficient algorithm based on the differential quadrature method for solving Navier–Stokes equations", Int. J. Numer. Meth. Fluids; 71,2013,422–445.
  13. Meral G. " Differential quadrature solution of heat- and mass-transfer equations " Appl. Math. Model. , 37, 2013, 4350-4359.
  14. Mukutmoni D.and Yang K.t." Wave number selection for rayleigh-Benard convection in a small aspect ratio box", Int. J. Heat Mass.Transfer.35(9),1992, 2145-2159.
  15. Muralidhar K., Verghese M., and Pilli K.M." Application operator splitting algorithm for advection-diffusion problem" Numer. Heat Trans. A, 23, 1993; 99-113
  16. Osthuizen PH.,and Panl JT., "Unsteady natural convective flow in an enclosure with a periodically varying side wall heat flux ", Adv. Comput. Meth. Heat TransferIV, Udine, Italy,1996, 53-62.
  17. Öğüt E. B." Magnetohydrodynamic natural convection flow in an enclosure with a finite length heater using the differential quadrature (DQ) method" Num. Heat Trans., Part A, ,58,2010, 900–921.
  18. Runchal A.K.."Convergence and accuracy of three finite difference schemes for a two- dimensional conduction-convection problem", Int. J. Num. Math. Engng. ,4,1972,541-550.
  19. Shu C."Differential Quadrature and its Application in Engineering" London(2000).
  20. Skiba Y.N." A non-iterative implicit algorithm for the solution of advection–diffusion equation on a sphere " Int. J. Numer. Meth. Fluids (2015) Published online in Wiley Online Library ( DOI: 10.1002/fld.4016
  21. Soong C.Y.,Tzeng P.Y.,and Hsieh C.D, "Numerical study of bottom-wall temperature modulation effects on thermal instability and oscillatory cellular convection in a rectangular enclosure ", Int. J. Heat Mass Transfers, 2001,44(20),3855-3868.
  22. Spalding D. R."A novel finite difference formulation for differential expressions involving both first and second derivatives", Int. J. Num. Math. Engng. ,4, 1972,551-559.
  23. Zekeriya G."Solution of the Blasius and Sakiadis equation by generalized iterative differential quadrature method", Int. J. Numer. Meth. Biomed. Engng. 2011; 27:1225–1234.
  24. Zong Z., and Lam K.Y." A localized differential quadrature method and its application to the 2D wave equation ", Comput. Mech., 29,2002, 382-391.
  25. Verma A., Jiwari R. and Koksal M.E." Analytic and numerical solutions of nonlinear diffusion equations via symmetry reductions"Advances in Difference Equations 2014, 2014:229

Article Tools
Follow on us
Science Publishing Group
NEW YORK, NY 10018
Tel: (001)347-688-8931