• Views 863
• Citations 2
• ePub 32
• PDF 462
`Journal of Applied MathematicsVolume 2014, Article ID 396738, 8 pageshttp://dx.doi.org/10.1155/2014/396738`
Research Article

## Accurate Simulation of Contaminant Transport Using High-Order Compact Finite Difference Schemes

Department of Civil Engineering, Faculty of Engineering, Pamukkale University, 20070 Pamukkale, Denizli, Turkey

Received 24 January 2014; Revised 27 March 2014; Accepted 10 April 2014; Published 29 April 2014

Copyright © 2014 Gurhan Gurarslan. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Numerical simulation of advective-dispersive contaminant transport is carried out by using high-order compact finite difference schemes combined with second-order MacCormack and fourth-order Runge-Kutta schemes. Both of the two schemes have accuracy of sixth-order in space. A sixth-order MacCormack scheme is proposed for the first time within this study. For the aim of demonstrating efficiency and high-order accuracy of the current methods, some numerical experiments have been done. The schemes are implemented to solve two test problems with known exact solutions. It has been exhibited that the methods are capable of succeeding high accuracy and efficiency with minimal computational effort, by comparisons of the computed results with exact solutions.

#### 1. Introduction

Transport of sediments and contaminants has long been one of the great concerns to hydraulic and environmental engineers. Sediment particles in alluvial rivers are subject to random and complex movement. Understanding the transport of sediment particles is of fundamental and practical importance to hydraulic engineering. Accurate simulation of suspended sediment transport is essential for water quality management, environmental impact assessment, and design of hydraulic structures. Among others, the advection-dispersion equation is crucial to the simulation of suspended sediment transport, contaminant transport in groundwater, and water quality in rivers. Therefore, improving the efficiency and accuracy of numerical schemes for the advection-dispersion equation has been a focus of research [1]. The analytical solutions of advection-dispersion equation can be obtained for limited number of initial and boundary conditions making some simplifying assumptions. But, the usage of analytical solutions in field applications is rather limited because ideal conditions could not generally be satisfied.

Remarkable research studies have been conducted in order to solve advection-dispersion equation numerically like method of characteristic with Galerkin method [2], finite difference method [35], high-order finite element techniques [6], high-order finite difference methods [720], green element method [21], cubic B-spline [22], cubic B-spline differential quadrature method [23], method of characteristics integrated with splines [2426], Galerkin method with cubic B-splines [27], Taylor collocation and Taylor-Galerkin methods [28], B-spline finite element method [29], least squares finite element method (FEMLSF and FEMQSF) [30], lattice Boltzman method [31], Taylor-Galerkin B-spline finite element method [32], and meshless method [33, 34].

Widely used discretization scheme for the numerical solution of hyperbolic partial differential equations is the MacCormack (MC) scheme [35] which is an explicit and two step predictor-corrector schemes. MC scheme is equivalent to the Lax-Wendroff scheme regarding linear equations. MC scheme does not give diffusive errors in the solution while first-order upwind scheme does. This procedure provides the reasonably accurate results and needs less CPU time. Several advantages of the MC scheme make the method a popular choice in computational hydraulics problems. Firstly, the scheme is a shock-capturing technique with second-order accuracy both in time and space. Secondly, the inclusion of the source terms is relatively simple. Thirdly, implementing it in an explicit time-marching algorithm is convenient [36]. This scheme has been successfully applied to unsteady open channel flows [3740], overland flows [41], and contaminant transport [12, 4244]. To be able to solve many problems accurately, using high-order numerical methods is necessary. The idea of using MC schemes with compact finite difference schemes was suggested for the first time by Hixon and Turkel [45]. In the corresponding study, two different fourth-order compact MC schemes were suggested. However in this study, a sixth-order compact MacCormack scheme (MC-CD6) which is structurally different than Hixon and Turkel schemes was proposed. MC-CD6 is applied to the contaminant transport problem in this study for the first time. Another scheme used in this study is RK4-CD6 scheme which is formed by combining a fourth-order Runge-Kutta (RK4) scheme and a sixth-order compact finite difference scheme (CD6) in space. This scheme was applied to the solution for one-dimensional contamination transport problem by Gurarslan et al. [19]. Gurarslan et al. [19] has declared that the RK4-CD6 scheme is very accurate solution approach in solving one-dimensional contaminant transport equation for low and moderate Peclet numbers, that is, . Using the related scheme for two-dimensional contaminant transport problem took place within this study for the first time. Examples of both one- and two-dimensional advection-dispersion problems will be used to investigate accuracy of the RK4-CD6 and MC-CD6 scheme. Numerical results obtained from these examples will be compared to available analytical and/or numerical results existing in the literature.

#### 2. Governing Equation

Two-dimensional advection-dispersion equation in the conservative form is given as follows: where is concentration of a tracer without deposition or degradation; and are space coordinates; and are depth-averaged horizontal fluid velocity components in - and -direction, respectively; and are dispersion coefficients in - and -direction, respectively; and is time. In case of applying classical finite difference schemes as a solution, fundamental difficulty is encountered with arises from the advection term which causes oscillations in solution [1]. Sixth-order compact finite-difference schemes are enhanced to overcome this existing problem. For approximating time derivative, MacCormack and Runge-Kutta schemes are used.

#### 3. Compact Finite Difference Schemes

In this section, compact schemes whose solution and various order derivatives are assumed as unknowns are introduced. An implicit equation including implicit derivatives and functions helps us to solve the derivatives at grid points. Two essential properties of compact schemes can be expressed as high spectral accuracy and relatively compact stencils which correlate the derivatives with function. Compact high-order schemes are closer to spectral methods and they maintain the freedom to retain accuracy in complex geometries, as well. Details about derivation of compact finite difference schemes can be obtained from [46, 47].

##### 3.1. Spatial Discretization

Compact finite difference schemes are used to evaluate spatial derivatives. For any scalar pointwise value , the derivatives of are reached by solving a linear equation system. For derivation of such a formula, great amount of work has been done [46]. When two-dimensional problem is considered, one needs to approximate both partial derivatives in and . The approximation is automatically carried out by using an equal number of grid points in both directions. If value is fixed, approximation of all partial derivatives with respect to is done by using the compact scheme. If value is fixed, approximation of all partial derivatives with respect to is done.

The formulation of first derivative with respect to at internal nodes can be expressed as follows [48]: where , , , , , and is grid size in x-direction. If , the scheme is fifth-order compact upwind scheme; if , it is reduced to sixth-order central compact scheme. The suggested value for is , and the corresponding fifth-order compact upwind scheme is [48]

The formulation of second derivative with respect to at internal nodes can be expressed as follows [46]:

Regarding the nodes close to boundary, approximation formulae of derivatives of nonperiodic problems can be derived by evaluating one-sided schemes. One can find further details about derivations for the first- and second-order derivatives in [46]. The derived formulae at boundary points are given as follows.

The third-order formulae at boundary point 1

The fourth-order formulae at boundary points 2 and

The third-order formulae at boundary point N,

Using formulae given above will result in following matrix equation: where , for all fixed . Here, resembles the number of grid points in each direction. Similarly, the formulae for -direction at boundary and internal points can be derived readily with all fixed .

##### 3.2. Temporal Discretization

In order to solve advection-dispersion equation, MC and RK4 schemes are used. Utility of the compact finite difference method to (1) gives rise to the following differential equation in time: where indicates a spatial differential operator. Compact finite difference formulae are used to approximate the spatial derivatives. Using the compact finite difference formulae enables us to obtain each spatial derivative on the right hand side of (9) and semidiscrete Equation (9) has been solved by the help of MC and RK4 schemes. Solution domain is discretized as to be equally spaced grids for numerical solutions of the problem with the taken boundary and initial conditions using the current scheme.

###### 3.2.1. MacCormack Scheme

MC scheme is a second-order accurate explicit scheme in both time and space, and composed of predictor and corrector steps. For approximating first-order spatial derivatives, first-order backward finite difference formula is used in the predictor step and first-order forward finite difference formula is being used in the corrector step. For approximating second-order spatial derivatives, second-order central finite difference formula is being used in both steps. The semidiscrete Equation (9) is solved by using MC scheme through the operations as follows:

In this study, for approximating first-order spatial derivatives, fifth-order backward compact finite difference formula is used in predictor step and fifth-order forward compact finite difference formula is used in the corrector step. For resolving second-order spatial derivatives, sixth-order central compact difference equations are used in both steps. An accurate finite difference scheme (MC-CD6) which is sixth-order in space and second-order in time is obtained.

###### 3.2.2. Runge-Kutta Scheme

Another time-integration scheme which was used in this study is RK4 scheme. In this scheme, a sixth-order central compact finite difference formula is used for approximating first-, and second-order spatial derivatives. Steps of RK4 scheme are given below:

#### 4. Numerical Applications

To be able to demonstrate behavior and capability of the present schemes, computational experiments were performed in this section. Checking accuracy of the methods was achieved by applying current methods for different grid size and time step values. Some codes produced in MATLAB 7.0 enabled us to carry out all computations.

Example 1. For solving the advection-dispersion equation, a straight prismatic channel in which the water flows at constant velocity was used. Channel length was taken as m and the channel is divided into intervals of constant length m. It is assumed in this example that flow velocity and dispersion coefficients are to be m/s and m2/s. These circumstances lead to the propagation of a steep front, that is, simultaneously subjected to the dispersion. Analytical solution of the advection-dispersion equation is given below [49]:

At the boundaries, the following conditions are used:

Initial conditions can be taken from exact solution. Table 1 exhibits the comparison between numerical solutions and exact solution. Table 1 apparently shows that solutions obtained for s with FEMLSF [30] and FEMQSF [30] do not sufficiently converge, for . It is proven by this status that selected time step of these methods is larger than what it needs to be. Because of the fact that solution for is not accurate enough, calculations have been done for situations of and , and corresponding results were compared with FEMLSF, FEMQSF, and RK4-CD6 [19]. As the results of the schemes for are considered on acceptable level, the results obtained by MC-CD6 and RK4-CD6 schemes for are same with exact solution. Errors of these two schemes ( norm error and norm error) are quite close to each other. As seen again in Table 1, the CPU time required for MC-CD6 scheme is less with respect to RK4-CD6 scheme. Thus, both RK4-CD6 and MC-CD6 schemes can be safely used in solving one-dimensional contaminant transport problems.

Table 1: Comparison between numerical solutions and the exact solution.

Figures 1 and 2 show comparison of exact solution and the numerical solution obtained by using MC-CD6 scheme for s and s (, ). Figures 1 and 2 prove that there arises an excellent agreement between MC-CD6 and exact solutions.

Figure 1: Comparison of the exact solution and the numerical solution obtained with MC-CD6 scheme for s.
Figure 2: Comparison of the exact solution and the numerical solution obtained with MC-CD6 scheme for s.

Example 2. Let (1) for, and domain , evaluated with initial condition presented below,

The exact solution is given by [16]and the appropriate boundary conditions can easily be obtained from the exact solution. Consider

Initial condition which is a Gaussian pulse and having a pulse height of 1 is centered at (0.5, 0.5). Figures 3 and 4 exhibit initial pulse and the pulse at obtained through the RK4-CD6 scheme. After a time period of 1.25 sec, Gaussian pulse moves to a position centered at (1.5, 1.5) with a pulse height of 1/6. Parameters are taken as , , and in Table 2. value is taken as 0.00625 in order to obtain average absolute and errors. Table 2 apparently exhibits that the errors obtained by using the RK4-CD6 are far smaller when it is compared to the literature. CPU time values required for RK4-CD6 and MC-CD6 schemes are found as 13.98 sec and 6.90 sec, respectively. Although MC-CD6 scheme requires less CPU time than RK4-CD6 scheme, it is apparently seen in Table 2 that MC-CD6 does not produce more accurate results than RK4-CD6. When is chosen, the value of average is obtained as and error is obtained as with MC-CD6. But in this case, CPU time required for MC-CD6 scheme is calculated as 65.40 sec. Therefore, using RK4-CD6 scheme is suggested for solution of two-dimensional contaminant transport problems.

Table 2: Comparison of average absolute and maximum absolute errors with the literature.
Figure 3: Initial Gaussian pulse of Example 2.
Figure 4: The RK4-CD6 solution of Example 2 with , , , and for , .

Table 3 presents the pulse height values obtained for the parameters , , and by using various time steps. Kalita et al. [16] have used three different compact schemes in their studies. Obtained results are compared with results of (9,5), (5,9), and (9,5) schemes of Kalita et al. [16]. Table 3 proves that pulse height values of the RK4-CD6 scheme is more accurate than the results of the (5,9), (9,5), and (9,9) schemes, despite the fact that the results of MC-CD6 scheme are accurate at acceptable level. Figure 5 shows contour lines of the RK4-CD6 solutions in the domain , with the parameters , , , and .

Table 3: Pulse height values of Example 2 for various values of with , , and .
Figure 5: Contour lines of the RK4-CD6 solution and absolute errors in the domain , with , , , and .

#### 5. Conclusions

Throughout this study, high-order compact finite difference schemes composed of second-order MacCormack and fourth-order Runge-Kutta time integration schemes have been used to be able to perform numerical simulation of one- and two-dimensional advective-dispersive contaminant transport. For demonstrating efficiency and high-order accuracy of the current methods, numerical experiments have been done. Then, the schemes are implemented for solving two test problems which have known exact solutions. It has been shown that the used methods are capable of succeeding high accuracy and efficiency with minimal computational effort, supported by comparisons of the computed results with exact solutions.

In solution for one-dimensional contaminant transport problem, it was seen that the error values obtained with RK4-CD6 and MC-CD6 schemes and the required CPU time values are close to each other. Whereas in solution for two-dimensional contaminant transport problem, it was observed that RK4-CD6 scheme is stable for great values and produces better results than MC-CD6 scheme. When value is decreased, it was determined that MC-CD6 scheme gives fine results but required CPU time value considerably increases. RK4-CD6 scheme has produced better results than the studies given in literature in solution for both one- and two-dimensional contaminant transport problem. The proposed schemes produce convergent approximations for the contaminant transport problems having low and moderate Pe number. Obtaining the solutions for contaminant transport problem in higher Peclet numbers by using compact upwind schemes was left to further studies.

#### Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

#### References

1. C. Man and C. W. Tsai, “A higher-order predictor-corrector scheme for two-dimensional advection-diffusion equation,” International Journal for Numerical Methods in Fluids, vol. 56, no. 4, pp. 401–418, 2008.
2. R. Szymkiewicz, “Solution of the advection-diffusion equation using the spline function and finite elements,” Communications in Numerical Methods in Engineering, vol. 9, no. 3, pp. 197–206, 1993.
3. H. Karahan, “Implicit finite difference techniques for the advection-diffusion equation using spreadsheets,” Advances in Engineering Software, vol. 37, no. 9, pp. 601–608, 2006.
4. H. Karahan, “Unconditional stable explicit finite difference technique for the advection-diffusion equation using spreadsheets,” Advances in Engineering Software, vol. 38, no. 2, pp. 80–86, 2007.
5. H. Karahan, “Solution of weighted finite difference techniques with the advection-diffusion equation using spreadsheets,” Computer Applications in Engineering Education, vol. 16, no. 2, pp. 147–156, 2008.
6. H. S. Price, J. C. Cavendish, and R. S. Varga, “Numerical methods of higher-order accuracy for diffusion-convection equations,” Society of Petroleum Engineers, vol. 8, no. 3, pp. 293–303, 1968.
7. K. W. Morton, “Stability of finite difference approximations to a diffusion-convection equation,” International Journal for Numerical Methods in Engineering, vol. 15, no. 5, pp. 677–683, 1980.
8. H. Karahan, “A third-order upwind scheme for the advection-diffusion equation using spreadsheets,” Advances in Engineering Software, vol. 38, no. 10, pp. 688–697, 2007.
9. Y. Chen and R. A. Falconer, “Advection-diffusion modelling using the modified QUICK scheme,” International Journal for Numerical Methods in Fluids, vol. 15, no. 10, pp. 1171–1196, 1992.
10. M. Dehghan, “Weighted finite difference techniques for the one-dimensional advection-diffusion equation,” Applied Mathematics and Computation, vol. 147, no. 2, pp. 307–319, 2004.
11. M. Dehghan and A. Mohebbi, “High-order compact boundary value method for the solution of unsteady convection-diffusion problems,” Mathematics and Computers in Simulation, vol. 79, no. 3, pp. 683–699, 2008.
12. G. Gurarslan and H. Karahan, “Numerical solution of advection-diffusion equation using a high-order MacCormack scheme,” in Proceedings of the 6th National Hydrology Congress, Denizli, Turkey, September 2011.
13. M. Sari, G. Gürarslan, and A. Zeytinoğlu, “High-order finite difference schemes for solving the advection-diffusion equation,” Mathematical and Computational Applications, vol. 15, no. 3, pp. 449–460, 2010.
14. A. Mohebbi and M. Dehghan, “High-order compact solution of the one-dimensional heat and advection-diffusion equations,” Applied Mathematical Modelling, vol. 34, no. 10, pp. 3071–3084, 2010.
15. A. A. Mohamad, “Spatially fourth-order-accurate scheme for unsteady-convection problems,” Numerical Heat Transfer B: Fundamentals, vol. 31, no. 3, pp. 373–385, 1997.
16. J. C. Kalita, D. C. Dalal, and A. K. Dass, “A class of higher order compact schemes for the unsteady two-dimensional convection-diffusion equation with variable convection coefficients,” International Journal for Numerical Methods in Fluids, vol. 38, no. 12, pp. 1111–1131, 2002.
17. S. Karaa and J. Zhang, “High order ADI method for solving unsteady convection-diffusion problems,” Journal of Computational Physics, vol. 198, no. 1, pp. 1–9, 2004.
18. Z. F. Tian and Y. B. Ge, “A fourth-order compact ADI method for solving two-dimensional unsteady convection-diffusion problems,” Journal of Computational and Applied Mathematics, vol. 198, no. 1, pp. 268–286, 2007.
19. G. Gurarslan, H. Karahan, D. Alkaya, M. Sari, and M. Yasar, “Numerical solution of advection-diffusion equation using a sixth-order compact finite difference method,” Mathematical Problems in Engineering, vol. 2013, Article ID 672936, 7 pages, 2013.
20. B. J. Noye and H. H. Tan, “Finite difference methods for solving the two-dimensional advection-diffusion equation,” International Journal for Numerical Methods in Fluids, vol. 9, no. 1, pp. 75–98, 1989.
21. A. E. Taigbenu and O. O. Onyejekwe, “Transient 1D transport equation simulated by a mixed green element formulation,” International Journal for Numerical Methods in Fluids, vol. 25, no. 4, pp. 437–454, 1997.
22. R. C. Mittal and R. K. Jain, “Numerical solution of convection-diffusion equation using cubic B-splines collocation methods with Neumann's boundary conditions,” International Journal of Applied Mathematics and Computation, vol. 4, no. 2, pp. 115–127, 2012.
23. A. Korkmaz and I. Dag, “Cubic B-spline differential quadrature methods for the advection-diffusion equation,” International Journal of Numerical Methods for Heat & Fluid Flow, vol. 22, no. 8, pp. 1021–1036, 2012.
24. F. M. Holly Jr. and A. Preissmann, “Accurate calculation of transport in two dimensions,” Journal of Hydraulic Division, vol. 103, no. 11, pp. 1259–1277, 1977.
25. T.-L. Tsai, J.-C. Yang, and L.-H. Huang, “Characteristics method using cubic-spline interpolation for advection-diffusion equation,” Journal of Hydraulic Engineering, vol. 130, no. 6, pp. 580–585, 2004.
26. T.-L. Tsai, S.-W. Chiang, and J.-C. Yang, “Examination of characteristics method with cubic interpolation for advection-diffusion equation,” Computers & Fluids, vol. 35, no. 10, pp. 1217–1227, 2006.
27. L. R. T. Gardner and I. Dag, “A numerical solution of the advection-diffusion equation using B-spline finite element,” in Proceedings of the International AMSE Conference “Systems Analysis, Control & Design”, vol. 1, pp. 109–116, Lyon, France, July1994.
28. I. Dag, A. Canivar, and A. Sahin, “Taylor-Galerkin method for advection-diffusion equation,” Kybernetes, vol. 40, no. 5-6, pp. 762–777, 2011.
29. S. Dhawan, S. Kapoor, and S. Kumar, “Numerical method for advection diffusion equation using FEM and B-splines,” Journal of Computational Science, vol. 3, no. 5, pp. 429–437, 2012.
30. I. Dag, D. Irk, and M. Tombul, “Least-squares finite element method for the advection-diffusion equation,” Applied Mathematics and Computation, vol. 173, no. 1, pp. 554–565, 2006.
31. B. Servan-Camas and F. T.-C. Tsai, “Lattice Boltzmann method with two relaxation times for advection-diffusion equation: third order analysis and stability analysis,” Advances in Water Resources, vol. 31, no. 8, pp. 1113–1126, 2008.
32. M. K. Kadalbajoo and P. Arora, “Taylor-Galerkin B-spline finite element method for the one-dimensional advection-diffusion equation,” Numerical Methods for Partial Differential Equations, vol. 26, no. 5, pp. 1206–1223, 2010.
33. M. Zerroukat, K. Djidjeli, and A. Charafi, “Explicit and implicit meshless methods for linear advection-diffusion-type partial differential equations,” International Journal for Numerical Methods in Engineering, vol. 48, no. 1, pp. 19–35, 2000.
34. J. Li, Y. Chen, and D. Pepper, “Radial basis function method for 1-D and 2-D groundwater contaminant transport modeling,” Computational Mechanics, vol. 32, no. 1-2, pp. 10–15, 2003.
35. R. W. MacCormack, “The effect of viscosity in hypervelocity impact cratering,” AIAA Paper 69-354, 1969.
36. M.-H. Tseng, “Kinematic wave computation using an efficient implicit method,” Journal of Hydroinformatics, vol. 12, no. 3, pp. 329–338, 2010.
37. R. J. Fennema and M. H. Chaudhry, “Explicit numerical schemes for unsteady free-surface flows with shocks,” Water Resources Research, vol. 22, no. 13, pp. 1923–1930, 1986.
38. R. J. Fennema and M. H. Chaudhry, “Simulation of one-dimensional dam-break flows,” Journal of Hydraulic Research, vol. 25, no. 1, pp. 41–51, 1987.
39. D. C. Dammuler, S. M. Bhallamudi, and M. H. Chaudhry, “Modeling of unsteady flow in curved channel,” Journal of Hydraulic Engineering, vol. 115, no. 11, pp. 1479–1495, 1989.
40. P. K. Mohapatra and M. H. Chaudhry, “Numerical solution of Boussinesq equations to simulate dam-break flows,” Journal of Hydraulic Engineering, vol. 130, no. 2, pp. 156–159, 2004.
41. C. M. Kazezyilmaz-Alhan, M. A. Medina Jr., and P. Rao, “On numerical modeling of overland flow,” Applied Mathematics and Computation, vol. 166, no. 3, pp. 724–740, 2005.
42. P. Rao and M. A. Medina Jr., “A multiple domain algorithm for modeling one-dimensional transient contaminant transport flows,” Applied Mathematics and Computation, vol. 167, no. 1, pp. 1–15, 2005.
43. P. Rao and M. A. Medina Jr., “A multiple domain algorithm for modeling two dimensional contaminant transport flows,” Applied Mathematics and Computation, vol. 174, no. 1, pp. 117–133, 2006.
44. P. Verma, K. S. H. Prasad, and C. S. P. Ojha, “MacCormack scheme based numerical solution of advection-dispersion equation,” ISH Journal of Hydraulic Engineering, vol. 12, no. 1, pp. 27–38, 2006.
45. R. Hixon and E. Turkel, “Compact implicit MacCormack-type schemes with high accuracy,” Journal of Computational Physics, vol. 158, no. 1, pp. 51–70, 2000.
46. S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” Journal of Computational Physics, vol. 103, no. 1, pp. 16–42, 1992.
47. T. K. Sengupta, Fundamentals of Computational Fluid Dynamics, Universities Press, Hyderabad, India, 2004.
48. X. Zhong, “High-order finite-difference schemes for numerical simulation of hypersonic boundary-layer transition,” Journal of Computational Physics, vol. 144, no. 2, pp. 662–709, 1998.
49. R. Szymkiewicz, Numerical Modeling in Open Channel Hydraulics, Springer, Amsterdam, The Netherlands, 2010.