Integral transforms for three-dimensional steady turbulent dispersion in rivers and channels

Applied Mathematical Modelling 31 (2007) 2719–2732 www.elsevier.com/locate/apm

Integral transforms for three-dimensional steady turbulent dispersion in rivers and channels Felipe P.J. de Barros, Renato M. Cotta

*

Center for Analysis and Simulations in Environmental Engineering – CASEE, Mechanical Engineering Department – DEM/POLI and PEM/COPPE, UFRJ Universidade Federal do Rio de Janeiro – Cidade Universita´ria, Cx. Postal 68503, Rio de Janeiro, RJ 21945-970, Brazil Received 1 November 2005; received in revised form 1 August 2006; accepted 24 October 2006 Available online 26 December 2006

Abstract A three-dimensional steady-state mathematical model is considered for predicting the fate of dissolved contaminants in rivers and channels under turbulent flow conditions. The model allows for variable velocity fields and non-uniform turbulent diffusivities within channels of rectangular cross sections. Making use of the Generalized Integral Transform Technique (GITT), a hybrid numerical–analytical solution for the concentration fields is then obtained. The solution convergence behavior is investigated and the criterion for reordering the terms in the infinite series is discussed, with the aim of reducing the computational effort associated with the multiple eigenfunction expansions. A test case is presented to illustrate the proposed approach.  2006 Elsevier Inc. All rights reserved. Keywords: Integral transforms; Contaminants dispersion; Mixing in rivers; Turbulent flow; Hybrid methods

1. Introduction Predicting the fate of pollutants in natural environments such as rivers and man-made channels has become a major concern over the last few years. In order to make accurate and more realistic predictions to aid decision makers, mathematical models must account for the variability of flow properties and chemical discharge in the source terms representation. Due to the difficulties in analytically solving models that incorporate such non-uniform parameters, many authors have preferred to employ numerical methods to solve the partial differential equations that govern the mass transfer processes for turbulent flow in channels and rivers, such as in Nokes and Hughes [1], Lin and Shiono [2], Croucher and O’Sullivan[3] and Mazumder and Dalal [4]. On the other hand, it is also widely accepted that depending on the complexity of the mathematical formulation, these methods tend to involve a considerable computational cost and to be less robust than the so desirable analytical solutions. *

Corresponding author. Fax: +55 21 2562 8383. E-mail addresses: [email protected] (F.P.J. de Barros), [email protected] (R.M. Cotta).

0307-904X/$ - see front matter  2006 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2006.10.024

2720

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732

Despite the difficulties in obtaining exact solutions for turbulent dispersion models, a few previous works attempted to solve this class of problems by making use of analytical techniques. For instance, several authors have used simplified one- and two-dimensional models that to some extent incorporate variable coefficients, Wang et al. [5], Yeh and Tsai [32], Zoppou and Knight [6] and Basha [7]. Approximate solutions were also obtained for two- and three-dimensional steady-state models with variable diffusivities and velocities, Nokes et al. [8] and Nokes and Wood [9] Zoppou and Knight [10] analytically solved the advection–diffusion equation with spatially variable coefficients in three dimensions, but with particular functional forms for the coefficients. Depending on the physical characteristics of the river, such as the aspect ratio and the type of inlet condition of the discharge (example: point, line and plane sources), one- and two-dimensional models fail to properly predict contaminant dispersion within the physical domain near the discharge location. Oneand two-dimensional models are however quite useful when vertical (and lateral) mixing is complete as predicted by Taylor’s dispersion theory, Taylor [11], Taylor [12] and Fischer et al. [13]. For these reasons, three-dimensional contaminant mixing models are of importance in realistically predicting concentration fields near the source discharge, where an environmentally sensitive target may be located and lateral and vertical mixing of the contaminant plume has not yet occurred, Fischer et al. [13]. As for the state of practice in engineering, the knowledge is fairly well consolidated for models that use constant coefficients such as for the velocity and turbulent diffusion coefficients, Lew et al. [14], including the availability of object oriented computer codes, for instance User’s Guide for RIVRISK [15]. Nevertheless, research findings are yet to be pursued in the analytical and robust solution of more generalized models. On the other hand, new advances within the so-called hybrid methods for convection–diffusion problems have been allowing for the treatment of more complex physical situations that could not a priori be handled by the classical analytical methods, such as in situations that involve complex geometries and nonlinearities, maintaining the simplicity, accuracy and relatively small computational cost typical of analytical approaches. Among such hybrid methods, the so-called Generalized Integral Transform Technique (GITT) is here pointed out, as reviewed in Cotta [16], Cotta and Mikhailov [17], Cotta [18], Santos et al. [19], Cotta and Mikhailov [20], as a widely employed approach for the solution of heat and fluid flow problems. The GITT borrows its principles from the Classical Integral Transform Technique, Mikhailov and Ozisik [21], and has inherent advantages due to its hybrid numerical–analytical nature as stated above. This hybrid technique has been largely used as a benchmarking tool and has been recently extended for applications in environmental sciences dealing with dispersion in the three paths, soil, air and water, Cotta et al. [22]. A number of applications dealing with environmental modeling in relation to this hybrid approach can be cited, such as Liu et al. [23], Cotta et al. [24], Almeida et al. [25], de Barros et al. [26]. This hybrid method is also easy to implement computationally, since there is no need for a discretization procedure and mesh adaptive generation, while its major computational task, related to the numerical solution of ordinary differential systems, is readily available as well-tested subroutines in both commercial and public domain subroutines libraries such as the IMSL Library [27] and the Mathematica system, Wolfram [28]. Therefore, the main objective of this paper is to solve a three-dimensional steady-state mathematical model for dispersion of contaminants in rivers and channels under turbulent flow through the use of this hybrid integral transform method. The proposed model includes non-uniform coefficients in any functional form that vary in the longitudinal, vertical and transverse directions. A hybrid numerical–analytical solution is obtained by making use of the GITT, yielding analytical expressions for the space dependence of the concentration fields at any channel cross-section. Since the final solution is expressed as a double summation associated with the eigenfunction expansions in the corresponding space directions (eliminated through the process of integral transforming the original partial differential equation, and analytically recovered by the inverse formula), a series reordering technique is proposed to convert the double summation into a single sum representation. This reordering scheme makes the computational procedure much more efficient once the terms in the eigenfunctions expansions are appropriately reordered according to their predicted individual importance to the final solution, Correa et al. [29], Cotta and Mikhailov [17]. The present paper also demonstrates this hybrid numerical–analytical solution for contaminant transport in rivers and streams as obtained employing mixed symbolic-numerical computations (Mathematica platform). The goal here is to improve and complement existing robust solution implementations to study the dispersion of waste materials in rivers, User’s Guide for RIVRISK [15], based on the integral transform method.

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732

2721

2. Problem formulation A steady-state three-dimensional reactive advection–diffusion equation is considered for mass transport in turbulent flow within rivers and channels. For the sake of illustrating the methodology, it is considered that the river banks, the river bed and surface are not dispersive and these boundaries will be treated with second type boundary conditions (zero flux). Another important aspect to be considered is that the density of the pollutant is approximately equal to the density of the receiving fluid. It is also assumed that the transport in the longitudinal direction is mainly due to advection. The problem formulation can be stated in dimensionless form as     oC oC oC o oC o oC U ðX ; Y ; ZÞ þ fV ðX ; Y ; ZÞ þ W ðX ; Y ; ZÞ ¼f ey ðX ; Y ; ZÞ ez ðX ; Y ; ZÞ þ  kC oX oY oZ oY oY oZ oZ ð1aÞ subjected to the following boundary conditions: oC ¼ 0 for Y ¼ 0; oY oC ¼ 0 for Y ¼ 1; oY oC ¼ 0 for Z ¼ 0; oZ oC ¼ 0 for Z ¼ 1; oZ Cð0; Y ; ZÞ ¼ F ðY ; ZÞ;

ð1bÞ ð1cÞ ð1dÞ ð1eÞ ð1fÞ

where the contamination inlet function F(Y, Z) can assume any general functional form. The dimensionless groups are defined as x X ¼ ; d

z y uðx; y; zÞ vðx; y; zÞ Z ¼ ; Y ¼ ; U ðX ; Y ; ZÞ ¼ ; V ðX ; Y ; ZÞ ¼ ; d x uav uav wðx; y; zÞ C Dyy ðx; y; zÞ ; ; C¼ ; ey ðX ; Y ; ZÞ ¼ W ðX ; Y ; ZÞ ¼ uav uav x C0 Dzz ðx; y; zÞ ked f ðy; zÞ d ; k¼ ez ðX ; Y ; ZÞ ¼ ; F ðY ; ZÞ ¼ ; f¼ ; uav d uav C0 x

ð2Þ

where d is the river’s average depth (L), x is the width (L), uav (L/T) is the average flow velocity, C0 is a reference concentration (M/L3), Dyy(x, y, z) and Dzz(x, y, z) are the dimensional turbulent diffusion coefficients (L2/T) in the transverse and vertical directions, f(y, z) is the inlet condition (at x = 0), ke is the dimensional chemical reaction parameter (T 1), u(x, y, z), v(x, y, z), w(x, y, z) are the variable velocity field components, f

Fig. 1. Geometry and coordinates system for dispersion of contaminants in rivers and channels for various types of inlet conditions (plane, linear and point sources).

2722

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732

is the aspect ratio, k is the dimensionless chemical reaction parameter, and C* is the dimensional concentration (M/L3). Note that x, y and z are the longitudinal, transverse and vertical directions, respectively (see Fig. 1). 3. Solution methodology Following the ideas in the solution methodology through the GITT approach [16], the following integral transform pair is proposed: Z 1Z 1 1 ffiffiffiffiffiffiffiffiffiffiffiffiffi p C in ðX Þ ¼ CðX ; Y ; ZÞwi ðY ÞXn ðZÞ dY dZ; transform; ð3aÞ N yi N zn 0 0 1 X 1 X C in ðX Þwi ðY ÞXn ðZÞ pffiffiffiffiffiffiffiffiffiffiffiffiffi CðX ; Y ; ZÞ ¼ ; inverse ð3bÞ N yi N zn i¼0 n¼0 and the appropriate eigenvalue problems, which shall provide the basis for the expansion of the original potential, C(X, Y, Z), are given below: d2 wi ðY Þ þ b2i wi ðY Þ ¼ 0; dY 2   dwi ðY Þ dwi ðY Þ ¼ ¼0 dY  dY  Y ¼0

ð4aÞ ð4b; cÞ

Y ¼1

d2 Xn ðZÞ þ c2n Xn ðZÞ ¼ 0; dZ 2   dXn ðZÞ dXn ðZÞ ¼ ¼ 0; dZ  dZ  Z¼0

ð5aÞ ð5b; cÞ

Z¼1

where the norms are defined as Z 1 Z 2 N yi ¼ ½wi ðY Þ dY ; N zn ¼ Y ¼0

1 2

½Xn ðZÞ dZ

ð6a; bÞ

Z¼0

and the eigenvalues and eigenfunctions for the vertical and transverse directions are readily obtained in analytic form as cn ¼ np;

Xn ðZÞ ¼ cosðcn ZÞ;

bi ¼ ip;

wi ðY Þ ¼ cosðbi Y Þ

for n; i ¼ 0; 1; 2; . . .

ð7a–dÞ

The next step is to operate the partial differential equation, Eq. (1a), with the following operator, R R natural 1 pffiffiffiffiffiffiffiffiffi ffi 01 01 Xn ðZÞwi ðY Þ dY dZ to yield N yi N zn

1 pffiffiffiffiffiffiffiffiffiffiffiffiffi N yi N zn

Z

1 0

1 þ pffiffiffiffiffiffiffiffiffiffiffiffiffi N yi N zn

Z

1

U ðX ; Y ; ZÞ

0

Z Z

1 0

Z

oC w ðY ÞXn ðZÞ dY dZ oX i

1

fV ðX ; Y ; ZÞ

0

1

Z

oC w ðY ÞXn ðZÞ dY dZ oY i

1

oC w ðY ÞXn ðZÞ dY dZ oZ i 0 0     Z 1 Z 1 1 o oC o oC ey ðX ; Y ; ZÞ ez ðX ; Y ; ZÞ f ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi þ wi ðY ÞXn ðZÞ dY dZ oY oY oZ oZ N yi N zn 0 0 Z 1Z 1 1 wi ðY ÞXn ðZÞkC dY dZ:  pffiffiffiffiffiffiffiffiffiffiffiffiffi N yi N zn 0 0

1 þ pffiffiffiffiffiffiffiffiffiffiffiffiffi N yi N zn

W ðX ; Y ; ZÞ

ð8Þ

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732

2723

Using the orthogonality property and the inverse formula, Eq. (3b), to transform even those terms that are not a priori transformable, we obtain the following system of ordinary differential equations: 1 X 1 X p¼0

1 X 1 dC pj ðX Þ X ¼ /ipnj ðX ÞC pj ðX Þ  kC in ðX Þ; dX p¼0 j¼0 Z 1Z 1 1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi wi ðY ÞXn ðZÞF ðY ; ZÞ dY dZ; N yi N zn 0 0

hipnj ðX Þ

j¼0

C in ð0Þ ¼ F in

ð9aÞ ð9bÞ

where hipnj ðX Þ /ipnj ðX Þ

1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N yi N zn N yp N zj

Z

1 0

Z

1

U ðX ; Y ; ZÞwp ðY ÞXj ðZÞwi ðY ÞXn ðZÞ dY dZ;

  dwp ðY Þ o ey ðX ; Y ; ZÞ wi ðY ÞXn ðZÞfXj ðZÞ dY dZ oY dY 0 0   Z 1Z 1 1 o dXj ðZÞ ez ðX ; Y ; ZÞ wi ðY ÞXn ðZÞwp ðY Þ þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dY dZ oZ dZ N yi N zn N yp N zj 0 0 Z 1Z 1 owp ðY Þ 1 dY dZ  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi fV ðX ; Y ; ZÞXj ðZÞwi ðY ÞXn ðZÞ oY N yi N zn N yp N zj 0 0 Z 1Z 1 1 oXj ðZÞ  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dY dZ: W ðX ; Y ; ZÞwp ðY Þwi ðY ÞXn ðZÞ oZ N yi N zn N yp N zj 0 0

1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N yi N zn N yp N zj

Z

ð10aÞ

0 1

Z

1

ð10bÞ

The coefficients given in Eqs. (10), can be analytically handled through symbolic manipulation packages such as Mathematica, Wolfram [28], and details in the determination of these coefficients can be found in de Barros [30]. Eq. (9) form an infinite system of coupled linear ordinary differential equations with variable coefficients to be solved for the transformed concentration field. Only a truncated version of this system can be solved and the numerical solution is obtained with the aid of a Mathematica subroutine, NDSolve. This Mathematica built in function is quite flexible and robust, allowing the user to specify the desired accuracy and precision, automatically switching to a specialized algorithm when required, whenever the system formed by the ordinary differential equations presents a stiff behavior, Wolfram [28]. However, it is worth noting that the plain truncation of the individual summations present in the double series final solution as represented by Eq. (3b), may not be efficient for computational purposes since this approach can include a number of terms that will not be important in improving the convergence criteria imposed by the user’s prescribed accuracy thus leading to increased computational PP Pcost. One way of avoiding this difficulty is to transform the double summation into a single one, p j ! l , according to an appropriate reordering scheme of the summation terms according to their contribution to the convergence of the final solution. It has been made clear in previous works that used this approach, that the appropriate way of performing multidimensional eigenfunction expansions computation must involve a reorganization of the multi-series represented in Eq. (3b) with an adequate ordering scheme that accounts for, in descending order of magnitude, the most important terms that contribute to the final converged numerical result. As described in [17], the eigenfunction expansion convergence is governed, mainly, by the behavior of the transformed potentials, C in ðX Þ, and these functions are expected to decrease in modulus as the orders i and/or n increase. Since the transformed potentials cannot be known a priori, being part of the numerical solution of the truncated version of the transformed system in Eq. (9) the ordering scheme should be based on well behaved and simple estimates of these quantities, [29]. One of the criteria adopted in the present work is borrowed from a simpler formulation of Eq. (1), with the same boundary conditions but adopting a constant velocity field and constant turbulent diffusion coefficients. For this case, only the longitudinal velocity was adopted as a preferential flow direction, while the others were

2724

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732

neglected. The exact solution for construction of the reordering scheme is then readily given by the method of separation of variables, Cotta [16] Z 1 Z 1  1 X 1 X 1 2 elm;n X Xn ðZÞwm ðY Þ F ðY ; ZÞXn ðZÞwm ðY Þ dY dZ ; ð11aÞ CðX ; Y ; ZÞ ¼ N zn N ym 0 0 m¼0 n¼0 where the global eigenvalue, l, is given by ey0 fb2m þ ez0 c2n þ k ¼ l2m;n ;

ð11bÞ

where ey0 and ez0 are the constant eddy diffusivities values. The convergence behavior of Eq. (11a) is essentially governed by the exponential term. It is of interest to reorganize the double series so as to have first those terms with the smallest arguments in the exponential, since these terms will have the largest contribution to the solution. So, the eigenvalues are reorganized in ascending order following the criteria in Eq. (11b). Notice, however, that we have used a much simpler formulation to obtain the reordering scheme, and this might not be the optimal solution for the highly variable coefficients case. Another straightforward procedure would be to follow the behavior of an approximate solution of the original problem, for instance, by holding only the main diagonal terms of the coefficients matrices in the transformed problem. To illustrate this socalled lowest order solution [16], we rewrite below the system of ordinary differential equations for the approximate transformed potentials, similar to Eq. (9a): dC lin ðX Þ ¼ /iinn ðX ÞC lin ðX Þ  kC lin ðX Þ; dX C lin ð0Þ ¼ F in :

hiinn ðX Þ

ð12aÞ ð12bÞ

Here the superscript l in the transformed potential is used to denote the lowest order solution. The system is now decoupled, and the transformed potentials are explicitly obtained, offering a direct formula for the reordering procedure  Z X   /iinn ðX 0 Þ  k l 0 dX : C in ðX Þ ¼ F in Exp  ð13Þ hiinn ðX 0 Þ 0 This lowest order solution is an approximation to the unknown solution of the original problem formulation, and allows for the grouping of the indices i and n so as to account for the most significant terms to the final unified sum in the appropriate sequence. The level of accuracy achieved by this solution is not a major concern in this context, but it is expected to present similar convergence trends as the unknown solution, and in this case it should allow even for a dynamic reordering procedure, modifiable along the longitudinal coordinate of the numerical integration procedure. In addition, the analytically iterated version of this solution, when the non-diagonal terms are approximately accounted for as a source term [16, 20], may be used for a more refined reordering, when required, for a cost effective convergence enhancement and system size reduction. Further discussion on reordering schemes may also be found in [17,20]. 4. Results and discussion First of all, in order to verify that the computer code has been appropriately developed, it is important to inspect if it can reproduce the solution of a simpler model, e.g., with constant coefficients. The numerical verification of the proposed model can be seen in Fig. 2. The input data are as follows: eyav = 0.002, ezav = 0.02, Uav = 1, Vav = 0, Wav = 0, the aspect ratio (f) is equal to 0.1, and the inlet condition (C0 = 1) is bounded by Y1 = 0.4, Y2 = 0.6, Z1 = 0.4, and Z2 = 0.6. The subscripts av indicate average values. The exact solution was obtained through the method of separation of variables, as in Eq. (11a), and the formulation is identical to Eq. (1a) but with constant coefficients (average values). As one may see, the two curves perfectly match, showing that the developed Mathematica notebook is able to recover the solution of a simpler formulation with known exact solution. We now proceed to analyze the behavior of the three-dimensional solution in comparison with a two-dimensional horizontal steady-state model [26]. Figs. 3 and 4 illustrate the dimensional concentration

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732

2725

1.0

0.8

Exact Solution GITT

0.6

Y0 = Z0= 0.5

C 0.4

0.2

0.0 0

50

100

150

X Fig. 2. Numerical verification of the hybrid numerical–analytical solution (GITT) with exact solution for Y0 = 0.5 and Z0 = 0.5.

0.0030

2D 3D: ε zav= 0.00001

0.0025

3D: ε zav= 0.00005

C*[mg/l]

0.0020

3D: ε zav= 0.0001 3D: ε zav= 0.0003

0.0015

0.0010

0.0005

0.0000 0

500

1000

X [m] Fig. 3. Concentration (mg/l) curves along the longitudinal direction (m) at y = 0 m for the two-dimensional model [26] and for the threedimensional model with an average vertical diffusivity (ezav = 105; 5 · 105; 104 and 3 · 104 m2/s).

distribution throughout the river/channel while progressively reducing the values for the vertical turbulent diffusion coefficient (Fig. 3) and the aspect ratio (Fig. 4). In Fig. 3, the aspect ratio was maintained constant and equal to 0.1. For the case in which the aspect ratio was decreased, the turbulent diffusion coefficient was kept uniform and equal to 3 · 104 m2/s. From Figs. 3 and 4, we can observe that the three-dimensional solution tends to the horizontal two-dimensional solution with both the decrease of the mixing coefficient in the vertical direction and with the reduction of the aspect ratio. As expected, the three-dimensional curves are located below the two-dimensional one, due to the vertical turbulent diffusion still present in the 3D cases, though increasingly less significant. In Fig. 4, as the aspect ratio is decreased, for instance, the channel depth is decreased while the width is maintained constant and the channel/river becomes shallower. The hypothesis that the contaminant is immediately dispersed towards the river bed and surface might then become reasonable, depending on the aspect ratio considered. This hypothesis is the basis of the two-dimensional solution presented in Fig. 4, more closely studied in [31]. It can be observed that as the channel aspect ratio is decreased, a longer longitudinal distance is required for the accumulated effects of vertical dispersion to play some role. Another interesting observation is that the shallower the river, the smaller will be the characteristic

2726

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732 0.0035 0.0030

2D 3D: ζ = 0.01 3D: ζ = 0.05 3D: ζ = 0.075 3D: ζ = 0.1 3D: ζ = 0.3 3D: ζ = 0.5

C* [mg/l]

0.0025 0.0020 0.0015 0.0010 0.0005 0.0000 0

500

1000

X [m] Fig. 4. Comparison of the three-dimensional and two-dimensional solutions with the reduction of the aspect ratio of the channel/river (ezav = 3 · 104 m2/s).

spatial scale in the vertical direction, thus limiting the size of the turbulent eddies and consequently diminishing the effect of vertical dispersion. It is important to state that the present analysis is only aimed at the numerical investigation of the proposed solution, since there is an actual relation between these two physical parameters (vertical eddy diffusion and aspect ratio). Based on a RIVRISK example found in User’s Guide for RIVRISK [15], together with velocity and diffusivity profiles taken from Yeh and Tsai [32], the test case presented in this section illustrates how the GITT is used to predict concentration fields within three-dimensional steady situations. Therefore, consider a river with the following longitudinal velocity and diffusivity profiles: uðzÞ ¼ a1 zq1 ; q2

Dzz ðzÞ ¼ a2 z ;

ð14aÞ ð14bÞ

where a1, a2, q1 and q2 are constant parameters that depend on the friction at the river bed and on the Reynolds number. The expressions for a1 and a2 are given as [32] uav dðq1 þ 1Þ ; d q1 þ1 Dzzav dðq2 þ 1Þ a2 ¼ ; d q2 þ1

a1 ¼

ð15aÞ ð15bÞ

Dzzav and Dyyav are taken as constant average values for the turbulent diffusion coefficients. For high Reynolds number one may take q1 = 1/7 and q2 = 6/7, Yeh and Tsai [32]. In the present example, the transversal turbulent diffusion coefficient is thus considered constant, following User’s Guide for RIVRISK [15]. It is also considered that the transversal and vertical velocity components are negligible. Figs. 5 and 6 illustrate the dimensionless velocity and vertical diffusivity profiles employed, as given by Eqs. (14a), (14b), (15a), (15b). The average velocity of the river is 0.3 m/s and the discharge is continuous and given in the form of a plane source limited by the dimensionless coordinates Y1 and Y2 in the horizontal plane and Z1 and Z2 in the vertical plane. The input data is consolidated in Table 1. The convergence behavior of the concentration field for two distinct aspect ratio (f = 0.05 and f = 0.2125) are illustrated in Tables 2a and 2b for a fixed position Y = 0.075 and Z = 0.65 in the Y–Z plane, but varying the longitudinal coordinate, with X = 1, 3, 5, and 8. The dimensionless concentration values are presented in terms of the reordered single expansion truncation order, M. The two sets of results serve to illustrate the robustness of the approach while dealing with two fairly different geometric configurations, and confirm that there is not a remarkable effect of the aspect ration on the expansions convergence behavior. The convergence

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732

2727

Fig. 5. Dimensionless longitudinal velocity profile for test case.

Fig. 6. Dimensionless vertical diffusivity profile for test case.

Table 1 Input data for simulations with f = 0.05 and f = 0.2125 w (m) d (m) f Z1 Z2 Y1 Y2 Dzzav (m2/s) Dyyav (m2/s) C0 (mg/l)

80 4 and 17 0.05 and 0.2125 0.58 0.705 0 0.125 0.017 0.05 1.0

rate, as in general expected for eigenfunction expansions, is slower for regions closer to the inlet, X = 1, when compared to the regions further away. With a truncation order M = 420, the dimensionless concentration value at X = 1 is converged to three significant digits while at X = 8, the concentration result is already converged to the fourth significant digit. On the other hand, it can be observed that the concentration is accurate to within 1% all over the inspected region for fairly low truncation orders, and thus very low computational cost. The dimensionless concentration distribution for the fixed coordinates, Y0 = 0.075 and Z0 = 0.65 and Y0 = 0 and Z0 = 1, along the longitudinal distance (dimensionless), are illustrated in Figs. 7a, 7b and 8. For a source in a centered position within the channel (Figs. 7a and 7b), one may observe the concentration

2728

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732

Table 2a Convergence behavior of the dimensionless concentration distribution along longitudinal coordinate for aspect ratio f = 0.2125 C(X, Y, Z)

Y = 0.075 and Z = 0.65

M

X=1

X=3

X=5

X=8

20 60 100 140 180 220 260 300 340 380 420

0.3094 0.4649 0.5370 0.5863 0.5885 0.5943 0.5953 0.5958 0.5956 0.5953 0.5952

0.2500 0.3112 0.3224 0.3244 0.3244 0.3244 0.3244 0.3244 0.3244 0.3244 0.3244

0.2079 0.2347 0.2363 0.2362 0.2362 0.2362 0.2362 0.2362 0.2362 0.2362 0.2362

0.1655 0.1738 0.1738 0.1738 0.1738 0.1738 0.1738 0.1738 0.1738 0.1738 0.1738

Table 2b Convergence behavior of the dimensionless concentration distribution along longitudinal coordinate for aspect ratio f = 0.05 C(X, Y, Z)

Y = 0.075 and Z = 0.65

M

X=1

X=3

X=5

X=8

20 60 100 140 180 220 260 300 340 380 420

0.2047 0.4356 0.5228 0.6011 0.6108 0.6252 0.6390 0.6377 0.6355 0.6393 0.6392

0.1908 0.3400 0.3730 0.3881 0.3888 0.3896 0.3891 0.3890 0.3890 0.3889 0.3889

0.1796 0.2762 0.2892 0.2917 0.2916 0.2912 0.2913 0.2913 0.2913 0.2913 0.2913

0.1662 0.2169 0.2197 0.2198 0.2198 0.2196 0.2197 0.2197 0.2197 0.2197 0.2197

ζ = 0.2125 1.0

Y0=0.075 and Z0=0.65

C(X, Y0, Z0)

0.8

0.6

0.4

0.2

0.0 0

2

4

6

8

X

Fig. 7a. Concentration distribution along the dimensionless longitudinal distance for f = 0.2125 with fixed dimensionless coordinates, Y0 = 0.075 and Z0 = 0.65.

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732

2729

ζ = 0.05 1.0

Y0=0.075 and Z0=0.65

C(X, Y0, Z0)

0.8

0.6

0.4

0.2

0.0

0

5

10

15

20

25

30

35

X

Fig. 7b. Concentration distribution along the dimensionless longitudinal distance for f = 0.05 with fixed dimensionless coordinates, Y0 = 0.075 and Z0 = 0.65.

0.10

Y0=0 and Z0=1

C(X, Y0, Z0)

0.08

0.06

0.04

0.02

0.00 0

2

4

6

8

X

Fig. 8. Concentration distribution along the dimensionless longitudinal distance with fixed dimensionless coordinates, Y0 = 0 and Z0 = 1 (at the river surface and left bank) for f = 0.2125.

value decay along the river length, while for the river surface and left bank (Fig. 8), the concentration build up is clearly observed as a result of the source contamination dispersion. Fig. 9 shows the different concentration profiles for various dimensionless longitudinal distances: X = 0.5; 1; 2; 3; 4 and 5 for Z0 = 0.65 (middle of the plane source). This figure allows to visualize the transversal dispersion of the contaminant plume for different values of X and a fixed vertical position Z0 = 0.65. The contour plots are in Figs. 10 and 11 to complement such analysis, by illustrating the contaminant plume dispersion in both the transversal and vertical planes. It can be noticed in Fig. 12 for the concentration profiles in the Y–Z plane, that the plume spreads more widely, and in an asymmetric fashion, towards the surface of the river. This behavior is due to the larger values of the velocity and of the turbulent diffusivity near the surface, as indicated by Eqs. (14a) and (14b) and Figs. 5 and 6. The three-dimensional nature of the proposed test case is clearly observable from these plots, and its treatment via a two-dimensional modeling, in either the X–Y or the X–Z planes choices, would certainly lead to substantial numerical differences and even physical deviations.

2730

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732 1.0

C(X, Y, Z0)

0.8

X = 0.5 X=1 X=2 X=3 X=4 X=5

0.6

0.4

Z 0 = 0.65 0.2

0.0 0.0

0.2

0.4

0.6

0.8

1.0

Y Fig. 9. Concentration distribution (dimensionless) versus transversal distance for X = 0.5; 1; 2; 3; 4 and 5 for f = 0.2125. 0.5 0.4

Y

0.3 0.2 0.1 0 0

0.5

1

1.5

2

2.5

3

X Fig. 10. Concentration isolines (0.75; 0.5; 0.25; 0.1 and 0.025) in the horizontal plane for Z0 = 0.65 using f = 0.2125.

Z

1 0.9 0.8 0.7 0.6 0.5 0.4 0

0.5

1

1.5

2

2.5

3

X Fig. 11. Concentration isolines (0.75; 0.5; 0.25; 0.1 and 0.025) in the vertical plane for Y0 = 0.075 using f = 0.2125.

0.6 1

C(1,Y, Z) 0.4 0.2

0.75

0 0

0.5 Z

0.25 0.5 Y

0.25 0.75 1

0

0.2 0.15 C(8,Y, Z) 0.1 0.05 0 0

1 0.75 0.5 Z

0.25 0.5 Y

0.25 0.75 1

0

Fig. 12. Three-dimensional plots of the dimensionless concentration field in the Y–Z plane for X = 1 and X = 8.

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732

2731

5. Conclusions The Generalized Integral Transform Technique (GITT) was advanced to obtain numerical–analytical solutions for contaminants dispersion in rivers and channels. The proposed model is a three-dimensional steadystate formulation with variable coefficients and within channels of rectangular cross section. As reviewed in the literature, three-dimensional models are very useful to characterize mass transport near source location where complete vertical or lateral mixing of the contaminant plume has not occurred. For these conditions, three-dimensional models yield more accurate results for early mixing times and consequently will lead to a more refined risk assessment. Although data acquisition for velocity and eddy fields are difficult to obtain for environmental flows, the use of approximate and estimated profiles found in the literature can provide important and less conservative results, thus providing a more realistic view of the actual physical process. By making use of an eigenfunction expansions reordering scheme, the computational cost of the GITT solution was reduced and numerical results of the concentration profiles, converged to the fourth significant digit, were presented in tabular form for reference purposes. This is clearly an attractive feature of the integral transform method, allowing us to obtain semi-analytical solutions for complicated models with a low computational cost. It is worth mentioning that the incorporation of time dependence and of the longitudinal turbulent diffusion term in the model, are not limiting factors for the integral transform method here proposed, and can be incorporated in future implementations of this methodology, following the generalization previously pursued in heat transfer and fluid flow problems [16,17,33,20]. As a final remark, we may conclude that the GITT can be used for benchmarking purposes in combination with other numerical methods, but also serve as an alternative solution method with robust analytical characteristics, for contaminant dispersion problems. Acknowledgements The authors would like to acknowledge the partial financial support provided by CNPq/Brazil, Tetra Tech and EPRI/USA, sponsors of CASEE, and other Brazilian agencies and programs, FAPERJ and PRONEX (‘‘Nu´cleo de Exceleˆncia em Turbuleˆncia’’). The authors also express their gratitude to Dr. William B. Mills from Tetra Tech, Inc., Lafayette, USA, author and coordinator of the RIVRISK software system. References [1] R.I. Nokes, G.O. Hughes, Turbulent mixing in uniform channels of irregular cross-section, J. Hydraul. Res. 32 (1) (1994) 67–86. [2] B. Lin, K. Shiono, Numerical modeling of solute transport in compound channel flows, J. Hydraul. Res. 33 (1) (1995) 773–788. [3] A.E. Croucher, M.J. O’Sullivan, Numerical methods for contaminant transport in rivers and estuaries, Comput. Fluids 27 (8) (1998) 861–878. [4] B.S. Mazumder, D.C. Dalal, Contaminant dispersion from an elevated time dependent source, J. Comput. Appl. Math. 126 (2000) 195–205. [5] S.T. Wang, A.F. McMillan, B.H. Chen, Dispersion of pollutants in channels with non-uniform velocity distribution, Water Res. 12 (1978) 389–394. [6] C. Zoppou, J.H. Knight, Analytical Solutions for Advection and Advection–Diffusion Equations with Spatially Variable Coefficients, J. Hydraul. Eng. February (1997) 144–148. [7] H.A. Basha, Analytical model of two dimensional dispersion in laterally nonuniform axial velocity, J. Hydraul. Eng. 123 (10) (1997) 853–862. [8] R.I. Nokes, A.J. McNulty, I.R. Wood, Turbulent dispersion from a steady two-dimensional horizontal source, J. Fluid Mech. 149 (1984) 147–159. [9] R.I. Nokes, I.R. Wood, Vertical and lateral turbulent dispersion: some experimental results, J. Fluid Mech. 187 (1988) 373–394. [10] C. Zoppou, J.H. Knight, Analytical solutions of a spatially variable coefficient advection–diffusion equation in up to three dimensions, Appl. Math. Model. 23 (1999) 667–685. [11] G.I. Taylor, Dispersion of soluble matter in solvent flowing slowly through a tube, Proc. Roy. Soc. A 219 (1953) 186–203. [12] G.I. Taylor, Dispersion of matter in turbulent flow through a pipe, Proc. Roy. Soc. A 223 (1954) 446–468. [13] H.B. Fischer, J.E. List, R.C.Y. Koh, J. Imberger, N.H. Brooks, Mixing in Inland and Coastal Waters, Academic Press, 1979. [14] C.S. Lew, W.B. Mills, J.Y. Loh, Power plant discharges of total residual chlorine and trihalomethanes into rivers: potential for human health and ecological risks, Hybrid Methods Eng. 1 (1999) 19–36. [15] User’s Guide for RIVRISK, 2000. Version 5.0: A Model to Assess Potential Human Health and Ecological Risks from Power Plant and Industrial facility Releases to Rivers, EPRI, Palo Alto, CA, 2000.

2732

F.P.J. de Barros, R.M. Cotta / Applied Mathematical Modelling 31 (2007) 2719–2732

[16] R.M. Cotta, Integral Transforms in Computational Heat and Fluid Flow, CRC Press, 1993. [17] R.M. Cotta, M.D. Mikhailov, Heat Conduction: Lumped Analysis, Integral Transforms, Symbolic Computation, John Wiley & Sons, 1997. [18] R.M. Cotta (Ed.), The Integral Transform Method in Thermal and Fluids Science and Engineering, Begell House, New York, 1998. [19] C.A.C. Santos, J.N.N. Quaresma, J.A. Lima (Eds.), Convective Heat Transfer in Ducts: The Integral Transform Approach, Mechanical Sciences Series, Brazilian Society of Mechanical Sciences (ABCM), E-Papers, Rio de Janeiro, 2001. [20] R.M. Cotta, M.D. Mikhailov, Hybrid Methods and Symbolic Computations, in: W.J. Minkowycz, E.M. Sparrow, J.Y. Murthy (Eds.), Handbook of Numerical Heat Transfer, second ed., John Wiley, New York, 2006, pp. 493–522 (Chapter 16). [21] M.D. Mikhailov, M.N. Ozisik, Unified Analysis and Solutions of Heat and Mass Diffusion, John Wiley, NY, 1984, also, Dover Publications, 1994.. [22] R.M. Cotta, M.J. Ungs, M.D. Mikhailov, H.R.B. Orlande, P.F.L. Heilbron, Hybrid methods and mixed computations for contaminant transport in environmental modeling, in: 3rd International Conference on Computational Heat and Mass Transfer, CHMT-2003, Invited Keynote Lecture, Banff, Canada, May, 2003. [23] C. Liu, J.E. Szecsody, J.M. Zachara, W.P. Ball, Use of the generalized integral transform method for solving equations of solute transport on porous media, Adv. Water Resour. 23 (2000) 483–492. [24] R.M. Cotta, M.D. Mikhailov, M.J. Ungs, Contaminant transport in finite fractured porous medium: integral transforms and lumpeddifferential formulations, Ann. Nucl. Energy 30 (3) (2003) 261–285. [25] G.L. Almeida, L.C.G. Pimentel, R.M. Cotta, Analysis of pollutants atmospheric dispersion from a point source using integral transforms, in: 3rd International Conference on Computational Heat and Mass Transfer, CHMT-2003, Banff, Canada, May, 2003. [26] F.P.J. de Barros, R.M. Cotta, Hybrid solutions for contaminants dispersion in rivers: application to biocides contamination from hydroelectric power plants, in: 6th Computational Modeling Meeting, IPRJ/UERJ, Friburgo, RJ, Brazil (CD-ROM), December, 2003. [27] IMSL Library, 1991, MATH/LIB, Houston, TX. [28] S. Wolfram, The Mathematica Book, fourth ed., Wolfram Media, Cambridge, 1999. [29] E.J. Correˆa, R.M. Cotta, H.R.B. Orlande, On the reduction of computational costs in eigenfunction expansions of multidimensional diffusion problems, Int. J. Numer. Methods Heat Fluid Flow 7 (7) (1997) 675–695. [30] F.P.J. de Barros, Multidimensional models for contaminants dispersion in rivers and channels: hybrid solutions via integral transforms, MSc Thesis, COPPE/UFRJ, Rio de Janeiro, Brazil, 2004 (in Portuguese). [31] F.P.J. de Barros, W.B. Mills, R.M. Cotta, Integral transform solution of two-dimensional model for contaminants dispersion in rivers and channels with spatially variable coefficients, Environ. Model. Software 21 (2006) 699–709. [32] G.T. Yeh, Y.J. Tsai, Dispersion of water pollutants in a turbulent shear flow, Water Resour. Res. 12 (6) (1976). [33] H. Luz Neto, J.N.N. Quaresma, R.M. Cotta, Natural convection in three-dimensional porous cavities: integral transform method, Int. J. Heat Mass Transfer 45 (14) (2002) 3013–3032.