AN ACCURATE INTEGRAL-BASED SCHEME FOR DISPERSION EQUATION

Tung-Lin Tsai ¹, Jinn-Chuang Yang ¹ and Liang-Hsiung Huang ²

¹Department of Civil Engineering, National Chiao-Tung University, Hsinchu, Taiwan 30010, China

²Department of Civil Engineering, National Taiwan University, Taipei, Taiwan 10617, China

Abstract: This paper proposes an accurate integral-based scheme for solving dispersion equation. In the newly-proposed scheme the dispersion equation is integrated over a computational element using the quadratic polynomial interpolation function. Then elements are connected by the continuity of first derivative at boundary point of adjacent elements. The newly-proposed scheme is unconditionally stable and preserves a tridiagonal system of equations which can be solved efficiently by the Thomas algorithm. Using the method of fractional steps, the newly-proposed scheme, originally developed for one-dimensional problems, can be extended straightforward to multi-dimensional problems. To investigate the computational performances of the newly-proposed scheme five numerical examples are considered : (i) dispersion of Gaussian concentration distribution in one-dimensional uniform flow, (ii)one-dimensional viscous burger’s equation, (iii) pure advection of Gaussian concentration distribution in two-dimensional uniform flow, (iv) pure advection of Gaussian concentration distribution in two-dimensional rigid-body rotating flow and (v) three-dimensional diffusion in a shear flow. In comparison not only with the QUICKEST scheme, fully time-centred implicit QUICK scheme and fully time-centred implicit TCSD scheme for one-dimensional problem but also with the ADI-QUICK scheme, ADI-TCSD scheme and MOSQUITO scheme for two-dimensional problems, the newly-proposed scheme shows convincing computational performances.

Keywords: quadratic polynomial function, dispersion equation, integral-based scheme

1    INTRODUCTION

The dispersion equation is one of the governing equations used for modelling solute transport processes and water quality in rivers, lakes, oceans and groundwater. The conventional Crank-Nicolson second order central difference scheme [1] is a relatively simple and convenient way. However, due to central discretization of advection processes, it suffers from severe numerical oscillation under large Peclect number. This numerical oscillation can be vanished by using a first-order upwind-type finite difference scheme, but it induces large excessive numerical damping. Thus, in order to tackle the numerical oscillation problem and reduce excessive numerical damping, several high order upwind-type finite difference methods have been proposed, such as, QUICK (Quadratic Upstream Interpolation for Convective Kinematics) scheme [2], QUICKEST (QUICK with Estimated Stream Terms) scheme [2] and the third-order convection second-order diffusion (TCSD) scheme³. These numerical schemes are all explicit formulation originally. Later on, some implicit forms of modified QUICK [4,5] and TCSD [6] schemes were proposed.

In this paper an accurate integral-based numerical scheme for solving dispersion equation is proposed. In the newly-proposed scheme the dispersion equation is integrated over a computational element using the quadratic polynomial interpolation function. Then elements are connected by the continuity of first derivative between boundary point of adjacent elements. By using the method of fractional steps and the technique of linearization, the newly-proposed scheme, originally developed for one-dimensional linear problems, can be extended straightforward to multi-dimensional and nonlinear problems. Five numerical examples, including (i) dispersion of Gaussian concentration distribution in one-dimensional uniform flow, (ii)one-dimensional viscous burger’s equation, (iii) pure advection of Gaussian concentration distribution in two-dimensional uniform flow, (iv) pure advection of Gaussian concentration distribution in two-dimensional rigid-body rotating flow and (v) three-dimensional diffusion in a shear flow, are used to investigate the computational performances of the newly-proposed scheme.

2    FRAMEWORK OF THE NEWLY-PROPOSED SCHEME

2.1    Discretization

The transient one-dimensional dispersion equation can be written as

                           (1)

where the scalar function  may represent, for example, temperature or concentration at position x and time t in a fluid moving with a speed  and  diffusion coefficient . First, equation (1) is integrated over an interval from position to position  and is expressed as

                     (2)

By defining  as

                               (3)

and adopting quadratic polynomial interpolation function in the interval , the first derivatives at point and  can be expressed as

                          (4)

                          (5)

under a uniform grid space ,i.e., ,as shown in Figure 1. Substituting equations (3), (4) and (5) into equation (2), one obtains

  (6)

Using the Crank-Nicholson technique to equation (6) for time discretization, the average quantity over the interval  at  time step can be expressed as

           (7)

where

   ,

, ,

with the Courant number 、 and Diffusion number 、 . Assuming the first derivative of dispersion equation is continuous at any grid point, the continuity of the solution at point and  time step is imposed by satisfying the following condition

                     (8)

By substituting equation (7) into (8), the discretization of dispersion equation becomes

  (9)

where

      , ,

      , , , ,

If the diffusion coefficient  are constant, equation (9) can be reduced to

             (10)

where

  , ,

, ,   

with the Diffusion number .

2.2    Solution procedure and stability analysis

From equation (9) the newly-proposed scheme would result in a tridiagonal system of algebraic equations which can be solved efficiently by the Thomas algorithm [7]. The boundary condition of Dirichlet type can be applied to the newly-proposed scheme directly, whereas for the Neumann boundary condition the discretisation technique of finite difference method may be introduced to remain the tridiagonal system. If the diffusion coefficient is constant, one can directly apply equation (10) to solve the dispersion equation without computing the average

quantity of every element. For the non-constant diffusion coefficient the solution procedure of the newly-proposed scheme is depicted as follows: 1.Specify initial value, , at every grid point and compute initial average value, , for every interval from the initial condition. 2. Use equation (9) along with the boundary conditions to solve for the unknown, , at every grid point in space at next time step. 3. Use equation (7) to compute explicitly the unknown at every space interval at next time step. 4. Repeat steps 3 and 4 to the end of time of simulation. In this study the von Neumann stability analysis is conducted by assuming that the velocity and diffusion coefficient are constant and positive. The von Neumann stability analysis shows that the newly-proposed scheme is an unconditionally stable scheme.

2.3    Extension to multi-dimensional problems

The newly-proposed scheme is originally developed for one-dimensional linear dispersion equation. Using the method of fractional steps [8], the newly-proposed scheme can be extended straightforward to multi-dimensional dispersion problems without complication and difficulty. In addition, using the technique of linearization, the nonlinear problems can also be solved interatively by the newly-proposed scheme.

3    NUMERICAL EXAMPLES

3.1    One-dimensional problem

Calculation of dispersion

To investigate the computational performance of the newly-proposed scheme, the dispersion of a Gaussian concentration distribution is considered with constant velocity u = 0.3 m s^-1 and diffusion coefficient  = 1.0 m² s^-1. The computed results are shown in Figure 2 for the simulation time of 20000 s with time step size = 100 s and grid space =100 m.  From Figure 2, one can find that the newly-proposed scheme induces the least numerical oscillation and diffusion in comparison with the QUICKEST scheme, fully time-centred implicit QUICK scheme and fully time-centred implicit TCSD scheme.

3.2    Two-dimensional problem

Calculation of pure advection in uniform flow

The Gaussian concentration distribution with the peak value of 10 and a standard deviation of 220 m is advected for 5000 s at constant velocity = 0.5 ms^-1 and = 0.5 ms^-1. The grid size of 100 m × 100 m and time step of 100 s are used. The computed results in terms of maximum, minimum values and RMS error by using the newly-proposed scheme, ADI-QUICK scheme, ADI-TCSD scheme and a modified second-order QUICKEST scheme (MOSQUITO)[9] are displayed in Table 1. From Table 1, one can find that the newly-proposed yields the best results among those schemes considered. The ADI-TCSD scheme induces the largest numerical oscillation and RMS error.

3.3    Three-dimensional problem

Calculation of diffusion in shear flow

In order to investigate the capability of the newly-proposed scheme for solving three-dimensional problem with complex flow field, a diffusion in the shear flow is considered. The velocity shear in the diffusion of a patch of passive contaminant from an instantaneous source plays an important role in natural streams such as ocean, lake and estuaries. The governing equation for shear diffusion can be expressed by

              (11)

where x, y and z is the coordinate system; represents the concentration of contaminant;  is the mean velocity in the x direction; and represent the horizontal and vertical shears. In addition, and  are the eddy diffusivities in the x, y and z directions, respectively. The analytical solution for an instantaneous point source of mass M released at x = y = z = 0 was obtained by Carter and Okubo [10] as following

  (12)

where

                 (13)

To allow numerical solution having an initial peak concentration of unity, calculation begins at time . Thus, the point source of mass can be expressed as

                 (14)

The following parameters are used in this numerical example: 1000 s, 0.5ms^-1, 0.0003 s^-1,  8.0 m²s^-1, time step size 100s and grid space  100 m. Figure 3 depicts the contour plots as simulated by the newly-proposed scheme and analytical solution at time = 3000 s and 5000 s on the plan , respectively. It is obvious that the newly-proposed scheme yields computational results which are almost in agree with the exact solution.

4    CONCLUSIONS

In this paper an accurate integral-based numerical scheme for solving dispersion equation is proposed. The newly-proposed scheme is unconditionally stable and preserves a tridiagonal system of equations which can be solved efficiently by the Thomas algorithm. By using the method of fractional steps and the technique of linearization, the newly-proposed scheme that originally developed for one-dimensional linear problem can be extended straightforward to multi-dimensional and nonlinear problems. Some one-, two- and three- dimensional numerical examples are used to investigate the computational performance of the newly-proposed scheme. In comparison with the QUICKEST scheme, fully time-centred implicit QUICK scheme, fully time-centred implicit TCSD scheme, ADI-QUICK scheme, ADI-TCSD scheme, MOSQUITO scheme and exact solutions, the newly-proposed scheme yields convincing simulated results.

References

[1] J. Crank and P. Nicolson, ‘A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type’, Porc. Camb. Phi. Soc., 43, 50-67 (1947).

[2] B. P. Leonard, ‘A stable and accurate convective modelling procedure based on quadratic upstream interpolation’, Comput. Methods Appl. Mech. Eng., 19, 59-98 (1979).

[3] D. Bradley and M. Missaghi, ‘ A Taylor-series approach to numerical accuracy and a third-order scheme for strong convective flows’, Comput. Methods Appl. Mech. Eng., 69, 133-151 (1988).

[4] B. P. Leonard and B. J. Noye, ‘Second and third order two level implicit FDM’s for unsteady one-dimensional convection-diffusion’, Computational Techniques and Applications: CTAC-89, Hemisphere, New York, 1990, pp. 311-317.

[5] Y, Chen and R.A. Falconer, ‘Advection-diffusion modelling using the modified QUICK scheme’, Int. j. numer. methods fluids, 15, 1171-1196 (1992).

[6] Y. Chen and R. A. Falconer, ‘ Modified forms of the third-order convection second-order diffusion scheme for the advection-diffusion equation’, Advances in Water Resources, 17, 147-170 (1994).

[7] L. H. Thomas, Elliptic Problems in Linear Difference Equations over a Network, Waston Scientific Computing Laboratory, Columbia, NY, 1949.

[8] N. N. Yanenko, The Method of Fractional Steps: The Solution of Problems of  Mathematical Physics in Several Variables, Springer-Verlag, NY, 1971.

[9] A. Balzano, ‘ MOSQUITO : An efficient finite difference scheme for numerical simulation of 2D advection.’ Int. j. numer. methods fluids, 31, 481-496 (1999).

[10] H. H. Cater, and A. Okubo, A study of the physical process of movement and dispersion in Cape Kennedy area. Final Rep. U.S. At. Energy Comm. Contract No. AT (30-1), 1965.

   Table1    Performances of various schemes in two-dimensional pure advection test

Fig.1    Sketch of grid representation of the newly-proposed scheme

Fig.2    Comparison of various schemes for dispersion of Gaussian concentration distribution

    Fig.3    Comparison of contour plots of diffusion in shear flow on plane z = 0 at time = 3000s,
5000s by using newly-proposed scheme and exact solution