Time dependent transport problems by BEM

Time dependent transport problems by BEM P. Skerget

Department of Mechanical Engineering, University of Maribor, PO Box 224, 62000 Maribor, Yugoslavia G. Kuhn

University of Niirnberg-Erlangen, FRG A. Alujevic

Institute Jozef Stefan, University of Ljubljana, Yugoslavia C. A. Brebbia

Computational Mechanics Institute, Southampton, UK

INTRODUCTION The paper deals with the application of BEM to time dependent energy and momentum transport in laminar flow of an isochoric fluid. The Boussinesq approximation is used to consider the buoyancy effects. Although the Navier-Stokes equations can be solved in terms of physical variables 8'21,22, or in terms of vorticitystream function 9, the velocity-vorticity formulation is used in the paper. Velocity-vorticity formulation was originated by Wu and his co-workers 25-2v, and later further developed by Skerget et al. ~3 19. The classical example of natural convection in closed cavity is studied and the numerical results obtained by BEM are compared with those obtained by finite differences (FDM) and finite elements (FEM).

conditions for the velocity vector and temperature are prescribed. The following sets of boundary and initial conditions can be associated with the above equations v= v

on F for t > t o

b'= ~o

in f~ at t = t o

and

T= T

on F l for t > to

8T 8n

1 ko

1. G O V E R N I N G E Q U A T I O N S

~T

h

The time dependent laminar motion of an incompressible viscous fluid and energy transport is governed by the Navier-Stokes momentum equation, mass and energy conservation laws

~n

ko

p[SF/at+(~.V)~=#AF-Vp+p~

(1)

V" ~= 0

(2)

(3,T/at + V({T) = aA T

(3)

formulated in terms of primitive variables: velocity F, pressure p and temperature T. Material properties, density p, dynamic viscosity #, and thermal diffusivity a = ko/cp (where k o is conductivity, and c specific heat), are assumed here to be constant parameters. Equations (1), (2) and (3) are sufficient for the determination of velocity, pressure and temperature fields, provided appropriate boundary and initial

(4)

T=T o

on F 2 for t > t o

(T-Ts)

on F 3 for t > t o

in ~ a t t = t o

(5)

where h is the heat transfer coefficient and Ts ambient temperature. The buoyancy body force can be approximated with Boussinesq relation, where the temperature influence on density is used only with body forces, while it is neglected in all other terms. The following relation between the fluid density and temperature may be used

p =pod -/~(V- To)]

(6)

where Po is some reference density at temperature To, and fl is the volume coefficient of thermal expansion. With equation (6) the momentum equation (1) can now be written as

Paper accepted November 1988. Discussion ends August 1989.

po[OF/& +(~'V)v-'] = # A { - - V P - p o f l ( T - - To)~

© Computational Mechanics Publications 1989

Adv. Water Resources, 1989, Volume 12, March

(7)

9

Time dependent transport problems by BEM: P. Skerget et al. P being the modified pressure

to the boundary are

P=P-Pog" F

(8)

with ~ = ( x , y , z ) being the position vector and ~= (g~, gy, g=) is the gravity acceleration.

1.1 Vorticity formulation It is convenient to introduce the vorticity variable 26

=Vx{

(9)

to divide the fluid flow computation into its kinematic and kinetic parts. With the following identity for the convective term

({- ¢){= ~v2/2 + ~ x

(1o)

the momentum equation (7) can be given in the form

a~/at+~x

F=vAF-VH-fl(T-To)~

(12)

which is the vorticity transport equation. For the plane motion the vorticity vector has just one component orthogonal to the plane of the flow, and equation (12) reduces to the scalar vorticity equation

aw/~t + V ( { w ) : v A w - f l ( g y , - g x ) V T

(13)

1.2 Vector potential The kinematics of the flow is described by equations (2) and (9). The vector potential of the solenoidal velocity field b~can be introduced by

F=Vx~,

V.~=O

=

vt(S)= {(S). ~(S)

(19)

The relations to velocity v~ and vy components are the following

v~(S) = v,(S)n~(S) - vt(S)ny(S)

vy(S)= v,(S)ny(S) + vdS)n~(S)

(20)

The global mass balance equation is

f ~ . ~(s)dn= f ~(s)~(s) dr =o

(21)

S denotes boundary points and s points inside of the domain.

(11)

where H=P/po+VZ/2 is the total pressure, and v the kinematic viscosity (v=l~/p). It is easy to derive the formulation without pressure by taking the curl on both sides of equation (11)

a~/at+({.V)~-(~'v.V){=vA~-Vxfl(T-To) ~

v.(s) ~'(s). i(s)

(14)

and the kinematic part of computation is represented by the elliptic Poisson's equation for vector

2. BOUNDARY INTEGRAL EQUATIONS FOR THE KINETICS

2.1 BIE of the flow kinetics ( v - w - T ) The vorticity transport equation describes the development of the vorticity field with time. It is a parabolic equation representing an initial boundary value problem. In addition to the initial vorticity distribution, boundary vorticity values at subsequent time levels have to be known. The boundary vorticity values are determined through the kinematic relation between the velocity and vorticity fields. Let us consider inhomogeneous parabolic equation for : the scalar function w(s, t) vAw(s, t) - ~w(s, O/at + b(s, t) = 0

in ~

(22)

in a finite time interval I(to,tv) with corresponding boundary conditions of the first and second kind

w(S, t) = ~(S, t)

on F1

aw(S, t) a~(S, t) an(S) an(S)

on F z

(23)

The initial conditions needed for the solution of equation (22) are which, for plane motion, reduces to the scalar elliptic equation for the stream function A~F= - w

(16)

where the velocity vector is defined with

~= (vx, vy) = (aud/ @, - aud/ ax)

(17)

The unit normal and the tangent to the boundary are given by

Adv. Water Resources, 1989, Volume 12, March

(24)

Body forces are represented by b(s, t) term. Using Green's theorem for scalars or by weighted residual statement the following corresponding BIE can be written 2

tF; S,t) c(Ow(Lt~)+v i f w(S,t). &~*(L ~5~(~S)

dF dt

= v ( (aw(S, t) u*(¢, tF s, t) d r dt J J an(S)

(18)

while the normal and the tangential velocity components

10

in fl

+ffb(s,t).*(¢,t~;s,t)dndt

n(S) = (nx, ny) t(S) = (t x, ty)= ( - ny, nx)

w(s, t) = #o(S, to)

+ fw(s, tF - 1)U*(~, tF ; S, tv_, ) d ~

(25)

Time dependent transport problems by BEM: P. Skerget et al. where u*(~, t~; S, t) is a parabolic fundamental solution, with (4, tr) being source space-time point, and (S, t) or (s, t) are field space-time points on the boundary and in the domain respectively

u*(4, tF ; S, t) = 1/4nvz - e x p [ - r2(4, S)/4vz] &*(4, t~; S, t) = d(~, S)/81rvEz2" e x p [ - r 2(4, S)/4vz] On(S) (26) where z = tF - t and d(~, S) = [x/( 0 - xi(S)]ni(S), while ni is the unit outward normal. The coefficient c(O has values 1.0 (when 4 is in the domain), 0.5 (if 4 lies on a smooth boundary), or ®/2re (when 4 is on a nonsmooth boundary, ® being the internal angle of the boundary at the 4 point)'. Taking the body force term b(s, t) in equation (22) equal to convective and buoyancy terms in equation (13) one can extend the integral formulation to the vorticity parabolic transport equation

[ ~ . . . . au*(g,t~; s,t)

C(OW(4, tF)+ v J J w ( 5 , t)

?~n~

a boundary problem only, described by the first two boundary integrals. The third boundary integral represents the convective vorticity flux through the boundary which vanishes for v,=0, and boundary vorticity generation due to buoyancy effects. The first domain integral is due to force and buoyancy convection effects, while the last domain integral expresses the influence of the initial vorticity field on the development of a new one.

2.2 BIE for the energy transport The energy transport equation (3) is a parabolic diffusion-convection equation for temperature. Owing to the formal similarity between equations (3) and (13) one can formulate the following BIE for temperature 3

c(~)T(4, tF)+af_fT(S, t)

OU*(¢,tF; S,t) On(S) dF dt

/" ~8T(4, t)

=ajJ~n(s)-U*(4, t.; S,t) dr dt -ffT(S,t)v.(S,t)u*(4,tF;S,t)dFdt

d r dt

+ffT(s,t)~(,,t)%*(4,t~;s,t)dndt

~"l'Ow(S, t) • =v | | ~ - , ~ , u (~,tr;S,t)dFdt j d cn~a)

+ iT(s, tF - ~)U*(4, tF ; S, te - , ) df~ J

-ff +fl(gr, --g~)VT(s, t)}u*(4, tF; s, t) dn dt

+ lw(s, tr- 1)u*(~, tF; S, t r _ ,) dr) d

(27)

The domain integral comprises the derivatives of the vorticity and temperature at domain points. The problem can be solved in different ways. The derivatives may be computed explicitly by derivation of the above equation and BIE for the temperature, or can be derived in a similar manner as in FEM by means of interpolation functions 9. An alternative way is to eliminate the derivatives using Gaussian divergence theorem on domain integral resulting in the following BIE ~5

C(OW(4'tF)+Vff w(S't)Su*(4'tF;S't)dFdtOn(s)

:

v jf ji gwf ~,,l u*(4, tF ; S, t) dF dt

(29)

The BIE (29) expresses the energy transport in the integral form. The energy diffusion is a boundary problem only, expressed in the first two boundary integrals. The third boundary integral describes the convective energy flux. The first domain integral is due to convective effect, while the second one represents the contribution of initial temperature field on the development of a new one. 3. B O U N D A R Y E L E M E N T A P P R O X I M A T I O N F O R T H E KINETICS

In order to obtain approximate numerical solutions of the velocity, vorticity and temperature fields, one has to discretize the boundary into E boundary elements with N e boundary nodes, and the domain into C cells with N c internal points. The values of the variables T, 8T/Sn, w, v, and the products wF and Tv" are assumed to vary within each element or cell and each time step according to the space {~} and time interpolation function {~}, i.e. for the boundary

w(,, t)= {'I'}~{V}{w}~.,

TO1, t)= {O}T{~P}{T}~,

- j)[w(S, t)v.(S, t) +fig, T(S, t)]u*(4, tF ; S, t) dF dt

(,. t)= {O}r{V} +.l-.[[w(s,t)F(s,t)+fl(gr-9x)T(s,t)]Vu*(~,tr; s,t) d ~ dt + fw(s, tF- 1)U*(~, tF; S, tr - , ) df~

(28)

where gt=O.f is tangential gravity acceleration component to the boundary. The BIE (28) in a clear manner describes the transport of vorticity phenomenon in time. The diffusion process is

m' wn.t,, t)= {¢}~{V}{wv.}~.

Tv,(,, t)-- {O}~{V}{wv,}~.

(30)

and for the domain

WV,(~,, 1"]2, t)= {(1)}T{klI}{WUi}nm,

w(.,,.2)-- {o}T{w}"

Tvi(rh , fix, t)= {(I)}r{~} {Tvi}~n,

TI,,,,2)={*}T{T}" (31)

Adv. Water Resources, 1989, Volume 12, March

11

Time dependent transport problems by BEM: P. Skerget et al. where q and ~l,q2 are homogeneous coordinates. The super-index n refers to the number of nodes within each element or cell, associated with function nodal values, while the sub-index m refers to the degree of variation of the function {~}, i.e. m = 1 i f ~ is constant, m = 1,2 i f ~ is linear, etc. Let us assume that the field functions remain constant (W= 1) within each time increment (Z=tF--tF-~). The time integrals involved in the discretized equations can be evaluated analytically as (where the integration limits are from t F l to tF)

The equation is first applied to all boundary nodes

--1 ([6]{wv. +,g,v}~ - [Dx] {wvx+/3g,V}~ V

- [Dr]{wvr-flgxT}F - [B] {W}F - 1 )

_

U* = v f u * (~, te; S, t) dt= 1/4n" E a(xr_

1)

aU*

(au*(~, t~,' S, t) dt = d( ~, S)/2nr2( ~, S) "e x p ( -

au*. =v

(au*(~, ~ tr; s't) dt

Ox~

J

xr-

1)

(37)

In general the proper boundary conditions are defined by velocity vector {(S)=(Vx, Vy) or with normal and tangent velocities (v,, vt) to the whole boundary. The only unknowns in the equation (37) are then boundary vorticity fluxes, while the boundary vorticity values are obtained in the kinematic part of the computation. Vorticity values computation at any internal points follows in an explicit manner from (36) for c(~)= 1.

3.2 Discretization of the energy equation

cxi

= [xi(~) - xi(s)]/2nr2(~, s) • exp( - x r - 1)

(32)

The application of the collocation method to the equation (29) leads to the following discretized formulation n

where E1 is the exponential integral function defined by

c(#)T(~, tr)+

{h}T{T}~ = e=l

El(x)=

(e-'/t)dt

(33)

--a

{g}T e=l

F

{g}r{Tv.}~ - E {dl}T{Tv,}~ - Z {b}V{T} "F-' e

1

c=l

c=l

(38)

and the argument XF-~ is XF - 1 = r 2 ( ~ , S)/4VT

(34)

The following space integrals are involved in the discretized BIE fO" 8U*

h"~=j ~ n d r ,

Equation (38) is first applied to all boundary nodes to compute unknown temperature and heat flux values [H]{T}F

= [G] 1

~n

F

-- -- ([G] { Tv.}F -

g"e= f o " u * dF,

a

[D ~]{ Tv,}r -

[ B ] { T}F - , ) (39)

d:'=f *" ~x/°v* d~,

f

b~ = v O"u*_, dr2

(35)

which are functions of the geometry, time increment and material property.

3.1 Discretization of the flow kinetics ( v - w - T ) Boundary element discretization of BIE for parabolic transport equation can be readily obtained from the equation (28) as follows E

c(Ow(~, try)+ Z {h}~{w}~ e--I

- L {g} flaw)"

[A]{X}= {F} + {FN(T, v)}

where {X} stands for unknown boundary temperatures or energy fluxes. It is easy to introduce the Cauchy type boundary conditions into the formulation 1s. An explicit computation of domain temperatures follows from (38) for c(O= 1.

1

4. B O U N D A R Y I N T E G R A L E Q U A T I O N S FOR THE KINEMATICS

Vc=l

4.1 BIE for the stream function

+- L {dx}r{wv~+flgyT}~ 1

(40)

1 E

-e:,

The kinematic relation between velocity and vorticity fields is given by inhomogeneous elliptic equation for scalar function U?(s)

1 c

+v ~:,L {dr}r{wvr-flg~T}~ +v ~=,: {b}T{w}7~-, (36)

12

where i= 1,2 for plane problems. After rearranging columns and rows due to the application of boundary conditions, the above system can be written as

Adv. Water Resources, 1989, Volume 12, March

AUd(s)+b(s)=O

in f~

(41)

Time dependent transport problems by BEM: P. Skerget et al. with Dirichlet's and Neumann's boundary conditions

The term b(s) represents body forces. It is identically zero for potential flow while for the rotational flow represents the vorticity field. Using Green's theorem for scalars, one can write the following B1E

the flow or the relation between velocity and vorticity fields at the same instant of time in an integral form. The boundary integrals contain the velocity boundary conditions, the domain integral represents the contribution of the vorticity field on the development of the velocity field, while the potential effect of the flow is incorporated in the boundary integrals. The equation allows for explicit computation of the velocity components in the domain. For the plane problems, the vector equation (46) simplifies to two scalar equations

c(¢)~P(~)+ f ~P(S)0u*(~, On(S)S) dF

(~, S) dF c(Ov~(O+ f vx(S) 0u* On(S)

~P(S)= q~(S) oo,,(s)

-

on F 1

oae(s)

On(S) On(S)

= j -O°?(S) ~) u*(¢, S) dr +

-

ft(S)

on F 2

fw(s)u*(~, s) df~

(42)

(43)

where u*(~, S) is an elliptic fundamental solution, with being source point, while S and s are field points on the boundary and in the domain respectively

C .s. OU*(~,S) fw(s)~df~ = Jvrt ) Ot(S)- d r - . ) e'yts) c(~)vy(~)+f vy(S)0u*(~, On(S)S) dF = - f vx(S)Ou*(~,s) &(S) d._l + J w("t s,,) ~Ou*(~,s) dfl

(47)

u*(~, S)= 1/2n" ln[ro/r(~, S)] Ou*(~,S)_ Vu*(~, S)B(S) = 1/2n "d(~, S)/r2(~, S)

5. B O U N D A R Y ELEMENT FOR THE KINEMATICS

On(S)

0u*(~, S) _ Vu*(~, S)I (S) = 1/2n" d,(~, S)/r2(~, S)

Ot(S)

(44)

The potential flow is incorporated in the boundary integrals while the domain integral represents the effect of the vorticity field on the development of the velocity field. Velocity components in the domain can be obtained explicitly by the derivatives of the equation (43)

0oe(¢)

F. . . .

o

FOu*(~,s)] ,~

Let us divide the boundary into E boundary elements with N e boundary nodes, and the domain into C cells with N c internal points. The field functions b" and w are approximated within each element and cell according to the space shape functions multiplied by n nodal values

vi(q)={~}r{vi} ", i=1,2, w(ql,~12)={~P}r{w}"

(48)

The following integrals which are functions of the geometry only are involved in the discretized BIE he"= Jf o"

~0°?(S) 0 u * ( ~ , S ) dF

APPROXIMATION

Ou*~nd r ,

he",= f ~"

dr, 0u*~4_

dT,= f ~" ~Ou* df~

(" , , 0u*(~,s)

(49) (45)

5.1 Discretization for the internal flows The kinematic part of the fluid flow computation is described by

4.2 BIE for the velocity potential Since each component of the velocity vector t~ satisfies the Poisson's equation (41 k an integral representation for is readily obtainable from (43) by substituting • with

E

c(~)vx(~)+ Z {h}r{Vx}" e=l

,P.

E

An alternative kinematic integral formulation can be derived using Green's theorem for vectors 25, resulting in the following equation ~3 {, c(~)~(~) + j[Vu*(~, S)~(S)] b~(S)dF

C

= Z {h,}~{H "- Z {d,}~{w} ° e=l

(50)

c-1 E

c(~)vr(~)+ 2 {h}r{vy}" e=l E

=-

e=l

(46) where x denotes the cross product. The equation (46) is equivalent to the continuity and the vorticity equations, and expresses the kinematics of

C

Z (h,}~("x} "+ Z {d~}'{w}"

(51)

c=l

Applying equations (50) and (51) to all boundary nodes the following 2N e matrix systems are obtained [c(¢)]{vx(~)} + [HI {Vx}= [Ht]{vr} - [D,] {w}

(52)

[c(~)]{vr(~)} + [/4] {vy}= - [nt]{Vx} + [Dx]{w}

(53)

Adv. Water Resources, 1989, Volume 12, March 13

Time dependent transport problems by BEM: P. Skerget et al. or

[H]

{Vx}=

[Hi] {vr} - [Dr] {w}

(54)

ENERGY PART Solve system of equations for unknown boundary temperatures or heat fluxes Explicit calculation of internal temperatures -

[/4] {vr} = - [/-/,] {Vx} +

[Ox]{W}

~55)

A system of N e equations for N e unknown boundary values can be readily derived by multiplying the above two systems of equations by the unit tangent t'(~)= (t 1, t2) and by summing up the equations ([HI1

+

[Ht]2){Vx} = ([Ht] I -- [H]2){Uy} -([or],

-

[DJ~){w}

(56)

KINETIC PART Solve system of equations for unknown boundary vorticity fluxes Explicit calculation of internal vorticity values The reversed sweep is used in computation of internal temperature and vorticity values. With this, the convergence of the solution is accelerated and greater stability of the iterative scheme thus achieved.

This is identical to apply the tangential velocity component vt(~) = v'(~)~(~)

(57)

at all boundary nodes (~ = Z, Ne), W u 26. If the velocity boundary conditions b~ are prescribed on the whole boundary, then the only N e unknows in equation (56) are boundary vorticity values ([Dr] ' - [Dx]E)r {w}r = - ([/4],

+ [/4,]~){Vx}

- ([Dr], -

+ ([~,],

- [n]~){v,}

[Dx]~)~{w}o

(58)

The boundary vorticity fluxes are then computed in the kinetic part of the computation (37). The other boundary conditions, such as symmetry conditions given by zero vorticity flux #w/t3n = 0 and vanishing normal velocity v, = 0, can be easily introduced in the formulation. Using relation (20) it follows from equation (56) [E]{v,} = [F]{vt}- ( [ D y ] , - [Dx]2){w }

(59)

In this part of the boundary the only unknown is the tangential velocity vt, while the boundary vorticity values are computed in the kinetic part by the equation (37). The velocity vector is calculated from the relation (20). Velocity components in the domain are then obtained explicitly for c(~)= 1 as

7. N U M E R I C A L EXAMPLES

7.1 Thermally driven cavity flow The pure buoyancy driven natural convection is considered. This problem has been proposed by De Vahl Davis et al. 6 as a standard example for comparing different numerical techniques in computational fluid mechanics. Various uniform meshes of linear boundary elements and linear triangular cells are used: (a) M = 10 x 10, (b) M = 20 × 20, (c) M = 30 x 30to determine the temperature and velocity fields. The motion is caused by buoyancy force due to the heated left wall subject to a unit temperature, while the right wall is kept at zero temperature. The rest of the boundary is adiabatic, Fig. 1. The motion is analysed for Rayleigh numbers in the range of 1 0 3 < R a < 106 for the Prandtl number Pr= 1. Fig. 2 summarizes the steady state solutions. One large circulation cell is predicted for Ra<104 with the cell strength increasing with Ra number. At greater Rayleigh numbers Ra > 105 two circulation cells are created. Time dependent case is studied for Ra=105 with time step z=0.1. Figs 3 and 4 present the velocity fields and temperature contours of the time dependent phenomenon at time instants t = 0.2, 0.8, 1.2, 3.2, 6.4 and 40 secs. The solution converges to the steady state solution.

{Vx(~)} = - [/]] {Vx}+ [Ht]{v r} - [Dy]{w}

(60)

Y

{Vr(~)} = - [/4] {vr} - [HJ{vx} + IDol{w}

(61)

1

q=O To=O

6. S O L U T I O N P R O C E D U R E Kinematic relations, vorticity and temperature kinetic equations are coupled in the set of the nonlinear equations. Iterative point under-relaxation solution procedure 23, has to be employed to compute new domain temperature and vorticity values at each time step. The solution requires the following steps at each time increment:

T=O

T=I

g

Start with some initial values for vorticity and temperature K I N E M A T I C PART Solve the system for boundary vorticity values Calculate domain velocity components in an explicit manner

14

Adv. Water Resources, 1989, Volume 12, March

0 I

q=0

1 !

X ---

Fig. 1. Definition of thermally driven cavity flow

Time dependent transport problems by B E M " P. Skerget et al.

/

/

./

/

/

/

I

/

/ ' - ' ' \

I

I

I

"/

,

\ \

\

\

\

.

,

"

-

~ .

i

"

~

.

\

\

/

I

,

,

/

"

\

I

/

/

/

/

/

-I"

/

/

~

.

I

/

II

/117

same order of the accuracy at all boundary nodes. The results obtained by BEM are compared with the numerical results obtained by F D M and F E M . The BEM results for the average Nusselt numbers are shown in Table 1. Three different uniform meshes of the first order linear boundary elements and linear triangular cells are employed in the range of 1 0 3 < R a < 1 0 6 and P r = 1. In the Table 2 De Vahl Davis 5'6 original solutions and bench m a r k solutions obtained with finite difference method are collected for the same Rayleigh number range and Pr = 0.71. The stream function-vorticity formulation was used and different uniform F D M meshes were employed. In the Table 3 the finite element result by Stevens 2° are shown. Again the stream function-vorticity formulation was taken to determine the average Nusselt number. The triangular Lagrangian finite elements with from first to fourth order interpolation functions were considered but mainly the third order elements are recommended for the greater Rayleigh numbers. The different uniform and nonuniform meshes are studied and generally the meshes were refined towards the boundaries. Comparison of the boundary element solutions with those obtained by finite difference or finite element methods shows reasonable good agreement. In general, better results were obtained using BEM comparing with F D M solution for the same meshes. One of the reasons is that the heat fluxes are part of the b o u n d a r y element solution, and the second one is due to that the boundary vorticity values are included in the kinematic integral

Table 1. Mean Nusselt number (by BEM) Pr= 1.0 Ra

i 1 ~

%%

.

.

.

.

.

.

.

.

Mesh

103

104

11 x 11 21 x21 31 ×31

1.113

2.274

105

106

4.677 9.102

.

Table 2. Mean Nusselt number (by FDM) Pr=0.71 i

i

•

-

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

•

i

•

•

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

. .

. .

Ra

~'-,,\\

\

\_

'

Mesh

103

104

105

106

11 x 11 21 x 21 41 x 41 61 x 61 81 x 81

1.105 1.113 1.116

2.307 2.255 2.242

4.767 4.716 4.564 4.531 4.523

6.790 9.502 9.270 9.035 8.928

1.117

2.238

4.509

8.817

bench mark solution

Fig. 2. Thermally driven cavity flow - steady state. Velocity fields .]'or Ra = 10 4, 10 5, 10 6 Table 3. Mean Nusselt number (by FEM) Pr=0.71 Ra The c o m p u t a t i o n of the average Nusselt number defined by

Nu o =

dy 0

represents no extra work, because the b o u n d a r y solutions provide both the temperatures and the heat fluxes of the

Mesh

103

104

13 x 13 19 x 19 19 x 19" 28 x 19" 28 x 19"

1.1136 1.1155 1.1165

2.222 2.232

*

2.241 2.242

105

l 06

4.497 4.512

8.773 8.769

Non-uniform mesh

Adv. Water Resources, 1989, Volume 12, March

15

Time dependent transport problems by BEM: P. Skerget et al.

i" •~. \

\

~

\ \

-

\

\

-~..,. . . . . . . . . . . .

\'-.

\

k,,

I

!

I I I /

I I !

/

/

/

/

l | l l S O S o o - - . r

l

,

l

,

,

#

'

.

¢

.

.

.

.

.

.

.

.

.

.

.

.

g ¢ ~ * .

.

.

.

.

.

.

.

*

*

.

.

e

t

.

*

t

l

.

.

'" . . ' / / / f / - /

...

55

55 ; .

/[

Fig. 3.

Thermally driven cavity flow - transient states. Velocity fields for Ra = 105 at t = 0.2, 0.8, 1.2, 3.2, 6.4, 40 secs

formulation and so there is no need to employ approximate formulae to define vorticity boundary values. The Fig. 5 compares the bench mark solutions with BEM results for the average Nusselt number.

7.2 Natural convection in a horizontal channel Benard's cellular flow caused by the natural convection

16

Adv. Water Resources, 1989, Volume 12, March

due to temperature differences from the heated bottom to the cooled ceiling in a horizontal channel has also been evaluated, where any disturbance causes convective motion. It has been found by experiments that the cellular flow does exist up to Ra = 47 000. In order to demonstrate the same effect by the numerical evaluation, a channel length has been selected

Time dependent transport problems by BEM: P. Skerget et al.

r

Fig. 4.

Thermally driven cavity flow - transient states. Isotherms for Ra = 105 at t = 0.2, 0.8, 1.2, 3.2, 6.4, 40 secs

five times larger than the height, while the third dimension is considered great enough to ensure the plane state. The bottom has been at T = 1.0 and the celing at T=0.0, while the left and right sides were adiabatic. Velocities at all boundaries have been assumed zero. Material properties were selected to yield Gr= 104 and Pr= 1.0, i.e. R a = 104.

The domain has been discretized by 40 isoparametric quadratic elements and 100 internal cells (9 point Lagrange), producing 441 nodes, 80 of them being at the boundary. Time step selected was 0.2 secs. Results demonstrate the development of velocities (Fig. 6), vorticity (Fig. 7) and temperatures (Fig. 8) for t = 0.2, 6.4 and 20 secs. The development of the cellular flow can

Adv. Water Resources, 1989, Volume 12, March

17

Time dependent transport problems by BEM: P. Skerget et al. 1o g ~ l

'

1

;

|

||,w,

I

~

|

,

i |l|,

I

1

,

,

i

~

g4

z

1/ 10 3

10 4

10 5

10 6

Royleigh number, Ro

Fig. 5.

Average Nusselt number

I ---

\ \

I

i

..--

f

,,.--

f

f

f

"-----

J

Fig. 6.

T h e r m a l l y driven cellular f l o w - velocity at t = 0.2, 6.4, 20 secs

be observed, which starts from the lower corners. At time t = 20 secs there is practically a steady state reached, with six circulation cells developed. 8. CONCLUSION Boundary element method is applied to the numerical simulation of plane fluid flow and energy transport

18

Adv. Water Resources, 1989, Volume 12, March

problem (thermally driven cavity flows). The buoyancy force is included using Boussinesq approximation. The vorticity-velocity (v-w) formulation is used to solve the fluid motion problem. Introducing the vorticity variable, the computation of the problem is partitioned on its kinematics and kinetics parts. The behaviour and the physical meanings of the different terms in the integral equations are discussed. Due to the fundamental

Time dependent transport problems by BEM: P. Skerget et al.

~ Fig. 7.

!

~

~

~

~.~

~. ,, ,.

Thermally driven cellular flow - vorticity at t = 0 . 2 , 6.4, 20 secs

I

I

I

__..__._.____

Fig. 8.

Thermally driven cellular flow - isotherms at t = 0.2, 6.4, 20 secs

Adv. Water Resources, 1989, Volume 12, March

19

Time dependent transport problems by BEM: P. Skerget et al. solution, a part of the transport mechanism is transferred to the boundary producing a very stable numerical scheme. In general, the boundary integral equations behave well when the source point coincides with the boundary nodal points. The problem can arise with the integral equations involving the tangential derivatives of the fundamental solution when in some cases the source point has not to coincide with the boundary nodes x2'24. Comparison of BEM results with those obtained by FDM or FEM shows good agreement. To improve BEM solution it is desirable to implement higher order elements, and the corner effects have to be considered with special care.

12

REFERENCES

18

1 2 3

4

5 6 7 8 9

10 11

20

Brebbia, C. A. The Boundary Element Method for Engineers, Pentech Press, London, Halstead Press, New York, 1978 Brebbia, C. A., Telles, J. C. F. and Wrobel, L. C. Boundary Element Methods - Theory and Applications, Springer-Verlag, Berlin, 1984 Brebbia, C. A. and Skerget, P. Time Dependent Diffusion Convection Problems Using Boundary Elements, 3rd Int. Conf. on Numerical Methods in Laminar and Turbulent Flows, Seattle, 1983 Betts, P. L. and Haroutunian, V. A. A Stream Function Finite Element Solution for Two-Dimensional Natural Convection with Accurate Representation of Nusselt Number Variations near a Corner, Int. J. Num. Meth. Fluids, 1983, 3, 605-622 De Vahl Davis, G. and Jones, I. P. Natural Convection in a Square Cavity - a Comparison Exercise, Int. J. Num. Meth. Fluids, 1983, 3, 227 248 De Vahl Davis, G. Natural Convection of Air in a Square Cavity a Bench Mark Numerical Solution, Int. J. Num. Meth. Fluids, 1983, 3, 249 264 Fusegi, T., Faruk, B. and Ball, K. S. Mixed-Convection Flows within a Horizontal Concentric Annulus with a Heated Rotating Inner Cylinder, Num. Heat Transfer, 1986, 9, 591-604 Kitagawa, K., Brebbia, C. A. and Wrobel, L. C. Boundary Element Analysis of Viscous Flow by Penalty Function Formulation, Engineering Analysis, 1986, 3(4), 194~200 Onishi, K.. Kuroki, T. and Tanaka, M. Boundary Element Method for Laminar Viscous Flow and Convective Diffusion Problems, Chapter 8, Topics in Boundary Element Research, (Ed. C. A. Brebbia), Springer-Verlag, Berlin, 1985, 2, 209-229 Rizk, Y. An Integral Representation Approach for Time Dependent Viscous Flows, PhD Thesis, Georgia Institute of Technology, 1980 Schnipke, R. J. and Rice, J. G. A Finite Element Method for Free and Forced Convection Heat Transfer, Int. J. Num. Meth. Eng., 1987, 24, 117 128

A d v . W a t e r R e s o u r c e s , 1989, Volume 12, M a r c h

13

14 15 16 17

19

20 21

22

23

24

25

26 27

Schmidt, G. On Spline Collocation for Singular Integral Equations, Math. Nachrichten, 1983, 111, 177 196 Skerget, P., Alujevic, A. and Brebbia, C. A. The Solution of Navier-Stokes Equations in Terms of Vorticity-Velocity Variables by Boundary Elements, 6th Int. Conf. on Boundary Element Method, QE2, Southampton New York, SpringerVerlag, Berlin, 1984 Skerget, P., Alujevic, A. and Brebbia, C. A. Analysis of Laminar Fluid Flows by Boundary Elements, 4th Int, Conf. on Numerical Methods in Laminar and Turbulent Flows, Swansea, 1985 Skerget, P. and Brebbia, C. A. The Solution of Convective Problems in Laminar Flow, 5th Int. Conf. on Boundary Element Method, Hiroshima, Springer-Verlag, Berlin, 1983 Skerget, P., Alujevic, A. and Brebbia, C. A. Vorticity-VelocityPressure Boundary Integral Formulation, 8th Int. Conf. on Boundary Element Method, Tokyo, Springer-Verlag, 1986 Skerget, P., Alujevic, A., Brebbia, C. A. and Kuhn, G. Boundary Elements for Mixed-Convection Flow Problems, 3rd Intl Conf. on Boundary Element Technology BETECH, Rio de Janeiro, Springer-Verlag, Berlin, 1987 Skerget, P., Alujevic, A. and Brebbia, C. A. Nonlinear Transport Problems in Fluids by Boundary Elements, 3rd Int. Conf. on Numerical Methods in Nonlinear Problems, Dubrovnik, 1986 Skerget, P., Alujevic, A., Kuhn, G. and Brebbia, C. A. Natural Convection Flow Problems by BEM, 9th Int. Conf. on Boundary Element Method, Stuttgart, Springer-Verlag, Berlin, 1987 Stevens, W. N. R. Finite Element Stream Function-Vorticity Solution of Steady Laminar Natural Convection, Int. J. Num. Meth. Fluids, 1982, 2, 349 366 Tosaka, N. and Kakuda, K. Numerical Simulation for Incompressible Viscous Flow Problems Using the Integral Equation Methods, 8th Int. Conf. on Boundary Element Method, Tokyo, Springer-Verlag, Berlin, 1986 Tosaka, N. and Fukushima, N. Integral Equation Analysis of Laminar Natural Convection Problems, 8th Int. Conf. on Boundary Element Method, Tokyo, Springer-Verlag, Berlin, 1986 Wahbah, M. M. Computation of Internal Flows with Arbitrary Boundaries Using the Integral Representation Method, School of Aerospace Engineering, Georgia Inst. of Technology, DAAG 29-75-G-0147, 1975 Wendland, W. L. Asymptotic Accuracy and Convergence for Point Collocation Methods, Topics in Boundary Element Research, (Ed C. A. Brebbia), Springer-Verlag, Berlin, 1985, 2, 230-250 Wu, J. C. and Thompson, J. F. Numerical Solution of Time Dependent Incompressible Navier-Stokes Equations Using an Integro-Differential Formulation, Computers and Fluids, 1973, 1, 197 215 Wu, J. C. Problems of General Viscous Flow, Developments m Boundary Element Method, Vol. 2, Ch. 2, Elsevier Applied Sci. Publ., London, 1982 Wu, J. C., Rizk, Y. M. and Sankar, N. L. Problems of TimeDependent Navier-Stokes Flow, Developments in Boundary Element Methods, Vol. 3, Ch. 6, Elsevier Applied Sci. Publ., London, 1984