An Efficient Scheme of Differential Quadrature Based on Upwind Difference for Solving Twodimensional Heat Transfer Problems
AbdulSattar Jaber Ali AlSaif
Department of Mathematics, College of Education for Pure Science, University of Basrah, Basrah, Iraq
Email address:
To cite this article:
AbdulSattar Jaber Ali AlSaif. An Efficient Scheme of Differential Quadrature Based on Upwind Difference for Solving Twodimensional Heat Transfer Problems. Applied and Computational Mathematics. Vol. 4, No. 4, 2015, pp. 275285. 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 twodimensional heat transfer (convectiondiffusion) 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 convectiondiffusion problems ((i) Steadystate 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, Twodimensional Convectiondiffusion, 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 BurgerHuxley equation, which is one of the type of convectionsdiffusion 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 Gausselimination 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 nondimensional governing equation as
(1)
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 differencedifferential quadrature method for solving twodimensional 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 differencedifferential 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
(2)
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 firstorder derivative are given as
, but (3)
where ,
The weighting coefficients for the secondorder and highorder derivatives are given as
For but, (4)
where
For, (5)
The weighting coefficient of derivatives with respect tocan also be obtained in similar form like those in formulae (35) .
We choose type of sampling points from extremely useful formulas, taken from Refs. [10,19]. This kind of sampling points is the ChebyshevGaussLobatto points; forspace given as,
, for .
3. Upwind Difference Differential Quadrature Method (UDDQM)
Here, a UDDQM method will be suggested to solve the twodimensional 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 operatorsplitting 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;
(6)
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
(7)
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;
(8)
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:
(9)
(10)
From Equations (9) and (10) we can computed by the flowing form
(11)
where,
,
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 coordinates is consider as;
(12)
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.
The problem selected to test the implication of new numerical method is that of the fluid state of solidbody rotation[18] (see Figure1), the nondimensional stream function at any point in the region of solution, given by;
(13)
where and are dimensionless coordinates. The analytic solution of Equation (12) is
(14)
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
(15)
For the second case, was assumed to vary as where is constant. the parabolic corresponding exact distribution of temperature is
(16)
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 twodimensional 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 nondimensional parameter, which is called Peclet number () and written as;
(17)
the numerical results were obtained for a squareshaped 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 GaussSeidel 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:
(18)
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 25 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 27, we see that UDDQM for all values of and is better in accuracy and convergence than the other methods that are used for comparison.
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 noslip condition on boundaries.
In order to facilitate the analysis, the following assumptions and dimensionless variables are considered:
1. The fluid is Newtonian and the flow is twodimensional 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 nondimensional governing equations for vorticity (), stream function, and energy can be written as
(19a)
(19b)
(19c)
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 ChebyshevGaussLobatto 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 flowvortex 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 thermogravitational 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 flowvortex (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 thermogravitational convection is weak due to the low Grashof number and low temperature, and become stronger gradually with increases in temperature.
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 fivepoint UDDQM agreement with existing results, threepoint FDM and DQM. It pointed out the secondary –flow vortex appears starting the wall of channel toward the whole domain consequently.
Table1 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 higherorder 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 21522157).
Table 1. Comparison of the present study with [2,21], for , .
Method  Grids 



QUICK DQM  81x21
 0.4579 0.4552  0.6882 0.5917  0 .7951 0.6904 
UDDQM 
 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 firstorder derivative () and the secondorder 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 firstorder 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 loworder 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 twodimensional convection and diffusion equation. In addition, the effects of a heated on natural convection of fluid flow in solidbody 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 loworder finite difference method, in which a large number of grid points must usually be used.
References