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.

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.

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 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 .

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.