CHEMICAL ENGINEERING TRANSACTIONS VOL. 63, 2018 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Jeng Shiun Lim, Wai Shin Ho, Jiří J. Klemeš Copyright © 2018, AIDIC Servizi S.r.l. ISBN 978-88-95608-61-7; ISSN 2283-9216 Achieving a Sustainable Environment using Numerical Method for the Solution of Advection Equation in Fluid Dynamics Olusegun Adeyemi Olaijua,*, Yeak Su Hoea, Ezekiel Babatunde Ogunbodeb, Jonathan Kehinde Fabic, Ernest Ituma Egbad a Department of Mathematics. Faculty of Science, Universiti Teknologi Malaysia, Skudai, 81300. Johor Bahru, Malaysia. b Department of Building. Federal University of Technology Minna, Niger State, Nigeria. c Department of Quantity Surveying. Federal Polytechnic Ilaro. Ogun State. Nigeria. d Department of Technology and Vocational Education, Ebonyi State University, Abakaliki. Ebonyi State. Nigeria. olaiju@yahoo.com The paper compares some explicit numerical methods for solving the advection equation that could be used to describe transportation of a pollutant in one direction. The one-dimensional advection equation is solved by using five different standard finite difference schemes (the Upwind, FTCS, Lax- Friedrichs, Lax wendroff and Leith’s methods) via C codes. An example is used for comparison; the numerical results are compared with analytical solution. It is found that, for the example studied, all the explicit finite difference methods give satisfactory solutions to the advection problem. 1. Introduction Advection can be referred to as a lateral or horizontal transfer of mass, heat, or other properties or movement of some material dissolved or suspended in fluid. This definition suggested winds that blow across earth's surface to represent air advection. In the case of ocean, the currents are example of advection. There are a lot of debate currently by geologists on the presence and role of substantial advective processes on earth. Meteorologist rely on accurate numerical approximations of the advection equation for weather forecasting. The advection equation also governs acoustic or elastic wave propagation, gas discharge problems in physics and the motion of shallow free surface flows. The main applications of advection-diffusion equation are in fluid dynamics, heat transfer, and mass transfer (Appadu, 2013). The linear advection equation is a one-dimensional linear Partial Differential Equation (PDE) that describes one-dimensional transmission of mass, heat or other properties from one part of the fluid to another part by the movement of the fluid itself. Since, the exact solutions of these differential equations are sometimes difficult to get and may not even exist in some cases, the necessity to approximate their solution numerically is expedient. Oosthuizen (2009) used an approximate numerical model of convective heat transfer on a window mounted to natural convective heater to study how the relative temperature of the heater affects the convective heat transfer rate to the window and the flow in the room. Numerical results were obtained by using a commercial CFD code. Vonda and Hajek (2009) conducted an experimental and numerical analysis on wall heat transfer in non-premixed gas combustor. They presented CFD simulations of the combustor using Reynolds-average navier-stokes models and validated them by the measured data and summarized the overall model configuration, boundary conditions and key sub models for turbulence and chemistry. They discussed the derivations from the experiment and offered possible explanations. This paper compares five explicit numerical methods for solving the advection equation, which can be used to describe fluid transportation in one direction. The one-dimensional advection equation is solved by using five different standard finite difference schemes via C++ codes. DOI: 10.3303/CET1863106 Please cite this article as: Olusegun Adeyemi Olaiju, Yeak Su Hoe, Ezekiel Babatunde Ogunbode, Jonathan Kehinde Fabi, Ernest Ituma Egba, 2018, Achieving a sustainable environment using numerical method for the solution of advection equation in fluid dynamics, Chemical Engineering Transactions, 63, 631-636 DOI:10.3303/CET1863106 631 1.1 Advection equation The one dimensional linear advection is given by Eq(1). u u a 0, 0 x 1, 0 t T t x ∂ ∂ + = < < < ≤ ∂ ∂ (1) With initial condition; ( )u x,0 f (x), 0 x 1= ≤ ≤ , and boundary conditions; ( ) 0u 0,t g (t),= 0 t T ,< ≤ ( ) 1u 1,t g (t), 0 t T ,= < ≤ The value a in Eq(1) is the velocity components of the fluid in the directions of x. The term ua x ∂ ∂ is called the advective velocity and u is the unknown function of x and t. Numerical methods are usually used to solve partial differential equation. These methods are not just simple tools for the purpose of application only, there are many important theoretical concepts associated with numerical methods (Thomas, 2013). Numerical methods such as the finite difference method, the finite element method, the finite volume method and the boundary value method are useful in solving physical problems. Out of all these methods, the finite difference method are straightforward as well as economical and easier to implement compared to the other numerical methods. This paper focuses on the finite difference method. The main idea behind finite difference method for approximating the solution of Eq(1), is to replace the continuous PDE by a discrete algebraic equation (William, 1992). There are two categories of Finite difference scheme, the implicit and explicit scheme. Implicit schemes enjoy the advantage of better stability properties over explicit schemes. Moreover, explicit schemes have better advantages over implicit via efficiency and ease of parallelization. The paper presents the features and characteristics of some finite difference schemes for the solution of linear advection equation. 1.2 Finite difference approximations Let u u( x, y ),= by dividing the first quadrant of the xy plane into uniform rectangular grid lines parallel to both the x and y-axis as in Figure 1, where the grid size is ∆x and ∆y on both x and y axis respectively. By using Eq(2) which is also referred to as the Taylor polynomials: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 3 i j i j x i j xx i j xxx i j 2 3 i j i j x i j xx i j xxx i j x x x u x x, y u x , y u x , y u x , y u x , y ... 1 ! 2 ! 3 ! x x x u x x, y u x , y u x , y u x , y u x , y ... 1 ! 2 ! 3 !  Δ Δ Δ + Δ = + + + +     Δ Δ Δ + Δ = + + + +   (2) and by considering Figure 1, the squares at the horizontal axis indicate the location of the initial values that are already known, while the squares along the vertical axis indicate the location of the boundary values, which are also known. The black dots indicate the position of the interior points where the finite difference approximation is to be computed. Figure 1: Discrete Grid Points 632 The approximations to x xx y yyu ,u u u about the point ( )i ju x , y are given by (Luciano, 2010) to be: Forward difference approximation for yu shown in Eq(3), Backward difference approximation for yu in Eq(4), ( ) ( ) ( )i j y i jy i j y u x , y u x , y u x , y + Δ − = Δ (3) ( ) ( ) ( )i j i j yy i j y u x , y u x , y u x , y − − Δ = Δ (4) Central difference approximation for yu shown in Eq(5), Central difference approximation yyu in Eq(6), ( ) ( ) ( )i j y i j yy i j y u x , y u x , y u x , y 2 + Δ − − Δ = Δ (5) ( ) ( ) ( ) ( ) 2 i j y i j i j y i j y u x , y 2u x , y u x , y x , y + Δ − + − Δ = Δ (6) Forward difference approximation for u is shown in Eq(7), while the Backward difference approximation for u is presented in Eq(8), ( ) ( ) ( )i j i jx i j u x x, y u x , y u x , y x + Δ − = Δ (7) ( ) ( ) ( )i j i jx i j u x , y u x x, y u x , y x − − Δ = Δ (8) Central difference approximation u is shown in Eq(9), ( ) ( ) ( )i j i jx i j u x x, y u x x, y u x , y 2 x + Δ − − Δ = Δ (9) Central difference approximation u is presented in Eq(10) ( ) ( ) ( ) ( )i j i i i jxx i j 2 u x x, y 2u x y u x x , y u x , y x + Δ − + − Δ = Δ (10) Eq(3), Eq(4), Eq(7) and Eq(8) are said to be first order accurate while Eq(5), Eq(6), Eq(9) and Eq(10) are said to be second order accurate. An initial condition and a boundary condition is required for the solution of Eq(1). The initial condition is one that prescribes the unknown function throughout a given region at some initial time t, usually taken as t = 0. The boundary condition on the other hand describes the function for all time t at the prescribed boundary given. There are two common boundary conditions. The Dirichlet boundary condition and the Neumann boundary condition. Dirichlet boundary condition takes account of the value of the function specified on the boundary. The Neumann boundary condition describes the value of the derivative normal to the boundary specified on the boundary (Thomas, 2013). Boundary conditions are required to be specified at each point on the boundary of the closed solution domain. Different portions of the boundary may have different types of boundary conditions specified (Morton and Mayers, 1994). 2. Consistency If in the limit that a grid spacing tends to zero, a system of finite difference equations generated by discretization process is equivalent to the original PDE at each point, then that system of finite difference equation is said to be consistent with the original PDE. Consistency is necessary for the approximate solution to converge to the solution of the PDE under consideration. Even though it can be quite tedious and messy, consistency is straightforward to demonstrate (Strikwerda, 2004). 3. Stability and Convergence A numerical scheme is said to be stable if errors from any source do not pile up as the calculation proceeds. A scheme is also considered stable if its solution remains a uniformly bounded function of the initial state for all sufficiently small ∆t. For problems that are dependent on time, stability ensures that the method produces a bounded solution if the exact solution is bounded in itself. Stability is very important in numerical analysis. (Burden and Faires, 1997). The finite difference scheme used to approximate a differential equation is said to 633 be convergent if the computed solution of the discretized equations tends to the exact solution of the differential equation as the grid and time spacing tend to zero (Louise, 2015). 4. Numerical dissipation and dispassion By considering Eq(1), with the associated initial and boundary conditions, assume the initial condition is a periodic function. The exact solution is a travelling wave at constant speed u. Numerical dissipation is a situation where the peaks of the wave are not maintained but diminish in time. Numerical dispersion means that packet of waves travel at speed different from that of u (Louise, 2015). 5. Courant, Friedrichs and Lewy stability criterion (CFL) Courant et al. (1967) in their study on the partial difference equations of mathematical physics laid the theoretical foundations for practical finite difference computations (Trefethen, 1996). The paper identified a fundamental necessary condition for convergence of any numerical scheme, which is known as the Courant, Friedrichs and Lewy (CFL) stability condition or criterion. The Courant number is t a x Δ Δ , where a is the characteristic wave speed of the system, ∆t is the time-step of the numerical model, and ∆x is the spacing of the grid in the numerical model. The Courant condition is: t a 1 x Δ ≤ Δ (11) This condition means that a fluid particle should not travel more than one spatial step size ∆ in one time step tΔ . Courant, Friedrichs and Lewy proved that the solution of the finite different system converges to that of the partial difference equation as xΔ and tΔ tend to zero provided the domain of dependence of the partial different equation lies inside that of the PDE (Arnold, 2015). Some of the explicit schemes for linear advection equation are the Upwind, FTCS, Lax- Friedrichs, Lax wendroff and Leith’s methods. Luciano (2011) gives details of their derivation, stability properties and truncation errors. Generally, explicit schemes are constrained by the CFL. For the one-dimensional linear advection equation, solutions ranging from the exact to highly unsatisfactory can be obtained using different discretization schemes. Depending on the choice of the parameter, even for the same discretization scheme, the nature of the solution may be satisfactory or unsatisfactory. In this case, the parameters are the values of ∆t and ∆x. Such a sensitivity of the solution to the discretization scheme used and the choice of tΔ and xΔ is deemed unacceptable when dealing with a general case whose exact solution may not be known (Recktenwald, 2011). Convectional explicit finite difference schemes for the advection equations are subject to the time step restrictions dictated by the CFL conditions. Time step sizes are mostly chosen to satisfy the CFL condition rather than to satisfy accuracy requirements (Luciano, 2010). 8. C++ code In scientific computing, the use of C++ has grown significantly recently. These is due to the fact that C++ supports both object-base and object-oriented programming and also template programming and attractive syntax via operator overloading (Acklam et al., 1998). These programming features make software development faster and more reliable. Although, in the scientific computing community, the computational efficiency of C++ has always been of some concern. According to Acklam et al. (1998), C++ can run very fast with careful coding. 9. Test example In this section, we apply some of the explicit numerical schemes to a test example in order to compare their relative accuracy via C++ codes. The hardware used was HP Notebook, Intel(R) Core(TM) i5-6200U CPU@ 2.40GHz, RAM 4.00GB, 64-bit operating system. x64-base processor and the C++. Software used was Code: Blocks version 16.01. Considering a one dimensional linear advection Eq(1) with initial condition ( ) ( )( )= − − ≤ ≤2u x,0 exp 100 x 0.5 , 0 x 1 (12) and boundary conditions ( )u 0,t 0,= 0 t T ,< ≤ ( )u 1,t 0 t T .= < ≤ Setting Lt 0.03,= Nt 1000,= j i Lt Lxt , Lx 1.0, Nx 100, x , t j t, 0 j Nt , x i x, 0 i Nx and a 1. Nt Nx Δ = = = Δ = = Δ ≤ ≤ = Δ ≤ ≤ = The exact solution of convection Eq(1) is; 634 ( ) ( ) ( )( )u x,t exp 100 x 0.5 2at x 0.5 .= − − + − (13) The 1D Upwind scheme Upwind schemes represent a class of numerical discretization approaches for resolving hyperbolic partial differential equations. Upwind schemes use an adaptive or solution-sensitive finite difference stencil to numerically simulate the course of dissemination of information in a flow field. ( ) ( )j 1 j j j j 1 j j ji i i i 1 i i i 1 it tu u a u u , u u a u u x x + + − + Δ Δ = − − = − − Δ Δ (14) This technique is stable for t 0 a 1 x Δ < ≤ Δ . (15) The 1D FTCS scheme ( )j 1 j j ji i i 1 i 1tu u a u u . 2 x + + − Δ = − − Δ (16) This scheme is unconditional instable (Luciano, 2010). The 1D Lax-Friedrichs scheme ( ) ( )+ + − + −Δ= + − −Δ j 1 j j j j i i 1 i 1 i 1 i 1 1 t U u u a u u 2 2 x (17) This technique is conditionally stable for Δ < ≤ Δ t 0 a 1 x , but suffers from a considerable dissipation (Luciano, 2010). The Leith's Scheme + + +      Δ Δ Δ Δ Δ     = + + − + − +               Δ Δ Δ Δ Δ           2 2 2 j 1 j j j i i 1 i i 1 1 t t t 1 t t U a a u 1 a u a a u 2 x x x 2 x x (18) This technique is stable for Δ < ≤ Δ t 0 a 1 x (Luciano, 2010) as given in Eq(11). The 1D Lax-Wendroff scheme ( ) ( )+ + − + − Δ Δ  Δ Δ= − − + − + 2 j 1 j j j j j j i i i 1 i 1 i 1 i i 1 i tt aa xxu u u u u 2u u u 2 2 (19) This technique is stable for Δ  ≤ Δ  2 t a 1 x (Luciano, 2010) as given in Eq(11). Considering the performance of the numerical schemes, it was observed that the five numerical schemes magnificently performed well in terms of computational time. Figure 2 shows the numerical solutions to the problem of one-dimensional advection equation. The figure indicates that almost all of the proposed numerical schemes solved the linear convection equation quite satisfactorily in exemption of Friedrichs scheme. Figure 3 shows the error in the numerical solutions from each of the methods when compared with the analytical solution, for N = 1000, dt = 0.00003 and dx = 0.01. Comparisons are made for the solutions at j = 50. The errors are all small except for the Friedrichs scheme. The error for the Friedrichs scheme is larger in magnitude than the errors for the other four methods. From Figure 2, it is shown that Lax-Friedrichs scheme is stable, but it suffers from a considerable dissipation. The truncation error is also expressed for the differential operator, the numerical algorithms will not be expressed in terms of the differential operators and will therefore have different (usually smaller) truncation errors. Accuracy is not the most important requirement in numerical analysis and a first-order but stable scheme is greatly preferable to one which is higher order (i.e., has a smaller truncation error) but is unstable. 635 Figure 2: Advection of a 1-d Gaussian pulse for all the numerical schemes at j=50 Figure 3: Errors of numerical calculation for all the numerical schemes at j=50 10. Conclusion Five explicit finite difference methods for one dimensional advection equation have been presented. For the test examples studied, it has been found that the methods gave good point-wise solutions to the test problem. The Upwind method yields satisfactory results only in the case in which the equations have an obvious transport character in one direction like the problem solved in this paper. In more general situations, the Upwind method may not be adequate. This work can be extended to the case of consideration of 1D nonlinear, 2D linear and 2D nonlinear convection diffusion problems and appropriate optimisation techniques can also be used to choose parameters ∆ and ∆ for accurate numerical solution and minimal numerical dissipation. References Acklam E., Jacobsen A., Langtangen H.P., Bruaset A.M., 1998, Optimizing C++ Code for explicit finite difference schemes, Oslo Scientific Computing Archive Report, Oslo, Norway. Appadu A.R., 2013, Numerical solution of the 1D advection-diffusion equation using standard and nonstandard finite difference schemes, Journal of Applied Mathematics, 13 (10), 1-14. Arnold Douglas N., 2015, Lecture notes on numerical analysis of partial differential equations accessed 15.09.2017. Burden R.L., Faires J.D., 1997, Numerical Analysis, Brooks-Cole Publishing Company, New York, United States. Courant R., Friedrichs K., Lewy H., 1967, On the partial difference equations of mathematical physics, IBM Journal of Research and Development, 11 (2), 215–234. Recktenwald G.W., 2004, Finite-difference approximations to the heat equation, Mechanical Engineering Department, Portland State University, Portland Oregon, USA. Trefethen L.N., 1996, Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations accessed 19.09.2017. Louise O., 2015, Numerical solution of partial differential equations accessed 14.09.2017. Luciano R., 2010, Numerical methods for the solution of partial differential equations, Lecture notes for the COMPSTAR, School on Computational Astrophysics, Caen, France. Morton K.W., Mayers D.F., 1994, Numerical Solution of Partial Differential Equations: An Introduction, Cambridge University Press, Cambridge, England. Oosthuizen P.H., 2009, Convective heater on the heat transfer rate from a cold recessed window, Journal of Chemical Engineering Transactions, 18, 69-74. Strikwerda J.C., 2004, Finite difference schemes and partial differential equations, SIAM, Philadelphia, United States. Thomas J.W., 2013, Numerical partial differential equations: finite difference methods, Springer Science and Business Media, 22, 1-437. Vonda J., Hajek J., 2009, Experimental and numerical analysis of wall heat transfer in non-premixed gas combustor, Journal of Chemical Engineering Transactions, 18, 587-592. William F., 1992, Numerical Methods for Partial Differential Equations, Academic Press, Inc., Boston, United States. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x U Advection of a 1-d Gaussian pulse. Analytic FTCS Friedrichs Leith Wendroff Upwind 0 0.2 0.4 0.6 0.8 1 -0.3 -0.2 -0.1 0 0.1 0.2 x E rr o r Errors of the Numerical methods FTCS Friedrichs Leith wendroff Upwind 636