Numerical modeling of three-dimensional flow and pollutant dispersion around structures

Environmental Software9 (1994)

175 187 1994 Elsevier Science Limited Printed in Great Britain. 0266-9838/94/$7.00

ELSEVIER

Numerical modeling of three-dimensional flow and pollutant dispersion around structures R o b e r t L. L e e Regional Atmospheric Sciences Division, Lawrence Livermore National Laboratory, P.O. Box 808, Livermore, California 94551, USA

(Received 13 July 1993; revised version received 19 January 1994; accepted 5 May 1994)

ABSTRACT

This paper describes a numerical modeling approach that can be used to provide estimates of emissions at industrial sites. In particular the models presented are capable of simulating the wind flow and dispersion of airborne pollutants around surface-mounted structures such as buildings or building complexes. The calculational procedure in this approach consists of two sequential steps, namely: (i) prediction of the mean flow via a turbulent flow model; and, (ii) employment of the calculated flow field to drive a particle-in-cell wansport and diffusion model. A benchmark simulation is performed in which numerical results from the flow model are compared with other numerical models and with experimental data for flow over a backward-facing step. Results from three-dimensional simulations of flow and dispersion over a two-building complex are also presented KEY WORDS

Pollutant dispersion, turbulent flow around buildings, finite element model, particle-in-cell technique. SOFTWARE AVAILABILITY Name of software

FEMTKE (Finite Element Flow Code) ADPIC (Particle-in-Cell Dispersion Code)

Developers

Robert L. Lee (FEMTKE) Rolf Lange (ADPIC)

Contact address

Lawrence Livermore National Laboratory P. O. Box 808 Livermore, CA 94550

Hardware required

CRAY YMP (FEMTKE) DEC VAX running VMS (ADPIC)

Program language

Fortran

Availability & cost

FEMTKE availability and cost is undetermined at this time ADPIC available at norminal cost from: Energy Scientific and Technology Software Center P. O. Box 1020 Oak Ridge, TN 37831-1020 175

176

R. L. Lee

INTRODUCTION Perhaps the most widely used models for estimating pollutant transport and diffusion are Gaussian-type dispersion models (Slade, 1968). They are simple to use and there are many years of experience supporting these modeling techniques. The development of Gaussian-type dispersion models has reached a level of sophistication such that they are now routinely used as assessment tools for regulatory agencies as well as for many industries. Although these simplified models can be applied with reasonable confidence to pollutant transport within unidirectional flows (such as those over relatively flat terrain), they are less useful for situations where the flows are more complicated. For example, flows over complex terrain or separated flows around obstructions and building wakes are situations where Gaussian dispersion models cannot be applied. In this paper, we propose a more physically-based approach for predicting transport and diffusion for this class of flows. Over the past years our research has focused on the development of finite-element-based numerical models to simulate three-dimensional spills of dense gases and atmospheric boundary layer flows over complex terrain. In addition, we have developed an operational emergency response capability by which a time-varying wind field can he used to drive a particle-in-cell-based transport and diffusion model to predict the fate of pollutants that are released into the atmosphere. The maturity of these numerical tools is at a level that, with modifications, they can be applied to problems associated with building wakes and stack heights as they influence pollutant transport (Lee & Leone, 1991). It is well known that atmospheric flows around surface-mounted obstacles exhibit high Reynolds number phenomena such as flow separations and reattachments, turbulent wakes and transient vortical motions. To capture the turbulent motion associated with these flows we employ a model based on the Navier-Stokes equations with two prediction equations for the turbulent kinetic energy (k) and viscous dissipation (e). This turbulence model, known as the k-e model, has been used extensively and has a good reputation for reliability. Moreover, it is computationaUy economical compared to other more complex models (Ferziger, 1990). We are developing, in parallel, more complex turbulence models that are based on an algebraic-stress formulation and another based on a large eddy formulation. Both models are in early stages of evolvement and will be discussed in future publications. The flow model that we have developed employs a modified finite element method to solve the Navier-Stokes equations of interest. This numerical solution technique has been used with good success in a variety of atmospheric-related applications such as heavier-than-air gas dispersion (Chan, et al., 1984). and boundary layer

flow over complex terrain (Lee & Leone, 1988). In the present model, we include two additional equations for the turbulence field variables and employ blocks of "solid" elements to represent building structures. The transport and dispersion of pollutants are simulated via a particle-in-ceU technique. In this approach the pollutants are represented by inert particles with fixed masses each of which are advected and diffused by advancing the particles with a Lagrangian algorithm. The "advective" velocity of each particle is interpolated from an Eulerian grid where the velocity field is defined. Our calculational procedure consists of two steps. We first compute the flow field by solving the Navier-Stokes equations. The resulting velocity field is then interpolated onto an appropriate Eulerian grid and is used to advance the position of the particles to determine the concentration field. We begin by briefly describing the Navier-Stokes model and the finite element procedure used to obtain the flow field. We then describe the key aspects of the particle-in-cell model that are used to simulate transport and diffusion. To evaluate the performance of the turbulent flow model we performed calculations for flow over a backward-facing step and compared our results with those from other numerical models and with experimental data. Finally we present results for flow around a two-building complex and the dispersion of pollutants from point sources that are located in the vicinity of the buildings. THE FLOW MODEL Most of the numerical models that have been developed for simulating turbulent flows over surface-mounted obstacles are based on finite volume or finite difference techniques (Murakami, 1990; Paterson & Apelt, 1990). In contrast, the Reynolds-averaged flow model that we are developing (FEMTKE) uses a modified finite element method (FEM). The finite element technique provides great flexibility in the arbitrary grading of the mesh and is known to exhibit good accuracy on reasonably distorted meshes. The Model Equations

The partial differential equations that are employed in the flow model are the incompressible, Reynolds-averaged, Navier-Stokes equations with turbulence modeled via a ke formulation. These equations, expressed in cartesian tensor notation, are given by:

J% a% _- o,

%

/

%)' (2)

Three-dimensional flow and pollutant dispersion around structures

ae, ae. ~ae,

--

. o. .,

%tok%) 't% a,,)%

a

'N-- a',/o.a',) ', 'tN -q¥,

(4)

where U~ is the i~ component of the velocity vector, k and are the turbulent kinetic energy and viscous dissipation, respectively. In the above equations, the "pressure" P, Kt, and turbulent diffusivity v, are defined by:

P - - t'+--2/c . x , - - v + v , p 3

, v,--

~

(s)

with p being the dynamic pressure, p the density, and v the laminar diffusivity. The model constants, c~, c~, c2, crk, and o, are set to 0.09, 1.44, 1.92, 1.0, and 1.3 respectively, as suggested in the literature (Rodi, 1980).

177

where U is a vector containing all velocity components, M, K, and N are the lumped-mass diffusion and advection matrices, C is the pressure gradient matrix with its transpose Cr denoting the divergence matrix, and the remaining vectors represent the various source terms. The subscript "s" denotes an appropriate submauix associated with the scalar variables. More details of the formulation, with the form of several of these matrices, can be found in Gresho, et al. (1984).

The Time-Integration Scheme The time-integration scheme used for integrating the ordinary differential equations is (basically) the explicit forward Euler method, except for the pressure, which must be computed implicitly. The pressure solutions are calculated by performing a LU factorization of the pressure matrix, CrM'IC, once and invoking a forwardreduction and back-substitution with a new right-hand-side whenever a new pressure is required. At each timestep "n", the variables are updated sequentially in the following order:

Fhe Finite Element Procedure Equations (1) through (4) are spatially discretized by the finite dement technique in conjunction with the Galerkin method of weighted residuals. The dependent variables are approximated as:

k.. 1 - - k . + A t M ~ t { ( f ~ ) . - IN(U.) + ( g ~ ) . ] k m} (11)

ell

(1 + A t e l - - ) e . d

-- e. + AtM,-t{ff').

~ ¢,(x,)f,(O with/= ej,~.ore,

/(xt,0 =

k=-I

- [N(U.) + (/~,),]e,} (12)

(6)

m

A. = M - l l ( f ' ) .

J'(x.0 : ~ **(x~)J'k(0,

- [N(Ua) + Kff]U~}

(13)

k=-I

where in the discretized domain there are n nodes for the prognostic variables (Ui, k and e) and m nodes for the diagnostic variable P. The basis functions 0k(X) and Wk(x) are piecewise polynomials. For simplicity, ~ is chosen to be tri-linear (in three dimensions) and V is piecewise constant. A coupled set of first-order ordinary differential equations is generated by premultiplying each prognostic equation with the appropriate basis functions and integrating the pressure terms by parts. For numerical convenience, the continuity equation (F_xl. 4) is replaced with a consistent Poisson equation for pressure (Eq.7). Thus the resulting system of matrix equations that has to be solved are given by: MO

-- - C P -

IN(t0 + X ' ] U +/',

(cr~-Ic)e -- crM-z{-[go)

+ x ' ] e + f ' },

(7)

(s)

M./~ = - [/Y (O) + K:]k + Sk - M , e ,

(9)

~re , i T~.

(10)

M,t

-- - [ N , ( U ) + .{~a ]e + $, -

( ¢ r M - t C ) p . -_ e r A .

(14)

U,.I = Un + At(A n = M-]Cp,,)

(15)

Rather than computing the pressure solution at every timestep, a calculational procedure has been developed whereby the pressure is only occasionally updated based on a preset time truncation accuracy requirement (Gresho, et al., 1984). This method of pressure subeycling is very cost-effective, especially for three-dimensional solutions that are converging toward steady-state conditions. The special implicit treatment of a source term in the eequation (12) will he discussed in the next section.

Numerical Aspects of the Turbulence Model It is well known that solutions of k and £ often exhibit negative values in conjunction with numerical oscillations. Such oscillations typically occur when the grid Reynolds number (UAx/K,) is large (resulting in steep field gradients) or when turbulence levels are low (when k and e values are near zero). Besides being unphysical, even

178

R. L. Lee

very small negative values of k and e can lead to an unstable solution in time since negative turbulence diffusivities may be generated. To avoid numerical instabilities, a "clipping" procedure that readjusts the values of k and e is introduced in the model. More specifically, the negative (or near zero) values are clipped whenever they are found to be less than one-thousandth of the maximum field values at each timestep. Also, unrealistically large Kt's are eliminated whenever a preset maximum value is exceeded by readjusting e. Finally, the time-integration algorithm is further stabilized by treating the e2/k source term in the e-eq-~tion implicitly. This term is made implicit by linearizing and transposing it to the left hand side of the discretized equation. The resulting form of the e-equation is still explicit and can be updated according to the standard explicit algorithm. This numerical procedure was first proposed (in a finite element context) by Betts and Haroutunian (1985) who developed a finiteelement-based k-e model to simulate heavy gas dispersions. Boundary Conditions

In typical problems simulated by this model, the boundary layer flow of interest is "semi-confined" whereby the lower boundary is a solid wall with obstacles protruding from the surface. The remaining lateral and top surfaces can exhibit inflow and/or outflow through their boundaries. Following the normal practice of not attempting to resolve the laminar wall layer, the computational domain is displaced a short distance, h, away from all solid walls. The "slip" tangential velocity along the wall boundaries is computed by applying the following conditions:

c, = 0

au,1

u,1

ox.

IUI

g,~----u., "

Ou,~_ IU--I u,z u.2 /6, Oz,

2

=d

for the velocities and

~h

for k and e ,

THE TRANSPORT AND DIFFUSION MODEL In the two-model approach we use the flow field from the Navier-Stokes flow model to provide advective transport for a dispersion model. The dispersion model called ADPIC (Atmospheric Diffusion Particle-in-Cell) that was chosen for our applications was developed as part of the U. S. Department of Energy's Accidental Release Advisory Capability (ARAC) project. The ADPIC model is a three-dimensional, cartesian, particle-diffusion code for simulating the time-dependent distribution of air pollutants. It solves the three-dimensional advectiondiffusion equation in flux-conservative form. The solution method is based on the Particle-in-Cell technique with the pollutants represented by Lagrangian marker particles that are moved within an Eulerian grid via the pseudo-velocity technique (Lange, 1978). The pseudo-velocity method is based on the following solution approach. Given the conservative form of the nonlinear transport-diffusion equation

~

+V-W~Z-nCVx) -- O,

anak = 0

(16)

where x,, x~l, x.a and U,, U,1, Ua are the coordinates and velocity components normal and tangential to the wall boundary (1< is the von Karman constant and IUI is the magnitude of the velocity). The friction velocity u. is obtained from the log-law velocity profile:

(17)

(lS)

where Z is a scalar concentration, K is the turbulent diffusion coefficient, and UA is the advection velocity derived from the dynamic model, a pseudo-transport velocity Ur can be def'med to be up - u A + v o where

(Cff2k):312 e =

The values of U, k, and e are specified at the inflow plane with zero normal gradient and/or no-penetration conditions applied where appropriate at the lateral boundaries. At outflow boundaries, a finite-element-based normal stress condition that allows the flow to leave the computational domain smoothly is used (Gartling, 1990). This type of outflow condition has been used extensively in finite element calculations for both two-dimensional and three dimensional flows.

u . - -/CVx is a diffusion velocity. X

(19)

Thus, the equivalent equation that has to be solved is given by

~. -~-

V'(xUp) = 0.

(20)

The ADPIC model employs an Eulerian grid of equalsized cells with the concentration, Z, defined at the cell centers and the velocities, UA, Up) and U r defined at the cell comers. It is important to note that the Eulerian grid is used to define the pseudo-velocity field (U r) and to facilitate interpolation of these velocities onto particle positions. The actual transport and diffusion process is

179

Three-dimensionalflow and pollutant dispersion around structures

a description of the terrain within this grid (ADPIC av,ats buildings as terrain), and wind components at each of these ADPIC grid points.

advanced via a Lagrangian algorithm that is gridindependent. Finally, the concentration field can be calculated by summing the masses of the particles within each cell. This particle representation of air pollutants has a major advantage in that point sources can be easily modeled without the problem of source resolution that has plagued Eulerian approaches.

The ADPIC grid is constructed by f'trst specifying a rectangular volume that is usnally a subdomain of the finite element computational domain. This volume is then divided into (n~ x ~ x n,) equal-sized cells with an appropriate subset of these being used to define the building. The cells representing the building are marked as impenetrable for the dispersion calculations.

Because of the differing grid structures in ADPIC, the dispersion model, and the finite element wind field model, we have developed a utility program to function as an interface between the two models. The purpose of this program is to read in output information from the finite element model and to construct the appropriate input data for ADPIC. The input data includes a def'mition of a grid,

Once the ADPIC grid is defined, the utility interpolates the FEM-generated wind fields onto the ADPIC grid. The standard finite element tri-linear basis

4 -

(a) ~

..-~---.~~

~

~

~

~

~

~..

,

,

-

-.

i

•

Y ,

,

,

4 ~

(b) [[!! .

[ ! .

.

.

] ! .

.

.

! .

[ ! ! [ .

.

.

.

.

[ .

[ .

.

.

t

T

!

i

2 Y 0

(c)

4 . . . . . . . . . . . . . . . .

2

- -._-:: : : : : : :..93..

:

:

:

~

:~

:

"''~

:

: : : :

:

!

!

:

:

!

t

[

.001:2 J

,,

Y 0 I

0

,

I

2

,

I

4

,

I

~

6

I

8

,

I

10

,

I

12

x

Figure 1: Steady-state results for turbulent flow over a twodimensional backward-facing step. (a) Velocity vectors; (b) Contours lines of k/U/; (c) Contours lines of d-l,/U/.

i

I.. 14.

180

R.

L.

Lee

functions are used for the interpolation calculation. Each above ground ADPIC grid point, (xp,yp,zv),is identified within the f'mite element grid in terms of the element Ep that contains it and its local coordinates (~31p,~p). Since FEMTKE is employed fkst to calculate the velocities at all finite element nodes within each element, it is then straightforward to calculate the velocity at each ADPIC grid point ff the corresponding local coordinates are known. The local coordinates are obtained by employing Newton's method to invert the tri-linear isoparametric mapping: | XF

=

$

E*I~II* l-'-|

~D

=

= 1(1+~, I~.)(l+vlp1q=)(l+~'p

~=)

and the subscript c¢ denotes evaluations at node (x. With the element number and local coordinates of each ADPIC grid point known, the interpolated velocity field (U.~ is calculated via: | Dr|(X.ii~p,~.l

)

;

~

( ~ l ( ( p l ~ t l ~ p )

lit,

E * l Y l l ~ n ~ d z ,

For consistency, all velocities at grid points that are below terrain are assigned to be zero.

E * I Z l I=|

=

m

(a) 4 t

~ ~ t

t~

tl

i

!

!

!

!

! !

~....~

•

.

.

= i

i

: :

.

~

.

. . .

!

.

!

.

.

•

!

!

!

!

!

:

.

.

I

l

!

!

!

!

:

.

.

.

.

.

.

.

.

. D

,

0

I

,

I

,

I

0

,

I

i

4

I

,

l

:

,

8

I

,

I

,

12

l

I

I

16

20

X

4

(b) -

=.,!.! !

!

!

!

-

~

.

!

!

!

!

!

!

!

•

. ,

.. i

. r

:

:

~

i

!

!

~

~

:

i

i

i

2 . . . . . . . . . . .

0

.

:

::

!

i

:

:

i':

1

I

0

"

,

(22)

I-I

I

I--'[

imT1,,

(21)

Z;

I

,

I

2

,

I

,

I

,

I

4

i

I

6

,

I

,

I

i

8

X

Figure 2: Streamline plots for the region downstream of step. (a) Steamlines; (b) Details of recirculation zone.

I

,

I

10

Three-dimensional flow and pollutant dispersion around structures N U M E R I C A L EXAMPLES

Benchmark Calculation in Two Dimensions In this section we present results for turbulent flow over a two-dimensional backward-facing step to benchmark the performance of the k-e model (FEMTKE). This flow has been extensively simulated over the past two decades and excellent experimental data sets are available for comparison with numerical computations. Although the resulting flow pattern is not overly complicated, it exhibits many features that are present within building wakes such as separation, reattachment, and recirculation. The flow geometry chosen for our calculations is based on the channel flow experiments of Kim (Kim et al., 1980) for a developing channel flow approaching a step. The ratio of the upstream channel width to step-height is 2.0 and the Reynolds number based on the step-height (I-I.) and inlet maximum velocity (Uf) is 4.5 × 104. The computational domain for our calculation extends from three step-heights upstream of the step to 21 step-heights downstream. All domain boundaries are displaced 0.025 step-heights away from the solid boundaries to allow for imposition of the wall models. A fully-developed turbulent profile with values that approximate those reported in Kim's experiments is assumed at the inlet plane. The model was calculated for 17500 timesteps at which time an approximate steady-state condition was reached. Results of the simulation are shown in Figures 1 and 2. Figure 1 depicts the important field variables in the form of velocity vector plot and contour plots of k and e for the portion of the computational domain that is downstream of the step. These are in close agreement with similar results from other numerical models (Betts & Haroutunian, 1985). The steady-state streamline pattern is

shown in Figure 2 with a closeup view of the recirculation zone behind the step. FEMTKE predicted a reattachment length of 7.1 step-heights that is in excellent agreement with the measured value of 7.0 + 0.5 step-heights as reported by Kim (Kim et al., 1980).

Computations of Flow and Dispersion around a Threedimensional Two-Building Complex Numerical simulations of the flow and dispersion around two rectangular buildings are presented here as examples of a practical application of the model. Two buildings, with dimensions (~ × ~ × ~) of 50m × 50m × 40m and 50m × 100m × 20m were staggered 50m apart within a computational domain of 500m × 300m × 150m as shown schematically in Figure 3. The finite element mesh of 49 × 39 × 23 cells was graded such that the finest resolution was located next to the solid boundaries. A steady wind that was calculated from a power-law prof'de given by u = u, (y/y.)'~ where u., y., and ot are 10 m/s, 40m, and 0.1, respectively, approached the buildings from the left. The free-stream Reynolds number based on the 40m building was approximately 2.8×107 and the atmosphere was assumed to be neutrally stable. At the inflow boundary, the profiles of k and e were prescribed by the formulae k = c~"1t2 [ Qm~u/~Y]2 and e = % ld [ ~ ~u/~y]"l where ~. = min [ ~-y, 1¢h,], with h, (the height of the surface layer) taken to be 50m. Slip-wall boundary conditions were used at all solid surfaces (with the grid standoff distance h = 2m), a no-penetration boundary condition was applied at the top boundary, and zero traction conditions were used at the lateral and downwind outflow boundaries. The initial flow field was calculated by "mass-adjusting" (Gresho, et al., 1984) the one-dimensional power-law velocity field within the computational domain.

.J 5

Z

181

.7.

X

Figure 3: Schematic o f computational domain f o r f l o w around a two-building complex.

182

R. L. Lee

The computation of the dispersion from a source requires the FEM velocity field to be interpolated onto an equally-spaced ADPIC mesh with the building represented as terrain. The ADPIC mesh employed here had 45 × 45 x 25 rectangular cells within a somewhat smaller computational domain of 405m x 270m x 125m. The interpolated ADPIC velocity field for the cross-sectional

•

,

.

=~

,

:i

planes at z = 10m and y = 174m is shown in Figure 4. For demonstration purposes, two calculations were performed with point sources placed at (100m, 200m, 40m), source A and (162m, 200m, 4m), source B. The plume associated with each point source was represented by a continuous release of particles throughout the period of integration. As a first approximation, the K's were

~ii

,

•

°

.

.

°

. o

. °

. °

°

° °

° °

° °

°

. °

. °

,

°

.

°

° .

~

. .

.

°

.

. .

. .

.

. .

~

~

. .

. .

~

~

.

~

.

~

i

200

Y

1O0 "-22~Z--2Z2

,

I

.............

,

I

1~

a

I

2~

~

I

300

n

400

X

-

•

v

-

•

•

-

100 Z -

-

.

.

.

.

.

.

.

.

.

.

.

.

50

0

*

I

I

100

200

n

I

*

300 X

Figure 4: Interpolated velocity field on ADPIC mesh. Velocity vector projections are shown for two cross-sectional planes: (a) z = lOm; (b) y = 174m. A and B are locations of point sources.

I

400

j

Three-dimensionalflow and pollutant dispersion around structures

183

(a)

100

~&;£.,.-.

Z

"°o-:-=.o

50

• .." . :., - ".:, :k.-," ,

0

. . . . .

, ....

,

100

~

. "

p . .q",

y

°~.,% • . ,

"

,',~ .

-, "

200

(b)

100

z

•

•

"

. " "....,

•

:

5O

F i

0

Ar

.

.

.

1

.

.

.

.

1

.

."

Z

;

" "

" •

"-

.

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

;.~.,,:zz-.:..,,.,..,

.

; ; . ' " ' i

•

.

I

.

•

l

":

:.

""

"'

~

"

-..."

"

. ' 2 . ' "

"" •

..o

-

.

.

.

.

1

"

"

" ' ' . i

"

"

'

"

'

l

'

'

.'..-

:~

.'

~ ,~. .,.,. ,. . ; ' J-

.

-

.-.

...'

• ..," ."

.,,'"

•

.51""

"

. . . .

.

'

.

"

. . . -

.

.

..I......! • .

-~.

. . . ~. . .

..

..."

:,-

I

'

.

.: : .

:: . ..:. . . . .

,.0.-

v

.

,

-

•

.

" . .,.

:

..

:

,°

•

. o.

•

..

'

"

: : . . . .;. .

..

.

. . /. .~,,

.

.. •

.

.

.

|

.

.

.

100

.

|

.

.

.

.

|

.

.

,

.

.

I

200

.

.

.

I

l

300

i

,

.

I

.

l

.

,

I

.

i

i

400

x

Figure 5:

IZ

'

;i .

•

..

...

100

.

'

......

. . . . ........

•

..

I

Y

.

~-

(c)

o

.. .

•

o •

. •

0

.,.

.

'

.,

..,..?."

"

:

""

X

•

200

,

Particle concentration pattern 4 minutes after particles are released at source A. (a) End view; (b) Side view; (c) Top view.

i

184

R. L. Lee

(a)

10

Z

50 ..

#

I

,

. . . .

*-

- ' ' - -

-'

. . . . . . . .

0

,

.

-,

. . . .

1 O0

•

.

.

.

,

--,

,

.

,

,

,

.

~,

"

i •

-

""

i

""

200

Y

.

,"

,~°, •

.

.

.

,

-

.

,

-

.

,

.

(b)

100 Z

50

0

~;

""

",',

,

I

*

•

.

.

.

.

.

-

•

.

I

. . . .

.

I

.

.

.

.

I

. . . .

l

X

.

i

E-

.

. I

. . . .

.• •

I

.

~

'.

'

.:

%1

"

.

w

,

|

,

,

,

,

I

'

1

(c)

°.

.

.d

200

Y

100

0

,

1 O0

200

X

300

400

Figure 6: Particle concentration pattern 1 minute after particles are released at source B. (a) End view; (b) Side view; (c) Top view•

,

I

Three-dimensionalflow and pollutant dispersion around structures

.

100

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

i,;l

.

~"

:]

.-:."!'I

.|

•

z

5O

•

..

• .

•

0

,,~....

-.

4:,

....

.{.!:~..

• : , ~ - % . T=Tlt~.~d ; . . .

1 O0

185

]

'...-.~" ,.,~1

200 Y

.

.

.

.

t

.

.

.

.

I

.

.

.

.

I

.

.

.

.

I

.

.

.

.

•

.

.

.

.

I

•

.

.

.

.

,

.

.

.

,

(b)

100

o

-'1" •

50

~:~ •

• •

".

" - " ,L

--.-'- :--.,.

0

'

o°

'B

.

•

,

. o o . O

¢..."

, ".:-~' ~ . , ; ~ ! ~ "

o

: •

• **=

o .

-."

"'.'.':'.'

"I'

°

•

~.o

~

"~*o

....

"'-" t :

....

•

"

°

o

t

°*~Q

,.

-'

"'t:",:.

~•

o

• •

• ....

•

...

.

.

I..

"*o

~. .

:

'"

:.

•

",

.oqt •

,

.

.o "

•

"-"

•

"° •

~.

."

.

•

•

~

°°o

"':'."

."

, ....

.o

"

'

"

"

"

'

".'".;

"

" .-"

"

"

.",

" ....

X

'

"

"

"

I'o"

;..~

•

'

;.--"

: o ° •

•

.~.. ."

. .

o ° o . . .

•

•

•

•

~

•

°~,°

.

?s".. :

.*

". ~.:

".'."

; . . . . . . . .

• "

•

•

'

• . .

...

•

(@ o

•

•

o ~o

o

•

o

200

,~

- ' .

-.

.. •

• .',."

.I."

. o.

.::

: - ,...,.')j.; .~

"..J.'..'.:

:

* •

•

o . o °

.9

~, ; . J.. :: :..;,

,

o

*

ow

°l

•

•'

*

o

o

•

F v;

•

•

• "i •

•

o •

•

•

• °

,

°

•

•

•

o

o~.

o

o

•

•o

"°

Y

,tl ~.~

°

".

1 O0

0

I

I

l

I

1 O0

l

I

I

|

I

i

•

.

|

,

•

,

,

I

¢

200

¢

,

•

I

300

l

i

i

I

I

i

i

I

I

I

J

400

x Figure 7: Particle concentration pattern 5 minutes after particles are released at source B. (a) End view; (b) Side view; (c) Top view•

,

~

,

i

186

R. L. Lee

derived from formulas that involve the horizontal (cO and vertical (a,) standard deviations of the plume dimensions with values based on the Pasquill-Gifford curves (Slade, 1968). In the absence of turbulence measurements, the Pasquill-Gifford prescription has been used extensively in many atmospheric diffusion models. The projections of the velocity vectors onto crosssectional planes at z = 40m and y = 174m depict the main features of the flow which consist of separation and recirculation zones above and aft of the dual structures. The flow behind the taller building is channelled toward the left by the shorter, but wider, building downstream. An unsymmetrical, left-tilting, building wake is generated behind the taller building. Note that the recircalation zones behind each building are characterized by relatively low velocities. This suggests that pollutants which are released within these inactive regions will not be appreciably dispersed by advective processes. Figure 5 shows three views of the particle dispersion pattern from point source A at the top of the 40m building. The end (a), side Co), and top (c) views of the particle concentrations at 4 minutes after initiation of the release is shown in this figure. The results suggest that the plume is immediately carried upwards from the source into the mainstream as a result of the deflection of the flow at the upwind edge of the building. As the plume spreads, some of the particles are entrained into the wake regions behind the buildings and are Wansported in a direction against the mean flow. The pollutants from this elevated release reach ground level just slightly downstream of the second building. The lateral extent of the plume is about 100m at distance of 3OOm from the source and the maximum plume height is about 70m. For the second release we selected a contrasting scenario with the point source located near ground level between the two buildings (source B in Figure 4). Figures 6 and 7 show the calculated particle dispersion patterns 1 minute and 5 minutes after the release occurred. At early time the particles tend to aggregate within an area close to the source because velocities within the recirculation cavity between the buildings are relatively low. Most of the particles are transported across the downwind face of the building and are entrained into the separated region around the corner along the adjacent face of the 40m building. Some of the particles diffuse into the mean flow and are carried downstream. The dispersion pattern shown at 5 minutes suggests that the particles have now migrated up from ground level along the downstream face of the building and continued to merge into the separation zone around the corner of the building. The plume now becomes widely dispersed but very high concentration levels are maintained within a region in the neighborhood of the source. The particle dispersion pattern for this release scenario is distinctly different from the former case in which the release occurred on the roof of the building.

CONCLUSIONS We have demonstrated a viable approach of coupling a finite-element flow model FEMTKE to a particle-in-cellbased ADPIC model for simulating the transport and dispersion of pollutants from point sources. The flow model appears to capture many features that are relevant to flow around structures, such as separation and recirculation zones and turbulent wakes. The transport and diffusion model was used to predict dispersion from a plume for two contrasting release scenarios. The particle dispersion patterns resulting from these releases correlate closely with the flow pattern around the buildings. ACKNOWLEDGEMENTS The author wishes to thank Dr. John Leone, Jr. and David Turner of LLNL for developing the interface utility and for assistance in the ADPIC calculations. Thanks also to Dr. Alan Ross of U. C. Davis who performed several flow calculations with FEMTKE. This work was performed under the auspices of the Department of Energy by the Lawrence Livermore National Laboratory under Contract No. W-7405-Eng-48.

REFERENCES Betts, P. L. and Haroutunian, V., K-e Modelling of Turbulent Flow over a Backward Facing Step by a Finite Element Method. Comparison with Finite Volume Solutions and Experiment, in Numerical Methods in Laminar and Turbulent Flows, C. Taylor et al., ed., Pineridge Press, Swansea, 1985, 574-585. Chan, S. T., Rodean H. C. and Ermak, D. L., Numerical Simulations of Atmospheric Releases of Heavy Gases over Variable Terrain, Air Pollution Modeling and its Application II1, Plenum Press, New York, 1984, 295-341. Ferziger, J. H., 1990, Approaches to Turbulent Flow Computation: Applications to Flow Over Obstacles, J. Wind Eng. Ind. Aerodyn., 1990, 35, 1-19. Gartling, D. K., A Test Problem for Outflow Boundary Conditions - Flow Over a Backward Facing Step, Int. J. Num. Meth. Eng., 1990, 7, 953-968. Gresho P. M., Chan S. T., Lee R. L. and Upson C. D., A Modified Finite Element Method for Solving the TimeDependent, Incompressible Navier-Stokes Equation. Part 1: Theory," Int. J. Num. Meth. Eng., 1984, 4, 557-598. Kim, J., Kline, S. J. and Johnston, J. P., Investigation od a Reattaching Turbulent Shear Layer:. Flow Over a Backward-Facing Step," J. Fluids Eng., ASME Trans., 1980, 102, 302-308

Three-dimensional flow and pollutant dispersion around structures Lange, R., ADPIC - A Three-Dimensional Particle-in-Cell Model for the Dispersal of Atmospheric Pollutants and Its Comparison to Regional Tracer Studies, J. Appl. Meteor., 1978, 17, 320-329. l.,ee, R. L. and Leone, Jr., J. M., 1988, A Modified Finite Element Model for Mesoscale Flows Over Complex Terrain," Comp. and Math. with Applic., 1988, 16, 41-56. Lee, R. L. and I.game, Jr., J. M., 1991, Numerical Modeling of Turbulent Dispersion Around Slructures Using a Particle-in-Cell Method, Air and Waste Management Association, 1991, pub. 91-83.13.

187

Murakami S., Computational Wind Engineering, J. Wind Eng. Ind. Aerodyn., 1990, 35, 517-528. Paterson D. A. and Apelt C. J., Simulation of Flow Past a Cube in a Turbulent Boundary Layer, J. Wind Eng. Ind. Aerodyn., 1990, 35, 149-176. Rodi, W, Turbulence Models and Their Applications in Hydraulics, Int. Assoc. Hydraulics Research Publication, Delft, The Netherlands, 1980. Slade D. H., Meteorology and Atomic Energy, TID24190, U. S. Atomic Energy Commission, Oak Ridge, 1968, pp. 65-116.