- Email: [email protected]
Numerical modeling of buoyancy dominated dispersal using a Lagrangian approach
ARlorplknL &mvir-r Printed in G-1 Britain.
Vol. 21. No. 6. pp. 1285-1293.
ooo4-6981/%7 s3.a0+0.00 Pcrymon Journals Ltd.
1987.
NUMERICAL MODELING OF BUOYANCY DOMINATED DISPERSAL USING A LAGRANGIAN APPROACH D. J. GAFFEN*,
C. BENOCCI and D. 0 LIVARI
von Karman Institute for Fluid Dynamics, Chat&e
de Waterloo, 72 B-1640 Rhode Saint G&se,
(First received 28 October 1985 rendreceioed
jar publication 17 November
Belgium
1986)
Abstract-A three-dimensional Lagrangian pollutant dispersal model for passive contaminants has been modified to allow simulation of dispersal of heavy gases and moist, buoyant plumes, e.g. cooling tower plumes. The model uses output from a twodimensional Navicr-Stokes solver for description of the atmospheric boundary layer. For heavy gas dispersion modeling, Lagrangian markers representing the emission are released into the flow field and are displaced by mean winds, turbulent velocities derived from the Langevin equation, and buoyancy forces resulting from density differences between the pollutant and the atmosphere. In the case of moist, buoyant plume dispersal, the Lagrangian markers represent water in both liquid and vapor phasesand carry both dry energy (sensibleheat) and latent heat of condensation/evaporation, and accelerations due to buoyancy are related to temperature differences. Model evaluations were made using experimental tests of passiveand heavy gas dispersion in simulated atmospheric surface layers. Predicted concentration fields are in accord with measurements. Model simulations ofcooling tower plume dispersal in the planetary boundary layer were made using field data from the single natural draft cooling tower at the Philippsburg site in West Germany and the French energy central at Bugey. Predicted plume rise and visible plume outline are in good agreement with observations for various meteorological and source conditions. Key word index: Buoyant dispersion, multiphase transport, Lagrangian model, moist plumes, plume rise, diffusion.
l.INTRODU0ION Reliable predictive models of contaminant
dispersion
in turbulent flows arc required to treat environmental problems ranging from accidental releases of heavy gas to microclimatic changes induced by moist cooling tower plumes. In this context, analytic dispersion models, such as the classical Gaussian model (Hanna et al., 1982), are of limited general value as they do not include the dynamics, thermodynamics, and microphysics of either heavy gas clouds or moist, buoyant plumes. On the other hand, integral methods (e.g. Wigley, 1975) still require important simplifying assumptions (cf. Carhart et al., 1982; Schatxmann and Policastro, 1984 for reviews) which limit their value in strongly perturbed flow fields, such as in the neighborhood of buildings. Differential models, based on numerical solutions of the governing transport and diffusion equations in a Eulerian frame (e.g. Hamza and Golay, 1981; Louis et al.. 1984), require fewer limiting assumptions; however, computer time and storage requirements make these codes more valuable as research tools than as practical predictive models. Better results, in particular the elimination of the problems due to numerical diffusion, can be obtained using a Lagrangian approach in which the behaviour l Present address: National Oceanic & Atmospheric Administration, 6010 Executive Blvd. Rockville, Maryland 20852, U.S.A.
of the pollutant is simulated by markers whose displacements reproduce the statistics of turbulent transport. The average history of a large number of markers then yields an ensemble solution of the dilTusion equation. The fundamental theory of the Lagrangian approach has been developed for passive contaminants, for which case excellent agreement with experimental results has been obtained (e.g. Durbin and Hunt, 1980; Wilson et al., 1981; Baerentsen and Berkowicx, 1984). Extension to the case of buoyant pollutants and to heat as well as mass transfer has. until recently, received little attention and requires introduction of a certain amount of empiricism. In the following sections we present an original modification of the classical Lagrangian model for cases of buoyant and multiphase transport. Comparison with a wide range of experimental data will underline the practical value of such an approach.
2.THEORYOFLACRANClANMODELS
Taylor (1921) was the first to consider turbulent dispersion of Lagrangian particles in the atmosphere and to develop the statistical treatment of the problem. As noted above, Lagrangian modeling involves tracing the trajectories of fluid markers in a turbulent flow field. The problem, then, reduces to finding an appropriate expression for the paths of the particles.
1285
D. J. GAFFEN et al.
1286
Many Lagrangian models calculate particle displacement or velocity using the Langevin equation, a complete derivation of which will not be given here. Interested readers are referred to Durbin (1983) who gives an excellent development of the approach, illuminating the analogies between the Eulerian and Lagrangian formulations. Legg and Raupach (1982) and Gifford (1982) provide concise summaries of the theory. Thomson (1984) discusses the applicability of the Langevin equation to atmospheric dispersion modeling. The simplest formulation of the Lagrangian model is the one proposed by Gifford, where the velocity vector, u, is decomposed into a part correlated through the function R and a random part n: u(r + At) = R(At) o(t) + n(r). Combining expansion
(1) with the following
Taylor
(1) series
c(t + A.t) = o(r) + 2 At + higher order terms
(2)
and dropping higher order terms yields the Langevin equation for velocity: (3) where
(91
u: is the fluctuating velocity component
in the X, direction of the vector velocity u, and +’ the fluctuating component of J/. It should be remarked that, in (6), the mean displacement is made a function only of the turbulent properties at the source, neglecting the effects that nonhomogeneous turbulence and the consequent influence of past history may have on the trajectories at each point within the flow field. Hence, strictly speaking, (6) can be used only when fi = constant. To avoid such a problem, recourse can be made to the formulation proposed by Durbin, where the above features are taken into account. In this context, and in view of the very large scales associated with atmospheric boundary layers, the assumption can be made, even for perturbed flow fields, of a slowly varying medium as suggested by Durbin. Therefore, the trajectories of the markers are computed as sum of random displacements obeying a Gaussian law of variance given by (7) and of a mean displacement, resulting from the non-homogeneity of turbulence, given by: S(l)
=
K !. (jr as p
-
31/B) 2s
] + e-b’)
(e
_
p’- 1 +/Ire-“‘)
(10)
(4)
and ‘l=x.
n(r)
(5)
Integrating the Langevin equation twice yields an integral equation for the particle displacement which can then be averaged to yield an expression for the mean turbulent displacement vector, S(t), as a function of the velocity at the source, v,,: S(t)=(c,/P)(1-e-8’).
(6)
Squaring and averaging the displacement equation provides an expression for the mean square displacement which, when averaged over all values of initial velocity, vc, yields the mean square displacement 2, as a function of the mean square velocity 2 a,Z=(2Ji/p2)(/lf-1+e-8’). As our main interest lies in continuous, sions, Equation (7) can be employed plume variance. As shown, for example, by Durbin Gifford (1982), the velocity scale (7)‘.’ scale l/p can be related to the large scale diffusivity by the relation: &/!I = K where K is the turbulent
u:$’ = K a*/ax,
+K-
1 - R(At)
8=,I
species $I defined as
(7) steady emisto model the (1983) or by and the time Eulerian eddy (8)
diffusivity of the transported
where s is the direction tangent to the local trajectory. K and fi are assumed to be functions of s, to be determined following the approach illustrated in section 3. It has to be remarked that, since Lagrangian particles are supposed to be fluid markers, their properties ought to be, in the strictest case, those of the fluid. Thus the simulation of buoyant particles is not totally consistent with the Lagrangian approach. However, Diehl et al. (1982) introduce a vertical particle velocity to account for plume buoyancy, but the term is simply an application of Briggs’ (1969) mean plume rise formulation. Zannetti and Al-Madani (1983) model buoyant plumes by either (1) tagging each Lagrangian particle with an initial buoyancy, a fraction of which is consumed each timestep, or (2) transforming the Tennessee Valley Authority plume rise formula to yield an expression for vertical particle motion. Cogan (1985) adds a buoyant velocity proportional to temperature difference between the particles and ambient air and modified by a drag coefficient and entrainment factor.
3. THE PRESENT MODEL
The current study (Gaffen et al.. 1985) involves the modification of an existing quasi three-dimensional Lagrangian dispersion model (Benocci et al., 1984; Benocci and Olivari, 1985).
Numerical modcljng of buoyancy dominated dispersal using a Lagrangian approach
The dtfinition quasi thr~-d~mensjonal is used because a threG~imensional dispersion is computed in a two-dimensional flow field, which is predicted by a separate code, namely a finite difference SOiVerfor the incompres~ble~ Reynolds averaged Navier-Stokes equations where turbulence is modeled with the classical eddy diffusivity concept (Foussat, 1981). This approach is a consequence of the limited computer power and time available, which made unrealistic any attempt to solve the Navier-Stokes equations in a three-dirnensjo~al flow field. The relevant equations and assumptions included in the ffow field predictor are briefly resumed in Appendix 1. It should be noted that the equations solved are the full Navier-Stokes onesand not the thin layer approximation. Such a choice makes the solution suitable for the prediction of perturbed two dimensional Aow fields in the atmospheric boundary layer and allows inclusion of obstacles with recirculation zones and jets simulating mamentum flux from a source. The effect of the vertical atmospheric lapse rate upon the turbulence is modeled defining two different momentum diffusivities for the horizontal and vertical directions (respectively qX and q,) and introducing temperature dependent terms in the eddy diffusivity model itself. The Lagrangian dispersion code then uses the mean velocity and eddy diffusivity values as resulting from the NaviGF-stokes solver. The eddy dj~usiviti~ for the transported scalars Ki are obtained from the computed momentum diflusivity values using turbulent Prandtl/Schmidt numbers (see Appendix lf. The iatera1 eddy djff~sivjty K, (not computed by the twodimensional solver) is assumed equal to the streamwise one K,. Apart from this assumption, a drawback of the present approach is that the momentum equations and the diffusion equation are uncoupled, and, therefore, the efFects of density changes and mass forces due to the pollutant are neglected in the calculation of the velocity field and the modeling of turbuhrnce. Such a simplification is acceptable as long as the non-uniformity in density due to the pollutant is much smaller than that due to the atmospherjc lapse rate. in the ~grangian random walk mode!, the dispiacemen1 of a particle during a time step, in each coordinate direction, is obtained adding up the effects of mean flow advection and turbulent transport, This fatter term is itself the sum of two components: 8 mean effect due to the change of turbulent properties in space (to account for non-homogGnejty), common to ail the particles, and a random contribution of spectfied properties. Following thGappr~~hpro~sed by Durbin (1983), the equation for the displacement along the direction xi can be written as: ASi c Vi At + (G/dC)At
+f($)
(11)
1287
where Viis the mean advection velocity along xi*S&f is the mean displacement due to turbulence and obtained applying (10) in the i direction and I($) the random component obeying a Gaussian law, of mean 0 and variance a?, as given by (7). The Lagrangian time scales for the vetocity Ructuation (i/~i~ and the spatial scales ei in the three directions are obtained from the turbulent mass diffusivities Ki using the dimensional relationships: Kj I: ;;Ii/Ji & t a.5
(121 oi
113)
and the velocity scales 3 are obtained from the definition of eddy diffusivity of momentum by apply ing the following relationships:
jj
L=
Kz.5
L
(15)
where ij is the averaged momentum diffusivity [(q; + $4 ?$)/3]*,‘. K, the turbulent kinetic energy and L a turbulent length scale. For an un~rtur~d atmospheric boundary layer (14) reduces to the form: li = t(.iu*r 0, *
(16)
where a, are coelficients dependent on the atmospheric conditions (Plate, 1982) and Y* the shear velocity. The current version of the model can be applied to dispersion of passive and heavy contaminants, dry heated air, or moist (saturated) air. 4. RESUtTS 4.1. Dkpersion of d pffssioe contaminilnl The capability of Lagrangian models to predict dispersion of passive contaminants is already well documented in the literature and the present code has been extensively validated for this class of problems (F&to& and Ohvari, lQg5). As an example we present here a simulation of fully three-dimensional dispersion from a elevated point source in a simulated almospheric boundary layer (Fackreii and Robins, t982). Despite resolution limitations inherent in a threedjmensional mesh, reasonable agreement is obtained for streamwise decrease in maximum concentration (Fig. l), increase in transverse and vrrtical plume length scales (defined as the distance from the centerline where C = 0.5 C,,) op and a, (Fig. 21, and concentration profiles at various downstr~m distances (Fig. 31, normalized with respect to the local C mx* 4.2. Heavy gas dispersion Extension of the dispersion model to buoyancydominated phenomena has been a~mpi~sh~ by considering the buoyant force asa local space property related to local concentration. Particles are assigned a
1288
D. J. GAFTEN er al. 1.0
1000
cuaJ
2 source Q 750
0/Z
source
0.1
500
0.01 0.1
250
0
1
X/Zsource
10
Fig. 1. Variation of maximum concentration of a passive contaminant normalized with respect to upstream velocity, source height, and source strength with distance from elevated point source. -Prediction from code PANACHE. 0 Experiment (Fackrell and Robins, 1982).
f
XQsource
lo
Fig. 2. Evolution of normalixed vertical and transversal plume length scales with normalized downstream distance from elevated point sources. Model prediction and experimental data from Fackrell and Robins (1982). -Prediction from ----code PANACHE. AV Experiment (Fackrcli and Robins,
1982).
described by ASi = ViAt + (dSi(t)/dt)At +f(~2)
- (gi(Pr- ~‘)l~‘)Ar’/2 density w&city ference particle
greater than that of air and are given a vertical component pro~rtioM1 to the density difbetween the particles and air and to the local concentration. Particle displacement is now
3 x/zs,,,,e= 0
Z’zsource 2
2.88
[3 1 0 ~ 0
0.5
1.0
C/Cmox
C/Cmax
3 ZGource oxI2 sourcez3.$g 2 1 0
0 0
0.5 1.0 c/c max
(17)
where pr is the density of the reference state atmosphere without the particles and p’ the density of the atmosphere including the effect of the Lagrangian markers.
i
-:1 0
0.5
1.0
C/cmax
Fig. 3. Con~ntrat;on profifes normalized with respect to maximum concentration downstream of an elevated point source. -Prediction from code PANACHE. 0 Experiment (Fackrell and Robins, 1982).
1289
Numerical modeling of buoyancy dominated dispersal using a Lagrangian approach
which is a function of the heat flux from the source (just as it is normally assigned a mass dependent on the source mass flux). While the particles are treated as water substance, they represent moist air, and the dry energy associated with each particle is that of the heated air. Each particle carries a fixed amount of water substance; however, the fraction of that water which is in liquid state is allowed to vary depending on the saturation of the environment. As liquid-vapor phase change occurs, latent heat is transferred, and the particles carry with them this additional moist energy. At each time step, the thermal energy content and tem~r8ture of each grid cell are computed. Four terms appear in the equation for total thermal energy:
of the present model is that the SOhitiOOS of the flow field and of the dispersion arc independent, and the influence of a heavy gas on the turbulence cannot be taken directly into account. Therefore, it was decided to follow an approach which is relatively common in the modeling of gas dispersal with K models (Ray Phani, 1983)and modify the eddy diffusivity in the vertical direction as: A w&nets
(18)
K = KzU-Clp,)+K,CIp,
where C is the pollutant concentration. K, is the value computed by the mean flow solver in the absence of the pollutant, and K,a heavy gas eddy viscosity. A widely used form for K, is: K,=
(1%
1.26uL31CffU-P&N
3.0
(21)
where C,is the specific heat at constant pressure and C is the con~ntratioo of moist air (m~s/volume) and is related to the concentration of water, CO, via C= Co (Q,,w'Q,,,J
(22)
and Q -and Qnprare the tot81 mass and vapor fluxes from the sources, respectively. (The current version of the model assumes the emission of saturated air by the source.) The factor (p,- C) implies that a fraction C/p, of the reference particles
state mass has been displaced
in the cell, i.e. the atmosphere
by the
is assumed
incompressible.
The dry energy is assigned to each particle at its generation and conserved throughout its trajectory.
4.0
CONCENTRATION(%)
Fig. 4. Concentration profiles downstream of an elevated point source ofa heavy -Prediciion from code PANACHE. 0 Experiment (Botrego and Riethmuller. 1980).
gas.
(20)
‘%ATENT
E,,, = c#,@,--c)
For simulation of buoyant, moist plume dispersal, such as that of cooling tower plumes, we have modified the buoyancy term in the trajectory Equation (17) to allow the density difference to be a function of the energy content of the pollutant. Thus the Lagrangian markers transport not only mass but also thermal energy. Each particle is assigned an initial dry energy,
2.0
EDRY+
The reference energy is simply that of the background, or reference state and d&ends on the reference state temperature 8,:
4.3. Buoyant moist plume dispersion
1.0
EREF+
where each term is an energy per unit volume.
where pp is the particle or heavy gas density (Ermak et ai., 1982). This approach deviates from the strict Lagrangian model in that (1) the poliutant markers are no longer simply fluid markers; and (2) they are not passive in that their motion is affected by the position of other particles via the local concentration. Model evaluation with experimental data from Borrego and Riethmuller (1980) for CBrFJ dispersion from an elevated source in a simulated atmospheric boundary layer. shows good agreement with experimental centerline concentration profiles (Fig. 4).
0
E TOT =
1290
D. J. GAFFEN
The magnitude of this term also depends on the local concentration via Eo,, = &j&/W,
(23)
where &is the dry energy per unit volume per particle and W, the water mass per unit volume per particle. Finally, the latent heat ofevaporation/condensation is computed for each grid cell and each time step by evaluating the total water content of the cell, including both the effect of the particles and that of the background relative humidity. The latent heat of a cell is computed at the beginning of a calculation by associating the existing liquid water in the cell to the energy release required to condense it. Thus
et al.
not always transferred, as this could, in the case 01 condensation, for example, heat the cell to a temperature at which it was no longer saturated. Thus we compute the heat release E,,,required to bring the cell to saturation temperature, where E max= @sat-
‘TOT)P&p
(31)
and take the minimum (in absolute value) of the two estimates from Equations (30) and (31). The model has been validated by simulating several cases for which experimental meteorological, cooling tower, and visible plume data were available, two of which we present below. 4.3.1. Philippsburg Case 18. Data from the 800 MWe natural draft cooling tower at the Philippsburg E’,,,,, = c ( WpfL) (24) power station in West Germany for 1 April 1980 where the summation is made over all particles in the (Policastro et al., 1982; Schatzmann and Policastro, cell, / is the liquid water fraction of the particle (mass 1982) were used as input to the boundary layer and of liquid/total mass of water), and L is the latent heat of Lagrangian dispersion models. condensation. Moderate winds (about 8 m s^ ’ at source height) Then given the total thermal energy of a cell, the and high relative humidity (75-100%) allowed formatotal temperature is computed via tion ofa visible bent-over plume about 800 m long in a neutrally stableatmosphere (Fig. 5). Due to the present l?TOT = ETOT/(prc,) (25) status of the planetary boundary layer model, input velocity and temperature profile were logarithmic and The water content is expressed as a mixing ratio (mass of water/mass of dry air) and is compared to the adiabatic, respectively. saturation mixing ratio associated with the total The modeled visible plume is depicted by all grid temperature of the cell. Actual mixing ratios, qm, are cells in all lateral planes in which condensed water exists at the final timestep. Thus the plume may be dry expressed as the sum of moisture due to the particles, in the core region but contain liquid water in the outer qvlpand q,ig, and that due to non-zero environmental layer and still be visible. relative humidities, q,: Model results for the visible plume are in reasonable (26) agreement with observations, especially considering 4act= q”ap+qllg+qr the non-stationarity of the phenomenon. Modeled where plume rise is slightly less than observed. Application of (27) Briggs’ (1969) plume centerline rise formulae yields an 4vs,+ qi,q = C(wr+ P,) . The saturation mixing ratio q,, is computed by first overprediction of plume rise, indicating the superior performance of the Lagrangian model in this situation. computing the saturation vapor pressure e,, (in mbar) 4.3.2. Bugep Case C. A second moist plume simufrom the Magnus equation: lation was made for the four 900 MWe natural draft log, 0e,, = - 2937.4/O,o,- 4.9283 log, oO,or+ 23.55 cooling towers at the Bugey site in the Rhone Valley in (28) France. Observations at Bugey from 4 March 1980 (Hodin el ol., 1982) show a highrising plume advected and determining the corresponding mixing ratio via only slightly by light (2-3 m s- ‘) winds. Environ(29) mental relative humidities and temperatures were 4saal= evJpr fairly low, about 2060% and 274 K, respectively, where the reference state pressure p,is computed from the local ideal gas law. In undersaturated cells, excess providing good contrast to the Philippsburg case. The four cooling towers were modeled as a single source. available liquid water, as represented by particles with Visible plume results are in excellent accord with non-zero liquid water fractions, is evaporated and heat is absorbed. In supersaturated cells, excess vapor is observation and again are more accurate than the condensed and latent heat is released. At the end of Briggs’ formulation (Fig. 6). Slight overestimation of vertical plume extent is possibly due to inability of the each time step, all particles in a given cell have equal boundary layer code to model the observed temperaliquid water fraction, and any cell containing particles ture inversion at 900 m. with non-zero liquid water fraction is considered part of the visible plume. The latent heat term is then recomputed based on the total amount of water undergoing phase change: 5. DISCUSSION E LATENT=
LPPq,iq.
(30)
In fact the full latent heat associated with the change is
The Lagrangian dispersion model, program PANACHE, is able to simulate dispersion of passive
Numerical modeling of buoyancy dominated dispersal using a Lagrangian approach
1291
800
600
Z(m) 400
01 0
I I
I
300
I 600
I
900
I 1200
1 1500
Xlm)
Fig. 5. Modeled (dashed line with stars)and observed (solid liae)visiblc plume outlines and plume centerline prediction from Brig& formulae (dashed line with squares) for
Philippsburg Case 18.
1500 Z(m1
1200
900
600
-** 300
-
PREDICTION ALGEBRAIC PREDICTION CENTERLINE RISE OBSERVATION
0
0
300
600
900
1200 X(m) 1500
Fig. 6. Modeled (dashed line with stars) and observed (solid line) visible plume outlines and plume centerline prediction (dashed line with squares) from Briggs’ formulae for Bugcy Cast C.
contaminants, heavy gases, and moist, buoyant plumes from cooling towers. Particle displacement is dcten&ted by advcctive velocities and turbulent displacements derived from the Langevin equation. Cases of buoyancy-dominated dispersion are handled by addition of a vertical displacement proportional to
density differences between the contaminated and un~n~minat~ air. For moist cfffuents, phase change is explicitly included. Results for concentration fields and visible plume outlines are in good agreement with cxpcrimcntal data. The limitations of the current study can be attri-
1292
D. J. GAFFEN et 01.
buted mainly to the planetary boundary layer model more sophisticated, three-dimensional model for all stability employed. Linkage of the dispersion mode1 to a
regimes and including explicit calculation of the moisture
field is envisaged and could be expected to lead to
improved
results for the dispersion calculation.
Acknowledgemenrs-This work was carried out at the von Karman Institute for Fluid Dynamics while one author (D. J. Gaffen) was a Fellow of the Belgian American Educational Foundation.
REFERENCES Baerentsen J. H. and Berkowia R. (1984) Monte Carlo simulation of plume dispersion in the convective boundary layer. Atmospheric Environment 18, 701-712. Benocci C. and Olivari D. (1985) Lagrangian modelling of turbulent dispersion. In Pollutant Dispersion, von Karman Institute LS 1985-02. Benocci C., Olivari D. and Vergison E. (1984) Modelling of turbulent dispersion of neutral and buoyant contaminants released from 2D sources.von Karman lnstitute TN 150. Borrego C. and Riethmuller M. L. (1980) Optical measurements of the diffusion of a heavy gas - in a turbulent boundary layer. In Measurements and Prediction of Complex Turbulent Flows. von Karman Institute LS 198003; also von Karman Institute Preprint 1980-09. Briggs Cl. A. (1969) Plume rise. U.S. Atomic Energy Commission. Carhart R. A., Policastro A. J. and Zeimer S. (1982) Evaluation of mathematical models for natural draft cooling tower plume dispersion. Atmospheric Enoironmenr 16, 67-83. Cogan J. L. (1985) Monte Carlo simulations of buoyant dispersion. Atmospheric Enuironmenf 19, 867-878. Diehl S. R., Smith D. T. and Sydor M. (1982) Random walk simulation of gradient transfer processesapplied to dispersion of stack emission from coalfired power plants. J. appl. Mer. 21, 69-83. Durbin P. A. (1983) Stochastic differential equations and turbulent dispersion. NASA RP 1103. Durbin P. A. and Hunt J. C. R. (1980) Dispersion from elevated sources in turbulent boundary layers. J. de Mkanique 19, 679-695 Ermak D. er al. (1982) A comparison of dense gas dispersion model simulations with Burro series LNG test results. J. Haz. Mar. 6, 129. Fackrell J. E. and Robins A. G. (1982) Concentration fluctuations and fluxes in plumes from point sourcesfrom a turbulent boundary layer. J. Fluid Mech. 117, l-26. Foussat A. (1981) Modele de dispersion atmosphtrique nonisotherme d’un polluant gazeux de densite quelconque en presence de non-uniformites orographiques. von Karman Institute TN 135. Galfen D. J.. Benocci C.and Olivari D. (1985)Application ofa Lagrangian dispersion model to environmental problems. von Karman Institute TM 38. Gifford F. A. (1982) Horizontal diffusion in the atmosphere:a Lagrangian dynamical theory. Atmospheric Environment 16,505-512. Hamza R. and Golay M. W. (1981) Behaviour of buoyant moist plumes in turbulent atmospheres. MIT, Energy Laboratory Report EL 81-020. Harma S. R..Br&s G. A. and Hosker R. P. (1982) Handbook on atmospheric diffusion. US Dept. Energy. __- DOE/TIC 11223. Hodin A. er a/. (1982) ModClisation numirique dcs panaches d’atrorkfrigirant. Comparison avec da don&s expCrimentales. ElectricitC de France.
Legg B. J. and Raupach M. R. (1982) Markov chain s~mulation of particle dispersion in inhomogeneous flows: the
mean drift velocity induced by a gradient in Eulerian velocity variance. Boundary-Layer Mer. 24, pp. 3-13. Louis F., Biscay P.and Saab A. (1984) Modtlisation tridimensionnelle de I’impact atmosphtrique d’atrortfrigerants de grande puissanceValidation du code Gedeon. Electricite de France. Plate E. (ed.) (1982) Engineering Mereorolocv. __ Elsevier Publishing. Policastro A. J. er al. (1982) Validation ofcooling tower plume models at Phihppsbura and Gevers sites. Proc. 3rd IA HR Cooling Tower korkshop, Budapest, October 1982. Raj P. K. (1982) Heavy gas dispersion. A state-of-the-art review of the experimental results and models. In Heavy Gas Dispersal, von Karman lnstitute LS 1982-03. Schatzmann M. and Policastro A. J. (1982) A fully Gaussian integral model and its performance with the Philippsburg data. Proc. 3rd IAHR Cooling Tower Workshop, Budapest 1982. Schatzmann M. and Policastro A. J. (1984) Plume rise from stacks with scrubbers: a state-of-the-art review. Bull. Am
Met. Sot. 65, 210-215. Taylor G. I. (1921) DilTusion by continuous movements. Proc. London Mathematical Society 20, 196-211. Thompson D. J. (1984) Random walk modeling of diffusion in inhomogeneous turbulence. Q. JI R. Met. Sot. 110. 1107-1120. Wigley T. M. L. (1975) Condensation in jets, industrial plumes and cooling tower plumes. J. appl. Me?. 14,78-86. Wilson J. D., Thurtell G. W. and Kidd G. E. (198 1) Numerical simulation of particle trajectories in inhomogeneous turbulence-11. Systems with variable turbulent velocity scale. Boundary-Layer Mel. 21, 423-441. Zanetti P. and AI-madani N. (1983) Simulation of transformation, buoyancy and removal processes by Lagrangian particle methods. Proc. 14rh Jnr. Techn. Meeting on Air Pollurion Modeling and its Applicafion, Copenhagen, Denmark, September.
APPENDIX 1 The atmospheric flow predictor solves the Reynoldsaveraged, incompressible Navier-Stokes equations in two dimensions (streamwise and vertical). Turbulent effects are modeled with a non-isotropic eddy diffusivity, to take into account the infhtence of the vertical temperature gradient. The equations for the mean flow are, therefore: Conrinuify
av,_ 0 axi -
(A.])
Momentum
av, av,v,
at+-8x1 Enera v L,
(A.3) where & is the gravity vector (0, 0; -g), p the pressure, Q a source term and q, the non-isotropic total (Irminar and turbulent) di!Tuaivity of momentum, related to the eddy diffusivity of heat/mass, through a turbulent prandtl number Pr,.
Numerical modeling of buoyancy dominated dispersal using a Lagrangian approach momentum itself
For the present application. the relationship K, = rjJ0.87
1293
(A.4)
has been employed for the horizontal direction, and (A.5) for the vertical one. The turbulent Prandtl number is here defined asa function of the Monin-Obukhov length H (Plate, 1982). The turbulence model is an extension of the one proposed by Nee and Kovasznay (Foussat, 1981; Renocci er 01.. 1984). and based upon transport equations for the total diffusivity of
where v is the laminar viscosity, A, B and C empirical constants of the model, La turbulent length scale,and Ri the gradient Richardson number. The above set of equations is solved with a SOLA type algorithm, whose main features are: explicit marching in time; space discretization over a staggered mesh with hybrid finite differences for the advection terms and centered finite difTerencesfor the diffusive ones.
Vol. 21. No. 6. pp. 1285-1293.
ooo4-6981/%7 s3.a0+0.00 Pcrymon Journals Ltd.
1987.
NUMERICAL MODELING OF BUOYANCY DOMINATED DISPERSAL USING A LAGRANGIAN APPROACH D. J. GAFFEN*,
C. BENOCCI and D. 0 LIVARI
von Karman Institute for Fluid Dynamics, Chat&e
de Waterloo, 72 B-1640 Rhode Saint G&se,
(First received 28 October 1985 rendreceioed
jar publication 17 November
Belgium
1986)
Abstract-A three-dimensional Lagrangian pollutant dispersal model for passive contaminants has been modified to allow simulation of dispersal of heavy gases and moist, buoyant plumes, e.g. cooling tower plumes. The model uses output from a twodimensional Navicr-Stokes solver for description of the atmospheric boundary layer. For heavy gas dispersion modeling, Lagrangian markers representing the emission are released into the flow field and are displaced by mean winds, turbulent velocities derived from the Langevin equation, and buoyancy forces resulting from density differences between the pollutant and the atmosphere. In the case of moist, buoyant plume dispersal, the Lagrangian markers represent water in both liquid and vapor phasesand carry both dry energy (sensibleheat) and latent heat of condensation/evaporation, and accelerations due to buoyancy are related to temperature differences. Model evaluations were made using experimental tests of passiveand heavy gas dispersion in simulated atmospheric surface layers. Predicted concentration fields are in accord with measurements. Model simulations ofcooling tower plume dispersal in the planetary boundary layer were made using field data from the single natural draft cooling tower at the Philippsburg site in West Germany and the French energy central at Bugey. Predicted plume rise and visible plume outline are in good agreement with observations for various meteorological and source conditions. Key word index: Buoyant dispersion, multiphase transport, Lagrangian model, moist plumes, plume rise, diffusion.
l.INTRODU0ION Reliable predictive models of contaminant
dispersion
in turbulent flows arc required to treat environmental problems ranging from accidental releases of heavy gas to microclimatic changes induced by moist cooling tower plumes. In this context, analytic dispersion models, such as the classical Gaussian model (Hanna et al., 1982), are of limited general value as they do not include the dynamics, thermodynamics, and microphysics of either heavy gas clouds or moist, buoyant plumes. On the other hand, integral methods (e.g. Wigley, 1975) still require important simplifying assumptions (cf. Carhart et al., 1982; Schatxmann and Policastro, 1984 for reviews) which limit their value in strongly perturbed flow fields, such as in the neighborhood of buildings. Differential models, based on numerical solutions of the governing transport and diffusion equations in a Eulerian frame (e.g. Hamza and Golay, 1981; Louis et al.. 1984), require fewer limiting assumptions; however, computer time and storage requirements make these codes more valuable as research tools than as practical predictive models. Better results, in particular the elimination of the problems due to numerical diffusion, can be obtained using a Lagrangian approach in which the behaviour l Present address: National Oceanic & Atmospheric Administration, 6010 Executive Blvd. Rockville, Maryland 20852, U.S.A.
of the pollutant is simulated by markers whose displacements reproduce the statistics of turbulent transport. The average history of a large number of markers then yields an ensemble solution of the dilTusion equation. The fundamental theory of the Lagrangian approach has been developed for passive contaminants, for which case excellent agreement with experimental results has been obtained (e.g. Durbin and Hunt, 1980; Wilson et al., 1981; Baerentsen and Berkowicx, 1984). Extension to the case of buoyant pollutants and to heat as well as mass transfer has. until recently, received little attention and requires introduction of a certain amount of empiricism. In the following sections we present an original modification of the classical Lagrangian model for cases of buoyant and multiphase transport. Comparison with a wide range of experimental data will underline the practical value of such an approach.
2.THEORYOFLACRANClANMODELS
Taylor (1921) was the first to consider turbulent dispersion of Lagrangian particles in the atmosphere and to develop the statistical treatment of the problem. As noted above, Lagrangian modeling involves tracing the trajectories of fluid markers in a turbulent flow field. The problem, then, reduces to finding an appropriate expression for the paths of the particles.
1285
D. J. GAFFEN et al.
1286
Many Lagrangian models calculate particle displacement or velocity using the Langevin equation, a complete derivation of which will not be given here. Interested readers are referred to Durbin (1983) who gives an excellent development of the approach, illuminating the analogies between the Eulerian and Lagrangian formulations. Legg and Raupach (1982) and Gifford (1982) provide concise summaries of the theory. Thomson (1984) discusses the applicability of the Langevin equation to atmospheric dispersion modeling. The simplest formulation of the Lagrangian model is the one proposed by Gifford, where the velocity vector, u, is decomposed into a part correlated through the function R and a random part n: u(r + At) = R(At) o(t) + n(r). Combining expansion
(1) with the following
Taylor
(1) series
c(t + A.t) = o(r) + 2 At + higher order terms
(2)
and dropping higher order terms yields the Langevin equation for velocity: (3) where
(91
u: is the fluctuating velocity component
in the X, direction of the vector velocity u, and +’ the fluctuating component of J/. It should be remarked that, in (6), the mean displacement is made a function only of the turbulent properties at the source, neglecting the effects that nonhomogeneous turbulence and the consequent influence of past history may have on the trajectories at each point within the flow field. Hence, strictly speaking, (6) can be used only when fi = constant. To avoid such a problem, recourse can be made to the formulation proposed by Durbin, where the above features are taken into account. In this context, and in view of the very large scales associated with atmospheric boundary layers, the assumption can be made, even for perturbed flow fields, of a slowly varying medium as suggested by Durbin. Therefore, the trajectories of the markers are computed as sum of random displacements obeying a Gaussian law of variance given by (7) and of a mean displacement, resulting from the non-homogeneity of turbulence, given by: S(l)
=
K !. (jr as p
-
31/B) 2s
] + e-b’)
(e
_
p’- 1 +/Ire-“‘)
(10)
(4)
and ‘l=x.
n(r)
(5)
Integrating the Langevin equation twice yields an integral equation for the particle displacement which can then be averaged to yield an expression for the mean turbulent displacement vector, S(t), as a function of the velocity at the source, v,,: S(t)=(c,/P)(1-e-8’).
(6)
Squaring and averaging the displacement equation provides an expression for the mean square displacement which, when averaged over all values of initial velocity, vc, yields the mean square displacement 2, as a function of the mean square velocity 2 a,Z=(2Ji/p2)(/lf-1+e-8’). As our main interest lies in continuous, sions, Equation (7) can be employed plume variance. As shown, for example, by Durbin Gifford (1982), the velocity scale (7)‘.’ scale l/p can be related to the large scale diffusivity by the relation: &/!I = K where K is the turbulent
u:$’ = K a*/ax,
+K-
1 - R(At)
8=,I
species $I defined as
(7) steady emisto model the (1983) or by and the time Eulerian eddy (8)
diffusivity of the transported
where s is the direction tangent to the local trajectory. K and fi are assumed to be functions of s, to be determined following the approach illustrated in section 3. It has to be remarked that, since Lagrangian particles are supposed to be fluid markers, their properties ought to be, in the strictest case, those of the fluid. Thus the simulation of buoyant particles is not totally consistent with the Lagrangian approach. However, Diehl et al. (1982) introduce a vertical particle velocity to account for plume buoyancy, but the term is simply an application of Briggs’ (1969) mean plume rise formulation. Zannetti and Al-Madani (1983) model buoyant plumes by either (1) tagging each Lagrangian particle with an initial buoyancy, a fraction of which is consumed each timestep, or (2) transforming the Tennessee Valley Authority plume rise formula to yield an expression for vertical particle motion. Cogan (1985) adds a buoyant velocity proportional to temperature difference between the particles and ambient air and modified by a drag coefficient and entrainment factor.
3. THE PRESENT MODEL
The current study (Gaffen et al.. 1985) involves the modification of an existing quasi three-dimensional Lagrangian dispersion model (Benocci et al., 1984; Benocci and Olivari, 1985).
Numerical modcljng of buoyancy dominated dispersal using a Lagrangian approach
The dtfinition quasi thr~-d~mensjonal is used because a threG~imensional dispersion is computed in a two-dimensional flow field, which is predicted by a separate code, namely a finite difference SOiVerfor the incompres~ble~ Reynolds averaged Navier-Stokes equations where turbulence is modeled with the classical eddy diffusivity concept (Foussat, 1981). This approach is a consequence of the limited computer power and time available, which made unrealistic any attempt to solve the Navier-Stokes equations in a three-dirnensjo~al flow field. The relevant equations and assumptions included in the ffow field predictor are briefly resumed in Appendix 1. It should be noted that the equations solved are the full Navier-Stokes onesand not the thin layer approximation. Such a choice makes the solution suitable for the prediction of perturbed two dimensional Aow fields in the atmospheric boundary layer and allows inclusion of obstacles with recirculation zones and jets simulating mamentum flux from a source. The effect of the vertical atmospheric lapse rate upon the turbulence is modeled defining two different momentum diffusivities for the horizontal and vertical directions (respectively qX and q,) and introducing temperature dependent terms in the eddy diffusivity model itself. The Lagrangian dispersion code then uses the mean velocity and eddy diffusivity values as resulting from the NaviGF-stokes solver. The eddy dj~usiviti~ for the transported scalars Ki are obtained from the computed momentum diflusivity values using turbulent Prandtl/Schmidt numbers (see Appendix lf. The iatera1 eddy djff~sivjty K, (not computed by the twodimensional solver) is assumed equal to the streamwise one K,. Apart from this assumption, a drawback of the present approach is that the momentum equations and the diffusion equation are uncoupled, and, therefore, the efFects of density changes and mass forces due to the pollutant are neglected in the calculation of the velocity field and the modeling of turbuhrnce. Such a simplification is acceptable as long as the non-uniformity in density due to the pollutant is much smaller than that due to the atmospherjc lapse rate. in the ~grangian random walk mode!, the dispiacemen1 of a particle during a time step, in each coordinate direction, is obtained adding up the effects of mean flow advection and turbulent transport, This fatter term is itself the sum of two components: 8 mean effect due to the change of turbulent properties in space (to account for non-homogGnejty), common to ail the particles, and a random contribution of spectfied properties. Following thGappr~~hpro~sed by Durbin (1983), the equation for the displacement along the direction xi can be written as: ASi c Vi At + (G/dC)At
+f($)
(11)
1287
where Viis the mean advection velocity along xi*S&f is the mean displacement due to turbulence and obtained applying (10) in the i direction and I($) the random component obeying a Gaussian law, of mean 0 and variance a?, as given by (7). The Lagrangian time scales for the vetocity Ructuation (i/~i~ and the spatial scales ei in the three directions are obtained from the turbulent mass diffusivities Ki using the dimensional relationships: Kj I: ;;Ii/Ji & t a.5
(121 oi
113)
and the velocity scales 3 are obtained from the definition of eddy diffusivity of momentum by apply ing the following relationships:
jj
L=
Kz.5
L
(15)
where ij is the averaged momentum diffusivity [(q; + $4 ?$)/3]*,‘. K, the turbulent kinetic energy and L a turbulent length scale. For an un~rtur~d atmospheric boundary layer (14) reduces to the form: li = t(.iu*r 0, *
(16)
where a, are coelficients dependent on the atmospheric conditions (Plate, 1982) and Y* the shear velocity. The current version of the model can be applied to dispersion of passive and heavy contaminants, dry heated air, or moist (saturated) air. 4. RESUtTS 4.1. Dkpersion of d pffssioe contaminilnl The capability of Lagrangian models to predict dispersion of passive contaminants is already well documented in the literature and the present code has been extensively validated for this class of problems (F&to& and Ohvari, lQg5). As an example we present here a simulation of fully three-dimensional dispersion from a elevated point source in a simulated almospheric boundary layer (Fackreii and Robins, t982). Despite resolution limitations inherent in a threedjmensional mesh, reasonable agreement is obtained for streamwise decrease in maximum concentration (Fig. l), increase in transverse and vrrtical plume length scales (defined as the distance from the centerline where C = 0.5 C,,) op and a, (Fig. 21, and concentration profiles at various downstr~m distances (Fig. 31, normalized with respect to the local C mx* 4.2. Heavy gas dispersion Extension of the dispersion model to buoyancydominated phenomena has been a~mpi~sh~ by considering the buoyant force asa local space property related to local concentration. Particles are assigned a
1288
D. J. GAFTEN er al. 1.0
1000
cuaJ
2 source Q 750
0/Z
source
0.1
500
0.01 0.1
250
0
1
X/Zsource
10
Fig. 1. Variation of maximum concentration of a passive contaminant normalized with respect to upstream velocity, source height, and source strength with distance from elevated point source. -Prediction from code PANACHE. 0 Experiment (Fackrell and Robins, 1982).
f
XQsource
lo
Fig. 2. Evolution of normalixed vertical and transversal plume length scales with normalized downstream distance from elevated point sources. Model prediction and experimental data from Fackrell and Robins (1982). -Prediction from ----code PANACHE. AV Experiment (Fackrcli and Robins,
1982).
described by ASi = ViAt + (dSi(t)/dt)At +f(~2)
- (gi(Pr- ~‘)l~‘)Ar’/2 density w&city ference particle
greater than that of air and are given a vertical component pro~rtioM1 to the density difbetween the particles and air and to the local concentration. Particle displacement is now
3 x/zs,,,,e= 0
Z’zsource 2
2.88
[3 1 0 ~ 0
0.5
1.0
C/Cmox
C/Cmax
3 ZGource oxI2 sourcez3.$g 2 1 0
0 0
0.5 1.0 c/c max
(17)
where pr is the density of the reference state atmosphere without the particles and p’ the density of the atmosphere including the effect of the Lagrangian markers.
i
-:1 0
0.5
1.0
C/cmax
Fig. 3. Con~ntrat;on profifes normalized with respect to maximum concentration downstream of an elevated point source. -Prediction from code PANACHE. 0 Experiment (Fackrell and Robins, 1982).
1289
Numerical modeling of buoyancy dominated dispersal using a Lagrangian approach
which is a function of the heat flux from the source (just as it is normally assigned a mass dependent on the source mass flux). While the particles are treated as water substance, they represent moist air, and the dry energy associated with each particle is that of the heated air. Each particle carries a fixed amount of water substance; however, the fraction of that water which is in liquid state is allowed to vary depending on the saturation of the environment. As liquid-vapor phase change occurs, latent heat is transferred, and the particles carry with them this additional moist energy. At each time step, the thermal energy content and tem~r8ture of each grid cell are computed. Four terms appear in the equation for total thermal energy:
of the present model is that the SOhitiOOS of the flow field and of the dispersion arc independent, and the influence of a heavy gas on the turbulence cannot be taken directly into account. Therefore, it was decided to follow an approach which is relatively common in the modeling of gas dispersal with K models (Ray Phani, 1983)and modify the eddy diffusivity in the vertical direction as: A w&nets
(18)
K = KzU-Clp,)+K,CIp,
where C is the pollutant concentration. K, is the value computed by the mean flow solver in the absence of the pollutant, and K,a heavy gas eddy viscosity. A widely used form for K, is: K,=
(1%
1.26uL31CffU-P&N
3.0
(21)
where C,is the specific heat at constant pressure and C is the con~ntratioo of moist air (m~s/volume) and is related to the concentration of water, CO, via C= Co (Q,,w'Q,,,J
(22)
and Q -and Qnprare the tot81 mass and vapor fluxes from the sources, respectively. (The current version of the model assumes the emission of saturated air by the source.) The factor (p,- C) implies that a fraction C/p, of the reference particles
state mass has been displaced
in the cell, i.e. the atmosphere
by the
is assumed
incompressible.
The dry energy is assigned to each particle at its generation and conserved throughout its trajectory.
4.0
CONCENTRATION(%)
Fig. 4. Concentration profiles downstream of an elevated point source ofa heavy -Prediciion from code PANACHE. 0 Experiment (Botrego and Riethmuller. 1980).
gas.
(20)
‘%ATENT
E,,, = c#,@,--c)
For simulation of buoyant, moist plume dispersal, such as that of cooling tower plumes, we have modified the buoyancy term in the trajectory Equation (17) to allow the density difference to be a function of the energy content of the pollutant. Thus the Lagrangian markers transport not only mass but also thermal energy. Each particle is assigned an initial dry energy,
2.0
EDRY+
The reference energy is simply that of the background, or reference state and d&ends on the reference state temperature 8,:
4.3. Buoyant moist plume dispersion
1.0
EREF+
where each term is an energy per unit volume.
where pp is the particle or heavy gas density (Ermak et ai., 1982). This approach deviates from the strict Lagrangian model in that (1) the poliutant markers are no longer simply fluid markers; and (2) they are not passive in that their motion is affected by the position of other particles via the local concentration. Model evaluation with experimental data from Borrego and Riethmuller (1980) for CBrFJ dispersion from an elevated source in a simulated atmospheric boundary layer. shows good agreement with experimental centerline concentration profiles (Fig. 4).
0
E TOT =
1290
D. J. GAFFEN
The magnitude of this term also depends on the local concentration via Eo,, = &j&/W,
(23)
where &is the dry energy per unit volume per particle and W, the water mass per unit volume per particle. Finally, the latent heat ofevaporation/condensation is computed for each grid cell and each time step by evaluating the total water content of the cell, including both the effect of the particles and that of the background relative humidity. The latent heat of a cell is computed at the beginning of a calculation by associating the existing liquid water in the cell to the energy release required to condense it. Thus
et al.
not always transferred, as this could, in the case 01 condensation, for example, heat the cell to a temperature at which it was no longer saturated. Thus we compute the heat release E,,,required to bring the cell to saturation temperature, where E max= @sat-
‘TOT)P&p
(31)
and take the minimum (in absolute value) of the two estimates from Equations (30) and (31). The model has been validated by simulating several cases for which experimental meteorological, cooling tower, and visible plume data were available, two of which we present below. 4.3.1. Philippsburg Case 18. Data from the 800 MWe natural draft cooling tower at the Philippsburg E’,,,,, = c ( WpfL) (24) power station in West Germany for 1 April 1980 where the summation is made over all particles in the (Policastro et al., 1982; Schatzmann and Policastro, cell, / is the liquid water fraction of the particle (mass 1982) were used as input to the boundary layer and of liquid/total mass of water), and L is the latent heat of Lagrangian dispersion models. condensation. Moderate winds (about 8 m s^ ’ at source height) Then given the total thermal energy of a cell, the and high relative humidity (75-100%) allowed formatotal temperature is computed via tion ofa visible bent-over plume about 800 m long in a neutrally stableatmosphere (Fig. 5). Due to the present l?TOT = ETOT/(prc,) (25) status of the planetary boundary layer model, input velocity and temperature profile were logarithmic and The water content is expressed as a mixing ratio (mass of water/mass of dry air) and is compared to the adiabatic, respectively. saturation mixing ratio associated with the total The modeled visible plume is depicted by all grid temperature of the cell. Actual mixing ratios, qm, are cells in all lateral planes in which condensed water exists at the final timestep. Thus the plume may be dry expressed as the sum of moisture due to the particles, in the core region but contain liquid water in the outer qvlpand q,ig, and that due to non-zero environmental layer and still be visible. relative humidities, q,: Model results for the visible plume are in reasonable (26) agreement with observations, especially considering 4act= q”ap+qllg+qr the non-stationarity of the phenomenon. Modeled where plume rise is slightly less than observed. Application of (27) Briggs’ (1969) plume centerline rise formulae yields an 4vs,+ qi,q = C(wr+ P,) . The saturation mixing ratio q,, is computed by first overprediction of plume rise, indicating the superior performance of the Lagrangian model in this situation. computing the saturation vapor pressure e,, (in mbar) 4.3.2. Bugep Case C. A second moist plume simufrom the Magnus equation: lation was made for the four 900 MWe natural draft log, 0e,, = - 2937.4/O,o,- 4.9283 log, oO,or+ 23.55 cooling towers at the Bugey site in the Rhone Valley in (28) France. Observations at Bugey from 4 March 1980 (Hodin el ol., 1982) show a highrising plume advected and determining the corresponding mixing ratio via only slightly by light (2-3 m s- ‘) winds. Environ(29) mental relative humidities and temperatures were 4saal= evJpr fairly low, about 2060% and 274 K, respectively, where the reference state pressure p,is computed from the local ideal gas law. In undersaturated cells, excess providing good contrast to the Philippsburg case. The four cooling towers were modeled as a single source. available liquid water, as represented by particles with Visible plume results are in excellent accord with non-zero liquid water fractions, is evaporated and heat is absorbed. In supersaturated cells, excess vapor is observation and again are more accurate than the condensed and latent heat is released. At the end of Briggs’ formulation (Fig. 6). Slight overestimation of vertical plume extent is possibly due to inability of the each time step, all particles in a given cell have equal boundary layer code to model the observed temperaliquid water fraction, and any cell containing particles ture inversion at 900 m. with non-zero liquid water fraction is considered part of the visible plume. The latent heat term is then recomputed based on the total amount of water undergoing phase change: 5. DISCUSSION E LATENT=
LPPq,iq.
(30)
In fact the full latent heat associated with the change is
The Lagrangian dispersion model, program PANACHE, is able to simulate dispersion of passive
Numerical modeling of buoyancy dominated dispersal using a Lagrangian approach
1291
800
600
Z(m) 400
01 0
I I
I
300
I 600
I
900
I 1200
1 1500
Xlm)
Fig. 5. Modeled (dashed line with stars)and observed (solid liae)visiblc plume outlines and plume centerline prediction from Brig& formulae (dashed line with squares) for
Philippsburg Case 18.
1500 Z(m1
1200
900
600
-** 300
-
PREDICTION ALGEBRAIC PREDICTION CENTERLINE RISE OBSERVATION
0
0
300
600
900
1200 X(m) 1500
Fig. 6. Modeled (dashed line with stars) and observed (solid line) visible plume outlines and plume centerline prediction (dashed line with squares) from Briggs’ formulae for Bugcy Cast C.
contaminants, heavy gases, and moist, buoyant plumes from cooling towers. Particle displacement is dcten&ted by advcctive velocities and turbulent displacements derived from the Langevin equation. Cases of buoyancy-dominated dispersion are handled by addition of a vertical displacement proportional to
density differences between the contaminated and un~n~minat~ air. For moist cfffuents, phase change is explicitly included. Results for concentration fields and visible plume outlines are in good agreement with cxpcrimcntal data. The limitations of the current study can be attri-
1292
D. J. GAFFEN et 01.
buted mainly to the planetary boundary layer model more sophisticated, three-dimensional model for all stability employed. Linkage of the dispersion mode1 to a
regimes and including explicit calculation of the moisture
field is envisaged and could be expected to lead to
improved
results for the dispersion calculation.
Acknowledgemenrs-This work was carried out at the von Karman Institute for Fluid Dynamics while one author (D. J. Gaffen) was a Fellow of the Belgian American Educational Foundation.
REFERENCES Baerentsen J. H. and Berkowia R. (1984) Monte Carlo simulation of plume dispersion in the convective boundary layer. Atmospheric Environment 18, 701-712. Benocci C. and Olivari D. (1985) Lagrangian modelling of turbulent dispersion. In Pollutant Dispersion, von Karman Institute LS 1985-02. Benocci C., Olivari D. and Vergison E. (1984) Modelling of turbulent dispersion of neutral and buoyant contaminants released from 2D sources.von Karman lnstitute TN 150. Borrego C. and Riethmuller M. L. (1980) Optical measurements of the diffusion of a heavy gas - in a turbulent boundary layer. In Measurements and Prediction of Complex Turbulent Flows. von Karman Institute LS 198003; also von Karman Institute Preprint 1980-09. Briggs Cl. A. (1969) Plume rise. U.S. Atomic Energy Commission. Carhart R. A., Policastro A. J. and Zeimer S. (1982) Evaluation of mathematical models for natural draft cooling tower plume dispersion. Atmospheric Enoironmenr 16, 67-83. Cogan J. L. (1985) Monte Carlo simulations of buoyant dispersion. Atmospheric Enuironmenf 19, 867-878. Diehl S. R., Smith D. T. and Sydor M. (1982) Random walk simulation of gradient transfer processesapplied to dispersion of stack emission from coalfired power plants. J. appl. Mer. 21, 69-83. Durbin P. A. (1983) Stochastic differential equations and turbulent dispersion. NASA RP 1103. Durbin P. A. and Hunt J. C. R. (1980) Dispersion from elevated sources in turbulent boundary layers. J. de Mkanique 19, 679-695 Ermak D. er al. (1982) A comparison of dense gas dispersion model simulations with Burro series LNG test results. J. Haz. Mar. 6, 129. Fackrell J. E. and Robins A. G. (1982) Concentration fluctuations and fluxes in plumes from point sourcesfrom a turbulent boundary layer. J. Fluid Mech. 117, l-26. Foussat A. (1981) Modele de dispersion atmosphtrique nonisotherme d’un polluant gazeux de densite quelconque en presence de non-uniformites orographiques. von Karman Institute TN 135. Galfen D. J.. Benocci C.and Olivari D. (1985)Application ofa Lagrangian dispersion model to environmental problems. von Karman Institute TM 38. Gifford F. A. (1982) Horizontal diffusion in the atmosphere:a Lagrangian dynamical theory. Atmospheric Environment 16,505-512. Hamza R. and Golay M. W. (1981) Behaviour of buoyant moist plumes in turbulent atmospheres. MIT, Energy Laboratory Report EL 81-020. Harma S. R..Br&s G. A. and Hosker R. P. (1982) Handbook on atmospheric diffusion. US Dept. Energy. __- DOE/TIC 11223. Hodin A. er a/. (1982) ModClisation numirique dcs panaches d’atrorkfrigirant. Comparison avec da don&s expCrimentales. ElectricitC de France.
Legg B. J. and Raupach M. R. (1982) Markov chain s~mulation of particle dispersion in inhomogeneous flows: the
mean drift velocity induced by a gradient in Eulerian velocity variance. Boundary-Layer Mer. 24, pp. 3-13. Louis F., Biscay P.and Saab A. (1984) Modtlisation tridimensionnelle de I’impact atmosphtrique d’atrortfrigerants de grande puissanceValidation du code Gedeon. Electricite de France. Plate E. (ed.) (1982) Engineering Mereorolocv. __ Elsevier Publishing. Policastro A. J. er al. (1982) Validation ofcooling tower plume models at Phihppsbura and Gevers sites. Proc. 3rd IA HR Cooling Tower korkshop, Budapest, October 1982. Raj P. K. (1982) Heavy gas dispersion. A state-of-the-art review of the experimental results and models. In Heavy Gas Dispersal, von Karman lnstitute LS 1982-03. Schatzmann M. and Policastro A. J. (1982) A fully Gaussian integral model and its performance with the Philippsburg data. Proc. 3rd IAHR Cooling Tower Workshop, Budapest 1982. Schatzmann M. and Policastro A. J. (1984) Plume rise from stacks with scrubbers: a state-of-the-art review. Bull. Am
Met. Sot. 65, 210-215. Taylor G. I. (1921) DilTusion by continuous movements. Proc. London Mathematical Society 20, 196-211. Thompson D. J. (1984) Random walk modeling of diffusion in inhomogeneous turbulence. Q. JI R. Met. Sot. 110. 1107-1120. Wigley T. M. L. (1975) Condensation in jets, industrial plumes and cooling tower plumes. J. appl. Me?. 14,78-86. Wilson J. D., Thurtell G. W. and Kidd G. E. (198 1) Numerical simulation of particle trajectories in inhomogeneous turbulence-11. Systems with variable turbulent velocity scale. Boundary-Layer Mel. 21, 423-441. Zanetti P. and AI-madani N. (1983) Simulation of transformation, buoyancy and removal processes by Lagrangian particle methods. Proc. 14rh Jnr. Techn. Meeting on Air Pollurion Modeling and its Applicafion, Copenhagen, Denmark, September.
APPENDIX 1 The atmospheric flow predictor solves the Reynoldsaveraged, incompressible Navier-Stokes equations in two dimensions (streamwise and vertical). Turbulent effects are modeled with a non-isotropic eddy diffusivity, to take into account the infhtence of the vertical temperature gradient. The equations for the mean flow are, therefore: Conrinuify
av,_ 0 axi -
(A.])
Momentum
av, av,v,
at+-8x1 Enera v L,
(A.3) where & is the gravity vector (0, 0; -g), p the pressure, Q a source term and q, the non-isotropic total (Irminar and turbulent) di!Tuaivity of momentum, related to the eddy diffusivity of heat/mass, through a turbulent prandtl number Pr,.
Numerical modeling of buoyancy dominated dispersal using a Lagrangian approach momentum itself
For the present application. the relationship K, = rjJ0.87
1293
(A.4)
has been employed for the horizontal direction, and (A.5) for the vertical one. The turbulent Prandtl number is here defined asa function of the Monin-Obukhov length H (Plate, 1982). The turbulence model is an extension of the one proposed by Nee and Kovasznay (Foussat, 1981; Renocci er 01.. 1984). and based upon transport equations for the total diffusivity of
where v is the laminar viscosity, A, B and C empirical constants of the model, La turbulent length scale,and Ri the gradient Richardson number. The above set of equations is solved with a SOLA type algorithm, whose main features are: explicit marching in time; space discretization over a staggered mesh with hybrid finite differences for the advection terms and centered finite difTerencesfor the diffusive ones.