Pressure correction DRBEM solution for heat transfer and fluid flow in incompressible viscous fluids

EngineeringAnalysi,~withBoundaryElementslg (1996)265-272

Pll: S0955- 7997(97)00003-9

ELSEVIER

© 1997El~vierScienceLtd All rightsreserved.Printedin GreatBritain 0955-7997/96/$15.00

Pressure correction DRBEM solution for heat transfer and fluid flow in incompressible viscous fluids Christopher P. Rahaim & Main J. K a s s a b Department of Mechanical, Materials, and Aerospace Engineering, Institute for Computational Engineering, University of Central Florida, Orlando, Florida 32816-2450, USA

(Received 15 August 1996; revised version received 8 January 1997; accepted 10 January 1997)

The Dual Reciprocity Boundary Element Method (DRBEM) is used to solve incompressible laminar viscous fluid flows and heat transfer. The DRBEM is extended to develop a pressure correction scheme to solve the incompressible Navier-Stokes equations. The velocity field is then used as input to the DRBEM solution of the energy transport equation, thereby retaining the boundary only discretization feature of the BEM for the solution of this problem. Numerical results for the proposed DRBEM solution for laminar flow and heat transfer in a channel are obtained for several Reynolds numbers and compare well with previously published data. © 1997 Elsevier Science Ltd Key words: DRBEM, viscous flow and heat transfer, pressure correction,

boundary element method. 1 INTRODUCTION

discretizes both the domain boundary as well as the domain interior. The latter poses a large computational burden for even the most simple problem. There is another class of numerical methods which efficiently solve field problems using only a boundary discretization, and these methods are consequently referred to as boundary methods. The interest in these techniques is that they reduce the computational effort involved in obtaining a solution to a field problem by reducing the problem dimension by one. The boundary element method (BEM) is a numerical method for the solution of field problems which finds its basis in the theory of integral equations. The BEM is classified as direct, if the unknowns appearing in the integral equations are related to physical variables, or indirect, if the unknowns appearing in the integral equations are not related to physical variables. The panel method is an example of indirect BEM which has found much success in the aerospace industry. Although the BEM has many advantages, one of its important features is the boundary only discretization required to solve finear problems. The other important feature of the BEM stems from its basis in integral equation formulation which leads to numerical integration (a stabilizing numerical procedure) rather than differentiation.

The governing equations for nonisothermal viscous fluid flow are the Navier-Stokes (NS) equations supplemented by the energy transport equation. Collectively, these are coupled nonlinear partial differential equations which are not amenable to analytical solutions except in a few idealized limiting cases. However, the solution of these equations for general configurations is of obvious technical interest as it provides important engineering quantities such as drag coefficients, pressure drops, and heat transfer rates. Consequently, there has been much work devoted to developing efficient and broadly applicable numerical algorithms to solve the NS equations, and several numerical techniques have shown success in this effort. Currently, the most popular methods are based on the finite element method (FEM), the finite difference method (FDM), and, more recently, the finite volume methods (FVM). However, these methods can be thought of as belonging to a broad class of whole domain numerical methods. This is due to the fact that they require the solution of a comparatively large set of sparse algebraic equations resulting from the discretization of the original NS equations on a computational grid (or mesh) which 265

266

C. P. Rahaim, A. J. Kassab

In this paper, we discuss a direct BEM for the solution of the NS equations and the energy equation. The use of Green's free space (or fundamental) solutions to arrive at analytical solutions of heat transfer problems dates back to Morse and Feshbach 1 and Carslaw and Jaeger. 2 Shaw 3 and Chang et al. 4 used the fundamental solution of the diffusion equation to develop a method now known as the BEM for solving transient heat transfer problems. Wrobel and Brebbia 5 subsequently extended the approach and developed the formalism of the BEM for heat transfer problems. The method has since reached a stage of maturity, and nonlinear heat conduction problems are now routinely solved using the BEM. A recent monograph edited by Wrobel and Brebbia 6 discusses some of the recent advances in BEM in heat transfer. A derived variable velocity-vorticity BEM formulation was pioneered by Wu 7 and his co-workers to solve the NS equations for viscous fluid flow. Brebbia and Skerget 8 developed a derived variable BEM for incompressible viscous flow and heat transfer problems in channels. Skerget et al.9 applied the BEM to bouyancy driven viscous cavity flows. These applications of the BEM to the NS lead to domain integrals arising from treatment of convective terms as pseudo-source terms. Although this detracts from the boundary only discretization feature of the BEM, results show that very accurate solutions can be obtained with a coarse discretization in comparison to that required by FDM, FVM and FEM. This is due to the incorporation of the Green's free space solution in the numerical algorithm. Recently Wrobel and DeFigueiredo 1° and Qiu et al. 11 proposed a dual reciprocity boundary element method (DRBEM) for the solution of the convection-diffusion equation in an attempt to eliminate domain integration. Results of these studies suggest that the DRBEM is a promising technique for the solution of viscous flow problems. The recent monographs edited by Power et al. 12 and by Power 13 report several applications of the BEM in viscous fluid flows including modeling of turbulent flows using primitive variables, a one equation turbulence model, and a relatively coarse domain discretization for the nonlinear convective terms. In this paper, we extend our previous work on the DRBEM solution of a modified nonlinear Burgers equation, 14 to develop a DRBEM primitive variable based pressure correction scheme to solve the NS equations and a DRBEM to solve the energy equation. DRBEM computed velocity profiles, Nusselt numbers and surface wall temperatures for channel flow are obtained for several Reynolds numbers. Comparison with previously reported data reveal that the DRBEM accurately predicts the investigated flow fields. 2 GOVERNING EQUATIONS The governing equations for incompressible nonisothermal

Newtonian fluid flow are the continuity equation V. 1,7= 0

(1)

the conservation of linear momentum p(lT. VIT) = - V p + #V2I7

(2)

and the energy transport equation

pCp( V. VT) = kV2T + #¢~

(3)

Here, p is the density, # is the dynamic viscosity, k is the thermal conductivity, Cp is the specific heat at constant pressure, p the pressure, V is the velocity vector, T is the temperature and • is the dissipation function. Equations (1) and (2) constitute the NS equations. FVM and FDM solutions of these equations are commonly obtained by an iterative pressure correction scheme. In this method, the discrete form of eqn (2) is solved with an assumed or current estimate of the pressure field to yield a velocity field. The velocity field will generally not satisfy the conservation of mass equation. Therefore, a Poisson equation is solved to update the pressure field. Subsequently, once the velocity field is resolved, it is used to solve the energy transport equation. Below we develop a pressure correction DRBEM solution of the NS equations and a DRBEM solution of the energy transport equation.

3 BEM SOLUTION OF THE NS AND ENERGY TRANSPORT EQUATIONS In an effort to develop a DRBEM pressure correction scheme, we recast the two dimensional NS equations as

v2 =! #

@+p( ou

Oy)J

(4)

it -~y + P UOx + Oy,]J

V2 v = 1 0 p

and adjoin the pressure correction equation

( Ou Ov VZp=2p -~x Oy

Ou Ov) Oy Ox

(5)

and the energy transport equation

Vz T

pCp I" OT

- --c

+~

+

2 ~-~

OT'~ )

(ov+ou 1

+2 ~yy + Oxx OyJ j

(6)

A close examination of eqns (4)-(6) reveals that they have the general form

(x Ou Ou Ov Ov Op Op OT 0_~) V 2 u = b , , 'Y'U'v'p'T'ox'Oy'Ox'Oy'Ox'Oy' Ox' v . r ]

(7)

Pressure correction D R B E M solution Here, all nonlinear terms have been lumped into the right hand side forcing function b following Refs 1-3. Thus, we first focus on eqn (7) for the development of the basic steps in the BEM. Introducing the fundamental (or Green's free space) solution to the Laplace equation, u*, and its normal derivative, q*, eqn (7) is recast into the following boundary integral equation using Green's second identity C(,)u(,) + Jr uq*dF = Iv qu*dr + In bu*dfl

(8)

where F is the domain boundary of the domain, f~, and C(~) is one when u is evaluated for an interior point, E f L and a fraction of one when u is evaluated for a boundary point, ( E F. 15 It is noted that the above equation provides the value of u at any interior point, ~, of the domain, fL if both u and its normal derivative, q, are known on the bounding curve, F, and if the forcing term, b, is known in fL As this is not the case in a well posed problem, a boundary value problem is formed to solve for the boundary unknowns by evaluating the above equation at the boundary. Using standard BEM formulation, a pattern of nodes is laid out on the boundary, F, and the boundary is discretized into boundary elements. The distributions of u and q on the boundary are then interpolated using their values at boundary nodes. In our code, we use a linear interpolation for the boundary geometry as well as for u and q. Consequently, eqn (8) is discretized in matrix

267

eqn (7) is expanded as ND

b(x, y) = E aJfJ (x, y)

(10)

j=i

The expansion functions are chosen to satisfy the following useful property:

fJ(x,y) = v2~J(x,y)

(11)

Using the above property and Green's second identity, the domain integral in eqn (8) is transformed into a series of contour integrals as

J.

ND

+ [ OfJu*dr.IrOn

fJq*dV

(12)

There are several choices for the expansion functions fY which satisfy the crucial relation that f J is a particular solution ofeqn (11). Fourier series as well as radial polynomial expansion functions have been successfully used) 5 A recent review of the state of the art in the treatment of the expansion functions is given in Goldberg) 6 We choose to use a radial interpolating function in this paper. In particular, we employ the following expansion function:

fJ(x,y) = 1 + r(x,y)

(13)

fornl as

[n]{u} = [Gl{q} + {B}

(9)

where {q} is the vector nodal normal derivatives of u, {u} is the vector of nodal values of u and vector {B} is obtained by evaluating the domain integral in eqn (8). The influence matrices, [HI and [G], are functions of the geometry, and they are readily evaluated by use of quadratures. IntrodUcing boundary conditions, eqn (9) is solved for the boundary unknowns. This is the standard approach taken in previous BEM studies of the NS.

4 DRBEM TREATMENT OF THE NS EQUATIONS

AND ENERGY TRANSPORT EQUATION In the DRBEM, a special series expansion of the forcing term b is employed to turn the domain integral in eqn (8) into a series of contour integrals, thereby avoiding domain discretization. In what follows, we closely follow the developments described in Partridge et al., 15 and the reader is directed to that reference for details of the DRBEM. As a first step, we let the forcing term, b, in eqn (7), be only a function of geometry, e.g. b = b(x, y) only. Using the global approximating functions f / ( x , y ) and ND dual reciprocity expansion points, 15 the forcing term in

where r(x,y) = ~/(xi - x) 2 + (Yi - y)2 is the final radial distance measured from the dual reciprocity expansion point (xi,yi) to any point of interest (x,y). It should be noted that we have tried the 1 + r + r 2 and 1 + r + r 2 + r 3 expansion functions. The 1 + r expansion function provided the best results in the problems we considered in this paper. The results obtained with the 1 + r expansion function were satisfactory, as will be shown in the example section, and thus did not warrant use of other expansion functions, for example, thin plate splines. However, there may very well be cases other than those investigated in this study which may warrant trial of other expansion functions. However, in this paper, we laid out the procedure to carry out a full pressure-correction solution of the viscous NavierStokes equations using DRBEM, and it is just a matter of implementation to change the expansion function in eqn (13). There is no current standard method to distribute the dual reciprocity expansion points. The dual reciprocity points may be chosen to simply coincide with the boundary element nodes. However, experience has shown that, in addition to boundary points, there must be some internal dual reciprocity points for accurate resolution of highly nonlinear problems such as the one we are addressing) 5 A typical distribution of dual reciprocity points and boundary element nodes

C. P. Rahaim, A. J. Kassab

268

for a channel problem is shown in Fig. 1 (later). Requiring that the series expansion in eqn (10) collocates the function b at the ND dual reciprocity points leads to an equation for the expansion coefficients ay as {a} = [F]-I{b}

(14)

where the vector {b} contains the values of the forcing term b at the ND dual reciprocity points. The square matrix IF] is formed by vectors of the interpolating functions ff(x,y) evaluated at the dual reciprocity points. Introducing the above result and the BEM diseretization into eqn (12), results in the matrix equation {a} = [S]{b}

(15)

0p

br'x(X'Y) = "~x

(21)

bp,y(X,y) = Op

(22)

oy

lead to the following DRBEM vectors

{ap,x(P)} :

(23)

IS] 00~ [F]-l{p}

{B,,y(p)} = [S] ~

(24)

[r]-I {p}

t/y

where the vector {p} contains the dual reciprocity nodal values of the pressure. 3. The pressure correction equation forcing term

(On Ov OU~x) b.c(x,y) = -~x Oy Oy

where the matrix [S] is given as [S] = ([HJ[U]- [GJ[Q])[F] -1

(25)

(16) leads to the following DRBEM vector

The matrices [H] and [G] are the same as those in eqn (9), while the matrices [U] and [Q] are geometry dependent matrices containing the function ~ and its normal derivative ~ (Ref. 7). It is noted that the matrix IS] is purely a function of geometry, and, once computed, this matrix does not change, regardless of the makeup of the forcing function b. Having established the basic approach to address an arbitrary forcing function of position, the basic form in eqn (15) is then used to address the right hand side forcing functions in eqns (4)-(6). Generalized expressions can be derived for various forcing terms appearing in the NS using the above outline approach. In particular, the following terms are converted using the DRBEM: I. The convective terms in the x- and y-momentum equations

( Ou yOU) b~,x(X,Y)= U-~x+ Oy

(17)

( ov

bc,y(x,y)= u~x+v

(18)

lead to the following DRBEM vectors 0[F]

(spc(u,v)}

o[F]

=

[S~ (__0_~_x [el\

O[Fl [F]-l[udl-~x]) [F]-' (26) 4. The convective term in the energy transport equation

OT + vOT'~ bc'r(U'v' T) = U-~x Oy J

(27)

leads to the following DRBEM vector

O[F] , IVd] -~-y} O[F]~

{ac,T(u, ~, r)} = Is] [v d] ~ . x [FI-t{T}

(28)

5. The dissipation function in the energy transport equation

br,~(u,v)=

20xx

[ Ov

0[F]'~

-1 d O[F] [vl N

+2 ~yy

2]

Ou'~

(29)

{ac,x(u, v)} = Is] (IV d]-~x + [vd] -bT-yJ × [F]-I {u} f,- d, O[F]

{Bc,y(u,v)} = [S]~ IN

(19)

leads to the following DRBEM vector

O[F]~

l--~-x + [Vd]-~-y j

x IF]-1 {v}

(20)

Here, the matrix [S] is that given in eqn (16). The diagonal matrices [U d] and IV d] contain the dual reciprocity nodal values of u and v respectively. The matrix [F] and its inverse have the same meaning as before, and the vectors {u} and {v} contain the dual reciprocity nodal values of u and v respectively. 2. The pressure terms in the x- and y-momentum equations

x [r]-I {u}

"V(~y ] [F]-l[Vd]O[F]'~JOy × [~]-1 {~}/

0IF]

-i

a

+ (--~--x [F] [V ]0[F]'~0xj x [r]-I {v}

Pressure correction DRBEM solution

-l

5 ITERATIVE SOLUTION METHOD

d 0[4

Resolving the flow field requires iteration due to the nonlinear terms of the momentum and pressure correction equation. An iterative algorithm is outlined below:

[F]-' {v}]

x

O[FI[FI[ _, vd l--~y) O[F] + @~-y

x [Fl-~{u}}

(1) Supply initial values for u, v and p. Solve eqns (31)-(33) with the right hand side set to zero. (2) Solve x-momentum for u. Update u value using relaxation. (3) Solve y-momentum for v. Use updated u value. Update v value using relaxation. (4) Solve pressure equation for p. The convective terms used in pressure equation are calculated using latest u and v values. Update p value using relaxation. (5) Continue iteration with step 2 until solution convergence. (6) Solve the energy transport equation using eqn

(30)

With these DRBEM vectors in hand, the NS and the pressure correction equations can now be cast in standard BEM form as [H]{u} = [G]{qu} + {Bx(u, v,p)}

(31)

[H](u} = [G]{q,} + {By(u, v,p)}

(32)

[Gl{q,} + 2p{Bvc(U , v)}

(33)

[HI{T} = [G]{q} + {Br(u, v, r)}

(34)

[Hl{P}

=

269

where

(34).

1

[O{Bc,~(u,,,)} + {n,,x(t,)}]

(35)

{By(u, v,p)} = -~ [p{Be,y(u, v)} + {Be,j,(p)}]

(36)

{Bx(u,,,,p)}

= ~

1

-

-

Iterative convergence for the pressure correction scheme is measured in the L~ norm. Convergence is achieved when the preset tolerance e is met by all dependent variables, u, v and p. In all cases discussed below, we set the convergence tolerance at e = 10-3.

OCp

{st(u, v, r)} - - - U {Bc,r(U, v, r ) } + k {Br,,(u, v)} (37)

6 NUMERICAL RESULTS

The vectors, {q,}, {qv}, {qp} and {q} in eqns (31)-(34) contain nodal values of the normal derivative of the u and v velocity components and the nodal values of the normal derivative of the pressure p and temperature T respectively. It is noted that the influence matrices [HI and [G] in these equations have been adjusted in this formulation following Partridge et al. 15 in order to account for the dual reciprocity points.

UffiO

ttfU o

ap=

v=O

L fT,

a~v

au --=0

o

o

o

o

o

o

o

o

o

o

o

o

o

o

o

o

o

i 2h '= 3.75 m m o o[ o o

o

o

o

o

o

o

o

o

o

o

o

o

o

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

vffi0

pffiP~

-t

In this section we report results obtained with the pressure correction DRBEM scheme outlined above for the case of laminar channel flow. In particular, we consider the channel configuration shown in Fig. 1. Air at room temperature flows in a l'0m long and 3"75mm wide channel. Inlet pressure is varied to achieve the desired exit Reynolds number. A 1 m/s slug flow enters the

8x

v=O O_..~P= 0 ax

TffiT~ 01 0

0

u=,O

vffiO

~y " "a'7 a2v

8"r=0

0x

r.,. To

1.5 m

Fig. 1. DRBEM discretization for the analysis of flow in a channel.

J "1

270

C. P. Rahaim, A. J. Kassab

1.0 0.8 c-

0.6 0.4

O I-

0.2

o

0

0.. -I UJ Z Z
. I

-0.2 -0.4

I..,~_.~xST==IE:!0E x : , :.O I

~

U/Urn

-0.6 -0.8

-1.0 Fig. 2. DRBEM computed velocity profiles at various axial locations down the channel for Re = 300.

1.0 0.8 0.6 0.4 o

0.2

o

gt. _1 UJ Z Z

0

I

0.2

I

0.6

1

0.8

I

[ -~-

-0.4

-0.6

MORIHARA-CHENG P R E S E N T x = 10 MORIHARA-CHENG

.U'Om

1.0

] [ .,,0[ [~ x = 20 ] I ~ PRESENT x = 2 0 [ I ~ MORIHARA-CHENGx = 220[ [ ~ PRESENT x = 220 ]

-0.2

I

O

/

0.4

I~

x = 10

L.

-0.8 -1.0 Fig. 3. DRBEM computed velocity profiles at various axial locations down the channel for Re = 2000.

. /

o PRESEN'~,

Re=2000

Table 1. N m b e r of iterations for convergence of the pressure correction DRBEM with e = 1 0 - 3

Reynolds number 100

~

~ - . PRESEN'~Re=300 ,

I-

E

Z

EXACT~ SOLUTION

200 300 2000 2300

Number ofiterations 10 12 19 21

~.~

10 ''~'F' ~ ' ; 'k'~" ~ ' : '~'~" ~ ' ;'~'~" 0.0001 0.001 0.01 O.1 X

Fig. 4. DRBEM computed Nusselt number for Re = 300 and Re = 2000 compared with exact solution from Ref. 19.

channel. The exit is assumed to be fully developed and channel wise gradients of u and p are set to zero, while the v component o f velocity is set to zero. The no slip condition is imposed at the top and bottom walls of the channel by setting u and v to zero there. The pressure

boundary condition at the top and bottom walls is obtained from the y-momentum equation. The top and bottom wall temperatures are set at 100°C while the entrance temperature is set at 200°C. In all, 58 linear boundary elements are used to discretize the boundary, with 11 elements distributed on each end. The dual reciprocity points are taken as the boundary element nodes plus the additional 30 internal points shown in the figure. Thus a total of 88 dual reciprocity points were used to expand the right hand side terms and arrive at the RHS convective vectors {Bx(u,v,p)} and {By(u,v,p)}, the pressure correction vector {Bw(u , v) }, and the energy equation vector {Br(u,v, T)} in eqns (31)-(34). There are then 58 unknowns used to resolve this flow field.

Pressure correction D R B E M

271

solution

i 0.0.24

° ! °'11

"

'I'""

I

°1 --

Fig. 5. DRBEM computed temperature profile for Re

=

300 compared with data from Ref. 19.

0.8 0.6

i

0.4 0.2 o

-o.2+

'

'

1.1

1.2

'

1.3

'

1.4

'

T~-~"

1.5 ~'~m

17 TITw

°'t Fig. 6. DRBEM computed temperature profile for Re = 2000 compared with data from Ref. 19.

1.0

I.,-

o=T-T. T,,

i

o

1.0 I

x=5

I

x-30

I

x,,140

Fig. 7. DRBEM computed temperature profiles at various axial locations down the channel for Re = 2300. DRBEM results were obtained for several Reynolds numbers and compared to the finite difference channel data reported by Morihara and Cheng, 17 Bodoia and Osterle is and Shah and Bhatia. 19 In Figs 2 and 3 and 5-7, the axial position, x, is nondimensionalized by the channel half-width, H, in order to be consistent with the data in Refs 17 and 18, and in Fig. 4, the axial

position, x, is nondimensionalized with the product of the hydraulic diameter and the Peclet number, Pe = ReDhPr, to be consistent with the data in Ref. 19. Reynolds numbers are evaluated at the channel exit and are based on the hydraulic diameter, i.e. R e = pUm,e~tDh/lZ. Results are reported for R e = 300, 2000 and 2300. Table 1 provides the number of

272

C. P. Rahaim, A. J. Kassab

iterations to convergence for the pressure correction scheme for several Reynolds numbers. In Figs 2 and 3, the comparison of DRBEM computed velocity profiles for Re = 300 and Re = 2000 with the data reported in Refs 17 and 18, shows agreement in trend and magnitude. The developing flow is reproduced faithfully by the proposed DRBEM scheme, and no oscillations were observed in the flow field. In Fig. 4 the DRBEM computed Nusselt number for Reynolds numbers of 300 and 2000 are plotted and compared with the data of Ref. 19. In Figs 5 and 6 the cross channel temperature distribution for Reynolds numbers of 300 and 2000 are plotted and compared with the data of Ref. 19. The computed values are in good agreement with the data. Finally, DRBEM computed temperature distributions for a Reynolds number of 2300 are shown in Fig. 7 at representative locations down the channel. The numerical verification of our scheme is thus successful for channel flows.

7 CONCLUSIONS An iterative D R B E M pressure correction scheme is developed for the solution of incompressible laminar fluid flow in conjunction with a D R B E M solution of the energy equation. The dual reciprocity method is developed to discretize the nonlinear forcing terms in the Navier-Stokes equations and the pressure correction equation. Results are obtained for flow in a channel under several Reynolds numbers. D R B E M computed velocities, temperatures and Nusselt numbers are in agreement with previously published data.

REFERENCES 1. Mores, P. & Feshbach, H. Methods of Theoretical Physics. McGraw-Hill, New York, 1953. 2. Carslaw, H. S. & Jaeger, J. C. Conduction of Heat in Solids, 2nd edn. Clarendon Press, Oxford, 1959. 3. Shaw, R. P. An integral equation approach to diffusion. Int. J. Heat and Mass Transf., 1973, 17, 693-9. 4. Chang, Y. P., Kang, C. S. &Chen, D. J. The use of fundamental Green's functions for the solution of heat conduction problems in anisotropic media. Int. J. Heat and Mass Transf, 1973, 16, 1905-18. 5. Wrobel, L. C. & Brebbia, C. A. A formulation of the boundary element method for axisymmetric transient heat conduction. Int. J. Heat and Mass Transf, 1981, 24, 843-50.

6. Wrobel, L. C. & Brebbia, C. A. (eels) Boundary Element Methods in Heat Transfer. Computational Mechanics Publications and Elsevier Applied Science, Southampton and London, 1992. 7. Wu, J. C. Boundary elements and viscous flows. In Proceedings of BETECH92, The 7th International Boundary Element Technology Conference, eds C. A. Brebbia & M. S. Ingber. Albuquerque, New Mexico. Computational Mechanics Publications, Boston, MA, 1992, pp. 1-18. 8. Brebbia, C. A. & Skerget, P. Diffusion-convection problem using boundary elements. Advances in Water Resources, 1984, 7, 50-7. 9. Skerget, P., Kuhn, G., Alujevic, A. & Brebbia, C. A. Time dependent transport by BEM. Advances in Water Resources, 1989, 12, 9-20. 10. Wrobel, L. C. & De Figueiredo, D. B. A dual reciprocity boundary element formulation for convection-diffusion problems with variable velocity fields. Engineering Analysis, 1991, 8(6), 312-19. 11. Qiu, Z. H., Wrobel, L. C. & Power, H. An evaluation of boundary element schemes for convection-diffusion problems. In Proceedings of BEM XV, the 15th International Conference on Boundary Elements, eds C. A. Brebbia & J. J. Rencis. Vol. I, Worcester, Massachusetts. Computational Mechanics Publications, Boston, MA, 1993, pp. 25-37. 12. Power, H., Brebbia, C. A. & Ingham, D. B. (eds) Boundary Element Methods in Fluid Dynamics, 11. Computational Mechanics Publications, Boston, MA, 1994. 13. Power, H. Boundary Element Applications in Fluid Mechanics. Computational Mechanics Publications, Boston, MA, 1995. 14. Rahaim, C. P. & Kassab, A. J. Application of the DRBEM to the solution of viscous fluid flow in a channel. In Proceedings of BETECH94, the 9th International Conference on Boundary Element Technology, eds C. A. Brebbia and A. J. Kassab. Orlando, Florida. Computational Mechanics Publications, Boston, MA, 1994, pp. 147-54. 15. Partridge, P. W., Brebbia, C. A. & Wrobel, L. C. The Dual Reciprocity Boundary Element Method. Computational Mechanics, Southampton, 1992. 16. Goldberg, M. A. The numerical evaluation of particular solutions in the BEM--a review. Boundary Element Communications, 1995, 6(3), 99-106. 17. Morihara, H. & Cheng, Ralph a-Shun. Numerical solution of the viscous flow in the entrance region of parallel planes. Journal of Computational Physics, 11, 1973, 550-72. 18. Bodoia, J. R. & Osterle, J. F. Finite difference analysis of plane Poiseuiile and Couette flow developments. Applied Scientific Research Section A, 1961, 10, 265-76. 19. Shah, R. K. & Bhatia, M. S. Laminar convective heat transfer in ducts. In Handbook of Single Phase Heat Transfer. eds S. Kakak, R. K. Shah & W. Aung. Wiley-Interscience, New York, 1987, Chapt. 3, pp. 3-1.13-1.137.