A Domain-Boundary Integral Treatment of Transient Scalar Transport with Memory

It is well known that the boundary element method (BEM) is capable of converting a boundaryvalue equation into its discrete analog by a judicious application of the Green’s identity and complementary equation. However, for many challenging problems, the fundamental solution is either not available in a cheaply computable form or does not exist at all. Even when the fundamental solution does exist, it appears in a form that is highly non-local which inadvertently leads to a system of equations with a fully populated matrix. In this paper, fundamental solution of an auxiliary form of a governing partial differential equation coupled with the Green identity is used to discretize and localize an integro-partial differential transport equation by conversion into a boundary-domain form amenable to a hybrid boundary integral numerical formulation. It is observed that the numerical technique applied herein is able to accurately represent numerical and closed form solutions available in literature.


Introduction
The governing equation for a convective-dispersive transport phenomenon is given by: where c is concentration; v is the flow velocity; D is the dispersive tensor; t is the time variable and R is the rate of solute mass exchange.Equation (1) yields satisfactory results in controlled environments in the absence of media heterogeneities.However, it exhibits certain limitations in modeling real life phenomenon.Its deviation from Fickian theory, results in an infinite speed of propagation of the scalar field (Cattaneo [1] Curtin and Pipkin1 [2], Ferreira and Oliveira [3], Hristov [4]).In order to overcome this unphysical property, Fick's law is replaced by a flux accompanied by a memory term ( ) , J x t τ = + .It was proposed (Cattaneo's [1]) that the flux at a certain point x and at time t should be a consequence of scalar variation at some point x but after an elapsed time.Based on this, Ficks equation is modified to read: where τ is a delay parameter.Applying a first order Taylor series expansion to the flux we obtain the so-called Cattaneo's flux term.
( ) Without any loss in generality, Equation ( 1) can be recast to read where ( ) ( ) , , f c q x t are the reaction and source/sink terms.The source/sink term accounts for the injection or extraction of mass into or out of the system., D v are the diffusion coefficient and velocity variables.The term inside the integral term accounts for the memory of the system.Equation ( 4) is a one-dimensional integro-partial differential equation which describes the transport of quantities such as heat, vorticity, energy and mass and is very important in many physical systems especially those involving the environment.It comes with the following constraints.
Fundamental solutions generated from "free-space" Greens functions are usually applied in boundary integral analysis of transport equations.Efforts to make them more applicable to boundary element method (BEM) formulation, especially for non-homogeneous material, can be found in (Shaw [5]) and (Shaw and Manolis [6]).These have not only met with mixed results but have also been found to be inadequate for a good number of problems of engineering interest.
For example, adopting the classical boundary element approach to viscoelastic materials or to cases which incorporate the memory term requires fundamental solutions that are yet not well known except for the simple Maxwell material whose analytical solution is available (Gaul and Schanz [7]).Numerical attempts to deal with this include that of Schanz and Antes [8].They studied a BEM viscoelastic formulation based on the so called convolution integral method proposed by Lubich [9] in 1988.In this work, a quadrature formula whose weights are determined by the Laplace transform of the fundamental solution is applied to evaluate a convolution integral and the entire formulation numerically handled by a multi-step procedure.
The hybrid boundary-integral-finite element procedure adopted herein employs the Green's function of the Laplace operator, and retains the time marching feature of classical BEM, but its spatial discretization enlists a local support feature that is similar to that of the finite element.This guarantees some attractive outcomes; it facilitates the hybridization process especially for schemes that are element or nodal based (Onyejekwe [10] [11]) and establishes a hybrid computational method which exploits the strengths of its composite parts.This is the key motivation for this study.

Numerical Formulation
A typical boundary integral discretization requires both the complementary differential equation and the freespace Green's function.

(
) ( ) where i x is the source node and x is the field node and by the same token both t and t represent the source and field nodes in the temporal coordinates.The flexibility offered by a hybridization procedure allows the integral kernel in equation ( 4) to be evaluated at each node.We consider ( ) and We adopt a composite weighted trapezoidal rule in the temporal coordinate for the memory term. ( ) Equation ( 4) is finally converted to its integral analog by the Green's second identity to yield: where , , , α β γ η are given and I is the "memory" or the hereditary' term; α can be a constant or a function of the dependent variable, This hybrid formulation unlike a typical boundary integral formulation requires that the computational domain be discretized into elements over which the scalar variable will assume some functional distributions.This is represented for the dependent variable as where 2 m = represents the current time level, while 1 m = denotes previous time level.

( )
Ω are interpolating functions in space and time respectively.The element equation can now be written as: Equation ( 14) is put in a compact matrix form

Numerical Test Examples
Example 1 Consider the following special case of Equation ( 4) with reaction term of the nonlinear variety (Khuri and Sayfy [12]), where ( ) ( ) ( ] , 0,1 0,1 x t ∈ × subject to the initial condition ( ) ( ) and boundary conditions π e sin π e sin π 0 0.5 , 1 0.5 1 The exact solution is given by e sin π , 0 0.5 , e , 0.5 1 The numerical solutions are obtained with 0.02, 0.005 respectively.The closeness of Figure 1 to that of Khuri and Sayfy [12] and the magnitude of absolute errors as a result of comparison between analytic and numerical results in Table 1 demonstrate the reliability of the numerical formulation developed herein.
where ( ) ( ) ( ] , 0,1 0,1 x t ∈ × , , v d are the velocity and diffusion coefficient s respectively.The initial and boundary conditions are given as: Equation ( 21) is a transport equation that involves a convective-diffusive process coupled with a memory term.It is one of the most important equations of mathematical physics and has been used to describe the transport of such quantities like mass, heat, vorticity, energy and momentum.In the absence of the memory term, the analytical solution can be determined by method of separation of variables.This is expressed as: The numerical solution has been found to develop a steep front or boundary layers at x = 1 (Feng and Tian 2006).Profiles of the numerical solution of equation ( 21) and the corresponding analytical solution without the memory term (equation ( 24)) are displayed in Figure 2 21) shows that we have two elliptic operators whose second order derivatives are multiplied by positive parameters.These derivatives model the diffusion process while the only first order derivative is associated with the convective process.Diffusion has always been known to play a significant role in the region of rapid changes (the boundary layers) where the scalar profile exhibits a steep profile.When a standard numerical procedure is applied to a C-D problem, if the diffusion term is relatively smaller than the convective term, the computed solution profile is often oscillatory.Whereas if there exists an "excess" (superfluous) diffusion term, the computed profile is smeared as demonstrated in Figure 2(a).It will also be interesting to compare the effects on the scalar profiles introduced by different magnitudes and discretizations of the convection term.However, the details of this will be the theme of a forthcoming paper.

Conclusion
The main thrust of this paper is a seamless incorporation of non-Fickian flux and memory into a localized element-based boundary element method.The advantage of this technique lies in the possibility of admitting realistic experimental information into a governing differential equation without complicating BEM formulation.For example, it offers a feasible alternative for applying BEM to describe transport processes in viscoelastic media where the inclusion of memory term is vital.
for the source term where

Figure 1 .
Figure 1.3D plot of scalar profiles for Ex. 1.

Table 1 .
Comparison of numerical and analytical results at time = 0.025.