A finite Volume Runge-Kutta Ellam Method for the Solution of Advection Diffusion Equations

ةصلاخ :  لاقتنلإا تلاداعم لحل ةزيمملا تاينحنملا ىلع دمتعت ةدودحملا موجحلاب ةيددع ةقيرط ميمصتب ثحبلا اذه صتخي ثولملا لاقتنإ لثمت يتلا راشتنلإاو ةيماسملا طاسولأا يف تا . تاينحنملل ةيناثلا ةبترلا نم اتوكو جنور بيرقت ةقيرطلا هذه مدختست قرط راطإ يف ةزيمملا " ملايإ "  .(Eulerian Lagrangian Localized Adjoint Method , ELLAM) ةقيرطلا زيمتت روص ىلإ ةمكاحلا تلاداعملا ليوحتو ةلتكلا ىلع ظافحلاب اهططخم يف ةحرتقملا ىتح ةقيقد ةيددع لولح ىلإ يدؤت امك ةلثامتم ة ةريبك ةينمز تاوطخ مادختسإ لظ يف . قرط ةدع ةنراقمل ةيددعلا براجتلا ضعب ضرعب موقن  حيضوتل يطمن لاثم لحل ةفورعم ةحرتقملا ةقيرطلا ةءافك .


Introduction
Many problems arise in the numerical simulation of subsurface contaminant transport and remediation within porous media.The model equations, which describe such processes, are advection-diffusion equations which are known to present many numerical difficulties especially when advection dominates the physical process.Standard centered difference methods and Galerkin finite element methods (FEM) then generate solutions which exhibit non-physical spurious oscillations.Standard upwinding methods on the other hand usually eliminate these types of oscillations in their solutions.However they suffer from excessive numerical diffusion which smears out sharp fronts of the solutions where important chemistry and physics take place.
Many specialized methods have been developed which aim at resolving the difficulties mentioned when applied to both linear and nonlinear problems.A large class of these methods, usually referred to as Eulerian methods, use some form of upstream weighting and improved techniques in their formulations over fixed spatial grids.This class includes the Petrov-Galerkin FEM methods (Bouloutas and Celia, 1991;Westerink and Shea, 1989), the streamline diffusion methods (Hughes and Mallet, 1986;Johnson, 1987), the flux corrected transport methods (Boris and Book, 1997;Tóth and Odstrčil, 1996), and the high resolution methods from fluid dynamics which include the essentially non-oscillatory method (ENO) and the weighted ENO method (Shu 1999), and the Godunov methods (Dawson, 1991;Van Leer, 1984), as well as many other methods.These methods which are characterized by ease of formulation and implementation, succeed in suppressing the artificial oscillations which are present in standard methods.However, their solutions tend to be dominated by time truncation errors.Moreover they impose restrictions on the size of the time step taken for reasons of stability and accuracy.
A second class of methods, usually referred to as characteristic methods makes use of the dual nature of the governing equation which includes hyperbolic as well as parabolic components (Douglas and Russell, 1982;Pironneau, 1982).These methods incorporate Eulerian grids with Lagrangian tracking along the characteristics to treat the advective part of the equation.This treatment allows for larger time steps to be used in the simulation.Moreover, it significantly reduces the time truncation errors when compared to Eulerian methods.However, these methods have difficulty in conserving mass and in treating general boundary conditions.
The Eulerian Lagrangian localized adjoint method was developed by Celia, Russell, Herrera, and Ewing (1990) as an improved extension of characteristic methods which maintains their advantages and enhances their performance by conserving mass and treating general boundary conditions naturally in its formulation.This first ELLAM formulation (Celia et al, 1990) was a finite element formulation for one-dimensional constant-coefficient advection diffusion equations.The strong potential that this formulation has shown, led to the development of formulations for variable coefficient equations (Russell and Trujillo, 1990;Wang et al, 1992), for non-linear equations (Dahle et al, 1995), as well as finite volume formulations (Healy and Russell, 1993).
Most of the ELLAM formulations developed use a backward Euler approximation in time along the characteristics due to its simplicity and stability.These formulations are therefore only first order accurate in time.The authors have recently developed second order in-time ELLAM methods for advection-diffusion equations (Al-Lawatia et al, 1999).The derived schemes were shown to have a higher order in-time convergence rate compared to the backward Euler ELLAM methods.
In general, ELLAM methods generate regularly structured systems which are symmetric and positive definite and thus can be easily solved numerically.A principal drawback, however, to the increased use of ELLAM methods has been the increased effort required to implement the method within existing simulators, since it is somewhat more involved than several of the other methods.In this paper we develop a finite volume ELLAM method for the solution of variable coefficient advection-diffusion equations in one spatial dimension which uses a second order Runge-Kutta approximation for the characteristics.This method uses finite volume test functions in the spacetime domain defined by the characteristics which permits implementation with little code change in many existing simulators, whether they be legacy codes or more modern modularized versions.Numerical experiments are presented to illustrate the performance of the method developed and to compare it with many widely used and well-perceived methods, such as the streamline diffusion method, monotone schemes, essentially non-oscillatory methods, and flux corrected schemes.

Development of the Runge-Kutta ELLAM Method
The model problem we consider is the following one-dimensional advection-diffusion equation and initial condition where V is the velocity field and is the positive diffusion coefficient.To close the system, we assume ( b)-periodic boundary conditions.

Partition and Characteristic Curves
Let I and N be two positive integers.Following Healy and Russell (1993), we define a partition of the space-time T b a × of equation (1) as follows: (2) .... : The partition x δ divides the spatial domain into grid blocks (finite volumes) of size .We denote the center of block , by which represents a grid point in our discretization.
Multiplying equation ( 1) by a test function that vanishes outside and integrating by parts we obtain a weak form of equation ( 1 where which takes into account the fact that is discontinuous in time at time t .We note from the periodicity of formulation that the third integral on the lefthand side of equation (3) vanishes.
In the ELLAM framework (Celia et al, 1990), the test functions in equation (3) are selected to satisfy (exactly or approximately within the tolerance of the accuracy desired) the homogeneous equation of the hyperbolic part of the adjoint equation of equation ( 1) which reflects the Lagrangian nature of equation ( 1).This implies that the test functions should be chosen to be constant along the characteristics.These characteristics are given by solutions of initial value problems of the ordinary differential equation ). , ( t y V dt dy = (5) Solving this equation analytically for a generic velocity field, however, is not possible, and so we must consider approximate solutions.In the formulation of our scheme, we approximate the characteristics by a second order Runge-Kutta method and define the characteristic curve emanating from a point using Heun's method: where θ determines the time position along the characteristic.

The Finite Volume Runge-Kutta ELLAM Scheme
In the finite volume Runge-Kutta ELLAM scheme, we define the test function associated with the grid point by and extend to be constant back along the approximate characteristics (6) into the space-time The Runge-Kutta ELLAM scheme can be formulated by evaluating the space-time integrals in equation ( 3) along the approximate characteristics.We evaluate the second (source) term on the right hand side of equation ( 3) by a trapezoidal quadrature where is the foot of the approximate characteristic emanating from ( and is the truncation error due to the application of the trapezoidal rule.
) ,  ; ( We evaluate the second (i.e., diffusion) term on the left hand side of equation ( 3) in a similar manner as was done with the source term but noting that w is a Dirac-delta function at the two end points and , where we used the identity in the third equality.
The finite volume Runge-Kutta ELLAM can then be formulated from the weak form (3) using the periodic boundary conditions to drop the third term on the left hand side, dropping the fourth term on that same side by our choice of test functions (main idea of ELLAM), and by substituting equations ( 8) and ( 9), dropping both error terms and .Finally, we choose as trial functions which are piecewise linear in space on adjacent grid points at time t and obtain the following Runge Kutta ELLAM scheme The first integral on the left hand side of equation ( 10) is a standard integral which can be computed as follows for a non-uniform partition x δ and as for a uniform one.The second integral on the right hand side of equation ( 10), which represents a source/sink integral, can be evaluated in a similar manner.
The remaining integrals in equation (10) which are defined at time t require more careful evaluation.The tracking approach we use to evaluate these terms is forward tracking, which avoids the grid deformation problems associated with backward tracking schemes and is more feasible for multiple dimensions.In this approach numerical integration rules are applied at time t where the solution (and the associated weight) at each quadrature point is forward tracked while the value of the test function is determined at the location to which the quadrature point is tracked at time .However, using the test functions given by equation ( 7) presents mass conservation difficulty when the Courant number is close to an integer, since forward tracked quadrature points could come very close to end points of the cell blocks which could result in mass not being accurately distributed among the nodes.Healy and Russell (1993) suggest using alternative test functions which approximate the functions given in equation ( 7) and try to distribute mass more accurately while maintaining local mass conservation.We use their approach and approximate the first integral on the leftside of equation ( 10) according to where the test function introduced is defined by

Numerical Experiments
In this section we present numerical results to illustrate the performance of the finite volume Runge-Kutta ELLAM (RK-ELLAM) method developed in this paper and compare it to other well known methods.

Convergence Rates of the Runge-Kutta ELLAM Scheme
We observe numerically the order of convergence of the finite volume RK-ELLAM method developed for equation (1).In order to make comparisons with earlier results in (Al-Lawatia et al, 1999), we choose as a test problem the standard transport of a Gaussian hill with initial configuration given by where the center is and the spread The spatial domain is [ and we simulate over a time interval of [0,1].The velocity field considered is V and the diffusion coefficient is set to .The right hand side in equation ( 1) is generated from this data and the analytical solution (see section 5.1 of Al-Lawatia et al, (1999) for details).
In order to obtain the order of convergence in space (respectively, time) in this experiment, we perform runs varying the mesh size (respectively, x ∆ t ∆ ) with the remaining mesh size being fixed with small value, so that the contribution to the residual error due to that variable is negligible.We then use a linear regression to fit the data and obtain the order of convergence in the chosen variable whose size was varied.The process is then repeated for the remaining variable.Table 1.contains the norms of the errors generated for our runs as well as the estimated order of convergence and corroborates that the scheme is second order in time.

Comparison of the Runge-Kutta ELLAM Method With Other Schemes
We perform numerical experiments which present a comparison of the finite volume Runge-Kutta method developed in this paper with some well known and widely used methods, including the Galerkin (GAL), Quadratic (QPG), and Cubic (CPG) Petrov-Galerkin finite element methods, the Streamline diffusion finite element method (SDM), the monotone upstream-centered scheme for conservation laws (MUSCL), the Minmod scheme, the essentially non-oscillatory method (ENO) and weighted ENO (WENO), and flux-corrected transport (FCT) methods.See Al-Lawatia et al, (1999) and the references cited there for a discussion of most of these methods, their implementation, and the advantages and disadvantages of their use for numerically solving advection-diffusion problems.The QPG and CPG methods (Bouloutas and Celia,1991;Westerink and Shea, 1989) incorporate some upwinding in the test space as compared to the standard finite element method (GAL).Using the same principle, the SDM incorporates some upwinding by adding multiple of the linearized hyperbolic operator to the test functions in the space-time finite element formulation of equation ( 1) (Hughes and Mallet 1986;Johnson 1987).This results in numerical diffusion being added only in the direction of the streamlines.The amount of diffusion added depends on the value taken for a constant which appears in the formulation.There is no clear optimal choice for this constant and is heavily problem dependent.Following our reasoning and explanations in (Al-Lawatia et al,1999), we consider three different values for C , namely 1.0, 0.1 and 0.0001.The MUSCL and Minmod schemes are Godunov methods which are known to work well for hyperbolic conservation laws, and are extended for the advection-diffusion equations by operator splitting of the equation.In our experiment we consider the formulation given by (Dawson, 1991;van Leer, 1984).In a similar manner, the ENO and WENO which are high-order methods can be extend to solve the advection-diffusion equations (Shu, 1999).The ENO formulation of order k uses a polynomial interpolation over the locally smoothest stencil from k-1 candidates.This treatment provides a non-oscillatory interpolation which accurately resolves sharp fronts of the solution.The WENO formulation uses a convex combination of all candidate stencils.This allows for higher order (2k-1) approximation of smooth data while maintaining the advantages of ENO.In our experiments we consider ENO and WENO of order k=3 while the Roe flux is used to guarantee upwinding.The last method we consider is the Flux Corrected Transport scheme (FCT) which adds a high-order (anti-diffusive) term to a low order solution and a limiter which controls the amount of anti-diffusion added thus avoiding the formation of new minima or maxima.In these experiments we consider two FCT formulations, The ETB-FCT which is due to Boris and Book (1997) and the YD-FCT which is due to Odstrčil (Tóth and Odstrčil, 1996).

C
Table 3: The performance of the galerkin and petrov-galerkin FEM method.

Summary
In this paper we develop a finite volume characteristic method for the solution of variable coefficient advection diffusion equations in one-space dimension.This method uses a second order Runge-Kutta approximation for the characteristics within the framework of the Eulerian Lagrangian localized adjoint method.The derived scheme conserves mass, fully utilizes the transient behavior of the governing equation, and generates accurate numerical approximations even for large values of the Courant number, and thus it permits large time steps in thesimulation.Numerical experiments are presented which show that, in the context of linear advection-diffusion equations, the finite volume Runge-Kutta ELLAM method outperforms many well perceived and widely used methods, in a standard test example, including the Galerkin and Petrov-Galerkin FEM, the streamline diffusion method, the Flux-Corrected Transport methods, and the high resolutions methods MUSCL, Minmod, Essentially non-oscillatory method (ENO) and weighted ENO.The reader is referred to (Al-Lawatia et al, 1999) for other experiments which compare two finite element formulations of ELLAM to a wide range of well perceived and widely used methods.

Table 1 :
Spatial and temporal convergence rates of the RK-ELLAM method.

Table 2 :
The performance of the RK-ELLAM method.

Table 4 :
errors for the various methods described in the previous paragraph and with groupings determined according to their common features.In the example runs, we use a base spatial grid size of =1/60 which is needed to resolve the analytical solution.In the RK-ELLAM simulation we use a time step of =1/10 which gives a Courant number of 7.08.The RK-ELLAM solution is presented in Figure1, which matches the analytical solution with no noticeable artifacts.The norm of the error is given in Table2.All other methods have a restriction on the methods to insure that the Courant number is less than one.We then vary t ∆ and x ∆ separately until we find a solution comparable to that of RK-ELLAM.For compactness of presentation, we only display representative results in Tables 2-7 for all methods The performance of the streamline diffusion method with C=1.0,0.1, 0.0001.

Table 5 :
The performance of the MUSCL and minmod methods.

Table 6 :
The performance of the essentially and weighted essentially non-oscillatory schemes.These tables include the errors in the norm, and the CPU time (measured on a Sun Workstation) used in the simulations.In addition we present in Figures2-6two representative plots for each of the groups of methods described above.One plot uses

Table 7 :
The performance of the flux-corrected transport methods ETB-FCT and YD-FCT.