Assessment of a non-traditional operator split algorithm for simulation of reactive transport

https://doi.org/10.1016/j.matcom.2005.03.019Get rights and content

Abstract

A non-traditional operator split (OS) scheme for the solution of the advection-diffusion-reaction (ADR) equation is proposed. The scheme is implemented with the recently published central scheme [A. Kurganov, E. Tadmor, New high-resolution central schemes for non-linear conservation laws and convection-diffusion equations, J. Comput. Phys. 160 (2000) 241–282] to accurately simulate advection-reaction processes. The governing partial differential equation (PDE) is split into two PDEs, which are solved sequentially within each time step. Unlike traditional methods, the proposed scheme provides a very efficient method to solve the ADR equation for any value of the grid-Péclet number. An analytical mass balance error analysis shows that the proposed non-traditional scheme incurs a splitting error, which behaves differently to the splitting error incurred in traditional OS schemes. Numerical results are presented to illustrate the robustness of the proposed scheme.

Introduction

The operator split (OS) technique provides a convenient means of solving the advection-diffusion-reaction (ADR) equation:ct=D2cx2vcx+f(c),where c [ML3] is the solute concentration, D [L2T1] the diffusion coefficient, v[LT1] the advective velocity, x [L] the spatial variable and t [T] is time. The function f(c) represents kinetic processes within the system; for some applications, non-linearities in the system are confined to f(c). For the purposes of this manuscript, we will interpret f(c) as a reaction.

The standard two-step OS scheme used to simulate (1) involves splitting the reaction part from the equation, yielding a simpler linear PDE and a linear or non-linear system of ordinary differential equations (ODE) [5], [9], [12], [20], [24], [25], [28], [30]:ct=D2cx2vcx,dcdt=f(c).Traditionally, the linear advection-diffusion equation (2) is solved using a standard finite difference or finite element approach [20], [30], while the reaction equation (3) is solved using an ODE solution technique that is appropriate for the particular form of f(c)[20], [30]. We refer to this scheme as the DA-R scheme.

The use of OS schemes have become widespread because of the flexibility offered by choosing different solution techniques to suit the various terms in the conservation equation [24]. OS schemes are also advantageous because the decoupling of reaction and transport discretisations permits different reaction types to be simulated with minimal recoding [13], [14]. Obtaining numerical solutions using an OS technique can increase the computational efficiency compared to solving the full system and also provides opportunities for parallel computation [16], [19]. OS schemes are particularly attractive for coupled multi-species problems or problems with non-linear kinetics [11], [18]. Possible extensions of the traditional OS technique, such as the alternating OS scheme [20], [25], [30], the iterative OS scheme [18] and specialised schemes that combine OS with LU factorisation [9] have been proposed. This work, however, focuses on a sequential OS scheme as this is the most widely used approach [8], [14], [32], [35]. Previous investigations have shown that solving the ADR equation with the traditional two-step OS scheme can introduce an error, which is minimised by choosing Δtsufficiently small [20], [25], [30].

Traditional OS schemes can provide several advantages over solving the full conservation equation [30]. However, if the transport part of the equation is solved with standard Eulerian techniques, then the overall accuracy is limited by the well-known grid Péclet number restriction. The grid Péclet number restriction for mass-lumped linear finite element or finite difference solutions of the advection-diffusion equation requires Pegrid=vΔxD2in order to prevent spurious oscillations and excessive numerical diffusion [8]. For advection-dominated problems, the grid Péclet number restriction can be quite severe, requiring excessively fine grids to prevent numerical errors.

In addition to the traditional DA-R scheme (2), (3), previous analyses have used a three-step splitting scheme, whereby the diffusion, advection and reaction steps are all separated (D-A-R) [13], [14], [32], [35]:ct=D2cx2,ct=vcx,dcdt=f(c).The D-A-R procedure has been successfully used to simulate the ADR equation, where the advection term is typically simulated using particle tracking or front tracking techniques [13], [14], [32], [35]. These methods have been invoked to help overcome the grid Péclet number restriction of traditional Eulerian approaches. The advantage of the D-A-R method is that the original PDE is split into two linear PDEs and a linear or non-linear ODE system. Each of these steps can be solved by a solution technique appropriate to the nature of each term.

Overcoming the traditional difficulties associated with solving the ADR equation under advection-dominated conditions has been an active area of research. Investigations of alternative Eulerian algorithms that minimise numerical error for advection-dominated conditions have yielded a wide range of improved schemes. A survey of the developments of alternative difference schemes can be found in Yeh [6], Leonard [34] and Christon et al. [33]; similar developments based on weighted residual strategies are surveyed by Yeh [6] and Donea and Huerta [2].

Recent advances in solving hyperbolic conservation laws with central schemes and flux limiters provide high accuracy alternatives to standard Eulerian solutions to the advective transport equation [33]. In this work, the recent central scheme proposed by Kurganov and Tadmor (KT) [21] is used to solve the advection-reaction part of the ADR equation. The KT discretisation permits a semi-discrete form, which is explicitly integrated with a Runge–Kutta method. It is therefore convenient to simply incorporate any reaction terms into this integration process thereby temporally integrating both the advection and reaction terms simultaneously. Following this argument, we propose a “non-traditional” OS scheme, whereby the ADR equation is split:ct=D2cx2,ct=vcx+f(c).This non-traditional OS scheme (D-AR) at first appears more complicated than the traditional OS schemes as the single linear or non-linear PDE has been split into two PDEs: the first is a linear diffusion equation and the second is a linear or non-linear advection-reaction equation. The reason for splitting a single PDE into two other PDEs is that the KT algorithm requires explicit temporal integration [21]. In general, temporal integration of the diffusion term cannot be completed explicitly because of the well-known stability criteria due to the stiffness of the discretised form of the diffusion operator [1], [8]. Therefore, to maintain generality, the diffusion equation is to be solved implicitly using a standard finite difference or finite element approach.

The impetus for splitting the ADR equation according to (7), (8) originated through simulating cell migration problems. For models of cell migration by haptotaxis [23], [26] and chemotaxis [22], the advection term is non-linear. In these cases, the advective velocity is proportional to the spatial gradient of a signalling (chemotaxis) or adhesive (haptotaxis) substance. The advective velocity can be temporally and spatially variable and is determined as a part of the solution. This complication results in a non-linear advection term [22]. Under these conditions, traditional splitting (2), (3) does not convert the non-linear PDE into a linear transport PDE and linear or non-linear ODE pair. This is an important point that distinguishes the application of OS strategies to non-linear advection problems from classical linear advection problems such as those encountered in contaminant transport analyses within aquifer systems [6], [8].

Previous work by Marchant et al. [23] investigated the influence of adding diffusion to the pure haptotactic tumour invasion problem described by Perumpanani et al. [26]. Marchant et al.’s numerical simulations were limited by stability constraints owing to the explicit temporal integration of the diffusion equation within the KT algorithm. Therefore, only small amounts of diffusion were added [23]. Recent work by Landman et al. [22] successfully used an OS scheme analogous to the proposed non-traditional splitting (7), (8) to perform numerical simulations of a combined diffusion and chemotaxis cell invasion problem. Landman et al. simulated any combination of diffusion or chemotaxis without any stability constraints associated with the diffusive migration [22]. This previous work successfully investigated the interactions of diffusive and chemotactic migration, and all simulations were performed with a non-traditional splitting scheme similar to that which is to be analysed here.

The objective of this work is to explore the robustness of D-AR splitting scheme for solving the ADR equation. The explicit KT central scheme is used to solve the advection-reaction part of the conservation equation. The diffusion term is split and solved implicitly. Because the D-AR OS scheme is different from all standard OS approaches, the convergence behaviour of the proposed OS scheme is investigated. In addition, the mass balance error (MBE) incurred by the proposed D-AR scheme is investigated analytically and compared to the MBE for standard DA-R and D-A-R approaches. This analytical error analysis is particularly important as previous analytical analyses of OS errors have been restricted to standard sequential and alternating forms of traditional two-step OS schemes [20], [25], [30].

Section snippets

Implicit finite difference solution of the diffusion equation

The first step in the D-AR OS scheme is to solve (7). The diffusion equation is solved using a standard implicit finite difference approach [1], [3]. The domain is uniformly discretised with nodes spaced Δxapart. The nodes are indexed i=1,,N. Temporal integration can be performed with a time-weighted scheme, which yields the discretised form of the diffusion equation at a central node:cij+1cijΔt=ωDci+1j+12cij+1+ci1j+1(Δx)2+(1ω)Dci+1j2cij+ci1j(Δx)2,i=2,,N1,j=0,,T,where the subscript i

Discussion and conclusion

This work has extended the applicability of the KT central scheme to be used in conjunction with an OS technique to simulate ADR systems. The algorithm presented here is likely to be of interest to researchers analysing reactive transport problems, which are ubiquitous throughout all applied science disciplines. The KT central scheme is typically implemented with an explicit method for the temporal integration of the discretised equations. Therefore, to develop a general algorithm for the

Acknowledgements

Matthew Simpson and Kerry Landman are supported by a National Health and Medical Research Council project grant ID237144. Kerry Landman also appreciates the support of the Particulate Fluids Processing Center, an Australian Research Council Special Research Center. In addition, we thank D. Andrew Barry for his thorough review comments and Anjani Kumar and for his initial comments.

References (35)

Cited by (38)

  • Local meshless method for PDEs arising from models of wound healing

    2017, Applied Mathematical Modelling
    Citation Excerpt :

    In general, the discretization of convection term in convection-dominated IBVPs must be tackled with care [20]. Conventional discretization approaches are frequently unsuitable for such types of problems [19] due to excessive numerical smearing and nonphysical oscillations near sharp gradients if not complemented with some stabilization procedure. Additional difficulties include strong couplings and nonlinearities of advection-dominated systems of PDEs [21].

View all citing articles on Scopus
View full text