Engineering, Technology & Applied Science Research Vol. 8, No. 2, 2018, 2680-2684 2680 www.etasr.com Abdullah: Formulation of Low Peclet Number Based Grid Expansion Factor for the Solution … Formulation of Low Peclet Number Based Grid Expansion Factor for the Solution of the Convection-Diffusion Equation Aslam Abdullah Department of Aeronautical Engineering Universiti Tun Hussein Onn Malaysia Batu Pahat, Malaysia aslam@uthm.edu.my Abstract—Convection-diffusion problems, due to its fundamental nature, are found in various science and engineering applications. In this research, the importance of the relationship between grid structure and flow parameters in such problems is emphasized. In particular, we propose a systematic technique in the selection of the grid expansion factor based on its logarithmic relationship with low Peclet number. Such linear mathematical connection between the two non-dimensional parameters serves as a guideline for more structured decision-making and improves the heuristic process in the determination of the computational domain grid for the numerical solution of convection-diffusion equations especially in the prediction of the concentration of the scalar. Results confirm the effectiveness of the new approach. Keywords-convection-diffusion equations; finite difference method; non-uniform grid; grid-expansion factor; Thomas algorithm I. INTRODUCTION In Cartesian coordinates and tensor notation following the Einstein convention, the generic conservation equation in the partial differential form is ( ) + − − = 0  where is the density, is the conserved property, are velocity components of the fluid in the axes directions at the point ( , , ) at time , is the diffusivity of , and is the source or sink of . Navier-Stokes equations, which have special features of mass and momentum conservations, are extensions of this equation. The convection-diffusion equation (CDE) takes the simplified form of (1)  ( ) − = 0  where zero source/sink is assumed. The first term in (2) is called the substantial derivative  ( ) = ( ) +   The substantial derivative ( ) is physically interpreted as the time rate of change in ( ) following a moving fluid element. The first and second terms on the RHS of (3) are called the local derivative ( ) (i.e. the physical change in ( ) with time at a fixed position), and the convective derivative (i.e. the physical change in ( ) with time due to the mass transfer and change in its properties from one spatial position to another), respectively. Substituting (3) into (2) we have  ( ) + − = 0  In case of fluids at rest, or of small velocity ( ≈ 0), or large diffusivity , as well as solids, (4) is further simplified into  ( ) − = 0, representing the pure diffusion process where the local derivative ( ) is balanced by the diffusive derivative . In this paper, we consider the steady one-dimensional convection-diffusion problem where (4) reduces to  ( ) − ( ) = 0,  involving the scalar whose concentration is denoted by . Such scalar is carried along with the moving fluid element (convection) and spreads due to diffusion. Given appropriate boundary conditions, it can be shown that at relatively high velocity , or low diffusivity , the scalar concentration initially grows slowly in space and then suddenly rises over a defined distance. The sudden growth of not only provides a severe test of the discretization method, but also in the selection of compatible grid structure over the computation domain. Note that the solution of (5) is linear in space when is negligible. We investigate the relationship between the flow parameter of interest (i.e. the Peclet number ) in CDE and the expansion factor based grid structure in finite difference Engineering, Technology & Applied Science Research Vol. 8, No. 2, 2018, 2680-2684 2681 www.etasr.com Abdullah: Formulation of Low Peclet Number Based Grid Expansion Factor for the Solution … numerical scheme, and formulate the mathematical relationship between and which is necessary in achieving numerical accuracy in the solution of the equation, thus unify the deduction of heuristic selections of grid expansion factor for solving the contaminated fluids problem that leads to less pre- computation time. Note that the relatively smaller grid expansion factor does not unconditionally lead to higher numerical accuracy. II. CONVECTION DIFFUSION PROBLEMS Various numerical methods for solving CDE are by now well formulated and many useful schemes can be found such as finite differences, finite elements, spectral procedures, and the method of lines [1-12]. For instance, [1] presented a comparative study between the two most popular Lattice Boltzmann (LB) models for CDE (i.e. those in two dimensions with five and nine discrete lattice velocities, respectively). Other variants include multiple-relaxation-time LB model for the axisymmetric, as well as isotropic and anisotropic diffusion processes whose both applicability and accuracies have been investigated by [2] and [3] respectively. For the latter case, [4] proposed a finite-difference LB model for nonlinear equations. In the problem where no scalar or flux jump exists, [5] introduced a numerical scheme for dealing with curved interfaces with second-order spatial accuracy in conjunction with the LB method. Authors in [6] summarized well-known a priori error estimates for the discontinuous Galerkin approximation which carry over to the subspace of the discontinuous piecewise-quadratic space, while authors in [7] proposed the approximation of high order alternating evolution. Both [8] and [9] considered compact difference scheme for solving CDE. Authors in [8] claimed that the fourth-order scheme requires only 15 grid points, while authors in [9] successfully proved that it is computationally more efficient than the standard second-order central difference scheme. Recent methods include those that solve nonlinear fractional CDE, as homotopy analysis transform and homotopy perturbation Sumudu transform methods whose reliability and efficiency were clearly demonstrated in [10], and that based on the operational matrices of shifted Jacobi polynomials of high accuracy [11]. Author in [12] introduced a Schwarz waveform relaxation algorithm for the CDE that converges without overlap of the subdomains. The choice of suitable computational grid to discretize the governing partial differential equations (e.g. by means of polynomial fitting, Taylor series expansion and compact scheme to obtain approximations to the derivatives of the variables with respect to the coordinates) is necessary at the onset of numerical modeling of the convection-diffusion problems as in [1-17]. It is worth noting that the variable values at locations other than the defined grid nodes can also be determined by interpolation. Another important aspect is the method to solve the discretized algebraic equations. The solution is obtained via either direct [18-20] or iterative [21-24] methods. In this paper, the steady one-dimensional CDE is discretized by finite difference techniques on non-uniform grids with a defined number of nodes, and the solutions are obtained by using Thomas’ algorithm (i.e. direct method). III. DISCRETIZATION AND SOLUTION The starting point is the CDE in differential form as given by (5). Defining the boundary conditions as  (0) = 0,(1) = 1  Here we define the Peclet number as  =  The influence of the Peclet number on the diffusivity coefficient can be found in [25]. The profiles for different ranges of are illustrated in Figure 1. Fig. 1. Boundary conditions and solution profiles as a function of the Peclet number. The corresponding solution domain is covered by a grid. We define the independent variables whose domain is discretized. The interval = 0, ( − 1) is subdivided into ( − 1)/ℎ subintervals where and ℎ are odd and even integers, respectively. The nodes are defined by  = + ∆  where 1 ≤ ≤ ( − 1) , ∈ ℤ , and is the grid expansion factor. Clearly ∑ ∆ = ( − 1) . The grid is shown in Figure 2. Fig. 2. Computational molecules. At each node, the governing equation is approximated by replacing the partial derivatives with nodal values. The result is an algebraic CDE per node, in which the variables at that and immediate nodes appear as unknown. The equation system of is expressed by  + ∑ =   where P signifies the nodes at which the equations are assigned and index runs over the immediate nodes. The corresponding matrix in (7) has non-zero terms only on its main diagonal (represented by ) and the diagonals immediately above and below it (represented by and , Engineering, Technology & Applied Science Research Vol. 8, No. 2, 2018, 2680-2684 2682 www.etasr.com Abdullah: Formulation of Low Peclet Number Based Grid Expansion Factor for the Solution … respectively). The matrix elements are stored as three × array. Using the three-point computational molecules, (7) becomes  + + =   Since the convection-diffusion differential equation is linear, then the approximation contains only linear terms, and the numerical solution will not require linearization. The central difference scheme (CDS) is used to discretize the diffusion term, both for the outer derivative  − ( ) ≈ ( ) ( )( )   and the inner derivative  ( ) ≈−( ) ≈   as well as the convection term  − ( ) ≈   The contributions of the diffusion and convection terms to the coefficients of the algebraic equation (8) are therefore; = + = − − 2( − )( − ) = + = − − − 2( − )( − ) = + = − +  Tridiagonal matrix algorithm is applied for solving linear system of the algebraic equation (8). We set  = 1.0, = 1.0, = 11.  IV. SEQUENCES OF THE PECLET NUMBERS AND THE GRID EXPANSION FACTORS The range of low Peclet numbers is 0,100 . The mathematical relationship between and grid expansion factors is represented by a set of ordered pairs , , =1,2, … , n. We define a sequence of by  , = ⁄ ,= ⁄ ,= ⁄ ,   . . . = ⁄ , where the constants , ∈ ℤ . Next, defining a sequence of by  , = + , = + , = + ,   .  .  . = + , where the constants , ∈ ℤ . Let W and X be the domain and the target of , respectively, where the function from W to X is a collection of ordered pairs of the form ( , ). Note that and are in W and X, respectively. The following conditions need to be satisfied by the collection: Condition a For each in W, there is an element in X such that ( , ) is one of the ordered pairs. In other words, each element in the domain of has a value ( ) under . Condition b If ( , ) and ( , ) are both among the ordered pairs that make up the function, then = . This means that every element of the domain has at most one value under . The function is therefore a mechanism that assigns to each element of the domain a unique element ( ) of the target. We write = ( ), indicating that the ordered pair ( , ) is in the collection of ordered pairs which define the function . Thus the set { ( ): is a real number in W} of values of is the image of . Let = 1, = 6, = 100, = 0.5, = 2 and = 0.1  such that the sequence in (13) and (14) become 100, 50, 25, 12.5, 6.25, 3.125; and 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 respectively. Proposition a The sequences’ elements in (13) and (14), whose boundary values and independent variables are given in (15), form the ordered pairs (Pe, r ) which satisfy Condition a and b such that; ( , ), ( , ), … . . , ( , )= (100, 0.5), (50, 0.6), (25, 0.7),(12.5,0.8),(6.25, 0.9), (3.125, 1.0)  Proposition b If 0 ≤ ≤ 3.125, then the ordered pair ( , ) = ( , 1). Engineering, Technology & Applied Science Research Vol. 8, No. 2, 2018, 2680-2684 2683 www.etasr.com Abdullah: Formulation of Low Peclet Number Based Grid Expansion Factor for the Solution … V. RESULTS OF CALCULATIONS The boundary conditions and other parameter settings for CDE (i.e. equation in (5)) are given in (6), (12), (13), and (15). The concentration profiles which are numerically calculated for of interest are plotted in Figure 3, and show good agreement with the exact solutions. This proves Proposition a. The results vary exponentially in x-direction, and the area under the curve represented by the integral (x) x is inversely proportional to . Note that since r = 1.0 matches Pe = 3.125, then r = 1.0 is appropriate for 0≤ Pe ≤ 3.125 where φ profile is close to linearity with respect to x . The uniform grid is therefore sufficient for the correct prediction of Pe within the range. Proposition b is thus proven. Fig. 3. Concentration profiles In this numerical calculation of a low Peclet number convection-diffusion flow, it is found that the expansion factor is inversely proportional to the logarithm of as shown in Figure 4. Fig. 4. Grid expansion factor as a function of logarithm of the Peclet number Theorem Let 0 ≤ Pe ≤ 100, the grid expansion factor r for solving the convection-diffusion equation in (5), with the flow conditions in (6) and (12), is expressed as a linear function of lg Pe; r = m lg Pe + b, for 3.125 ≤ Pe ≤ 100, and as a constant; r = 1 for 0 ≤ Pe ≤ 3.125 , where m and b are curve slope and a constant, respectively. VI. FINAL REMARKS A new technique in the determination of grid expansion factor which represents a quantitative guideline for the numerical solution of the convection-diffusion equations is proposed. The understanding on the influence of the Peclet number on the grid expansion factor formes a basis for a more effective approach in the selection of grid type for the computational procedure. The key aspect in this research is the formulation of the special function as a collection of ordered pairs of the form , , = 1,2, …, n for the given CDE. This sheds light on the possibility of a more general framework for the selection of grid type in computational fluid dynamics, the relationship between the flow parameter/s and the grid quality in finite difference numerical scheme, as well as the influence of (e.g. low, transition, high) on the numerical error pattern. ACKNOWLEDGMENT Author would like to thank Universiti Tun Hussein Onn Malaysia (UTHM) and Ministry of Higher Education of Malaysia (MoHE) for the research facilities. REFERENCES [1] L. Li, R. Mei, J. F. Klausner, “Lattice Boltzmann models for the convection-diffusion equation: D2Q5 vs D2Q9”, International Journal of Heat and Mass Transfer, Vol. 108, Part A, pp. 41–62, 2017 [2] L. Li, R. Mei, J. F. Klausner, “Multiple-relaxation-time lattice Boltzmann model for the axisymmetric convection diffusion equation”, International Journal of Heat and Mass Transfer, Vol. 67, pp. 338–351, 2013 Engineering, Technology & Applied Science Research Vol. 8, No. 2, 2018, 2680-2684 2684 www.etasr.com Abdullah: Formulation of Low Peclet Number Based Grid Expansion Factor for the Solution … [3] Y. Hu, D. Li, S. Shu, X. Niu, “Lattice Boltzmann flux scheme for the convection-diffusion equation and its applications”, Computers & Mathematics with Applications, Vol. 72, No. 1, pp. 48–63, 2016 [4] H. Wang, B. Shi, H. Liang, Z. Chai, “Finite-difference lattice Boltzmann model for nonlinear convection-diffusion equations”, Applied Mathematics and Computation, Vol. 309, pp. 334–349, 2017 [5] Z. X. Hu, J. Huang, W. X. Huang, G. X. Cui, “Second-order curved interface treatments of the lattice Boltzmann method for convection- diffusion equations with conjugate interfacial conditions”, Computers & Fluids, Vol. 144, pp. 60–73, 2017 [6] M. Bittl, D. Kuzmin, R. Becker, “The CG1-DG2 method for convection- diffusion equations in 2D”, Journal of Computational and Applied Mathematics, Vol. 270, pp. 21–31, 2014 [7] H. Liu, M. Pollack, “Alternating evolution discontinuous Galerkin methods for convection-diffusion equations”, Journal of Computational Physics, Vol. 307, pp. 574–592, 2016 [8] J. Zhang, L. Ge, J. Kouatchou, “A two colorable fourth-order compact difference scheme and parallel iterative solution of the 3D convection diffusion equation”, Mathematics and Computers in Simulation, Vol. 54, No. 1–3, pp. 65–80, 2000 [9] H. Sun, N. Kang, J. Zhang, E. S. Carlson, “A fourth-order compact difference scheme on face centered cubic grids with multigrid method for solving 2D convection diffusion equation”, Mathematics and Computers in Simulation, Vol. 63, No. 6, pp. 651–661, 2003 [10] J. Singh, R. Swroop, D. Kumar, “A computational approach for fractional convection-diffusion equation via integral transforms”, Ain Shams Engineering Journal, in press [11] M. Behroozifar, A. Sazmand, “An approximate solution based on Jacobi polynomials for time-fractional convection-diffusion equation”, Applied Mathematics and Computation, Vol. 296, pp. 1–17, 2017 [12] V. Martin, “An optimized Schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions”, Computers & Fluids, Vol. 33, No. 5-6, pp. 829-837, 2005 [13] Z. F. Tian, P. X. Yu, “A high-order exponential scheme for solving 1D unsteady convectiondiffusion equations”, Journal of Computational and Applied Mathematics, Vol. 235, No. 8, pp. 2477–2491, 2011 [14] R. C. Mittal, R. K. Jain, “Redefined cubic B-splines collocation method for solving convection-diffusion equations”, Applied Mathematical Modelling, Vol. 36, No. 11, pp. 5555–5573, 2012 [15] H. H. Cao, L. Bin Liu, Y. Zhang, S. M. Fu, “A fourth-order method of the convection-diffusion equations with Neumann boundary conditions”, Applied Mathematics and Computation, Vol. 217, No. 22, pp. 9133– 9141, 2011 [16] M. Dehghan, “On the numerical solution of the one-dimensional convection-diffusion equation”, Mathematical Problems in Engineering, Vol. 2005, No. 1, pp. 61–74, 2005 [17] H. S. Shukla, M. Tamsir, “An exponential cubic B-spline algorithm for multi-dimensional convection-diffusion equations”, Alexandria Engineering Journal, in press [18] S. Biringen, “A note on the numerical stability of the convection- diffusion equation”, Journal of Computational and Applied Mathematics, Vol. 7, No. 1, pp. 17–20, 1981 [19] S. Zhai, X. Feng, Y. He, “An unconditionally stable compact ADI method for three-dimensional time-fractional convection-diffusion equation”, Journal of Computational Physics, Vol. 269, pp. 138–155, 2014 [20] Z. Zhou, D. Liang, “The mass-preserving and modified-upwind splitting DDM scheme for time-dependent convection-diffusion equations”, Journal of Computational and Applied Mathematics., Vol. 317, pp. 247– 273, 2017 [21] L. Ge, J. Zhang, “High Accuracy Iterative Solution of Convection Diffusion Equation with Boundary Layers on Nonuniform Grids”, Journal of Computational Physics, Vol. 171, No. 2, pp. 560–578, 2001 [22] Y. Ma, Y. Ge, “A high order finite difference method with Richardsonextrapolation for 3D convection diffusion equation”, Applied Mathematics and Computation, Vol. 215, No. 9, pp. 3408–3417, 2010 [23] S. A. Mohamed, N. A. Mohamed, A. F. Abdel Gawad, M. S. Matbuly, “A modified diffusion coefficient technique for the convection diffusion equation”, Applied Mathematics and Computation, Vol. 219, No. 17, pp. 9317–9330, 2013 [24] S. E. Buitrago Boret, O. J. Jimenez P., “Integrated framework for solving the convection diffusion equation on 2D Quad mesh relying on internal boundaries”, Computers & Mathematics with Applications, Vol. 74, No. 1, pp. 218–228, 2017 [25] N. W. Han, J. Bhakta, R. G. Carbonell, “Longitudinal and lateral dispersion in packed beds: Effect of column length and particle size distribution”, AIChE Journal, Vol. 31, No. 2, pp. 277–288, 1985