Numerical simulation of discharged waste heat and contaminants into the South Estuary of the Yangtze River

Pergamon

Mathl. Comput. Modelling Vol. 18, No. 12, pp. 107-123, 1993 Copyright 0 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0895-7177193 $6.00 + 0.00

Numerical Simulation of Discharged Waste Heat and Contaminants into the South Estuary of the Yangtze River L. Yul Shanghai Institute of Electric Power Shanghai, 200090, P.R. China S.-P. ZHU Department of Mathematics, The University of Wollongong Wollongong, NSW 2500, Australia (Received and accepted May 1993) this paper, a two-dimensional depth-averaged mathematical model based on the finite volume approach is formed, which can be used to compute the flow fields and contaminant movements driven by tidal flows in estuaries with variations in the river bed. The paper also presents a set of appropriate outflow boundary conditions for immersed outlets in water, when these outlets are important to the entire flow region. As an example, the distributions and variations of velocity, temperature, and concentration caused by discharging contaminants from seven sources, two of which are submerged, near the south estuary of the Yangtze river have been simulated over a full tidal cycle. The results of the simulation and prediction are presented, and the effects of the outflow boundary conditions are also discussed.

Abstract-In

1. INTRODUCTION It is significant, and very interesting, to compute and predict the pollutant transport and movement caused by the discharge of industrial and domestic effluent into natural water bodies, such as estuaries, where tidal flows play a dominant role in redistributing pollutants and waste heat. Correctly modelling pollutant transport and movement in estuaries is especially important to environmental engineers due to the strong impact of the process of discharging waste through estuaries to the near-shore oceanic environment. In order to correctly simulate pollutant transport, it is vital to model the current field correctly. While most three-dimensional models are used for modelling wind-driven currents in coastal areas [1,2], in the past two decades two-dimensional depth-averaged models have been successfully used in simulating and predicting flows and contaminant transport in natural water bodies with different time and length scales [3-S]. Since waters in estuaries are usually well-mixed and tidal flows are predominantly pressure-driven, the contaminants are distributed fairly uniformly over a water column. Consequently, the two-dimensional depth-averaged models have been shown to give reasonable results for simulating and predicting tidal motions [3,4,9]. They have been widely used to calculate the velocity and transport fields in wide, shallow water bodies for engineering planning, evaluation, and design [lO,ll]. Although the flat-bottom approximation was The authors wish to thank The Centre of Shanghai Atmospheric Environmental Sciences for kindly providing the data used in this paper, both for the field and physical model experiments. The authors would also like to express their gratitude to Bruce Cathers for the valuable comments he made after reading the manuscript. ‘The following address will take effect after the end of 1993: Dr. Liren Yu, Shanghai Institute of Electrical Power, Boom 203, No. 15, Lane 1040 Hejiang Rd., Yangpu Dist., Shanghai 200090, P.R. China.

Typeset by d@-Tf$X KM 18:12-H

1137

108

L. YIJ AND

S.-P.

ZHU

adopted in some earlier depth-averaged modelling exercises for the simplicity of numerical calculations (e.g., in [12,13]), depth-averaged models with both the variations of river beds and water elevations being taken into consideration have become dominant since the 1980’s. More sophisticated two-dimensional models have been developed recently. For instance, some of tidal models can even deal with the “wetting and drying” of tidal flats [14]. Maa [15] developed an efficient two-dimensional model which included the terms of atmospheric pressure gradients, wind stresses, tidal force, bottom and lateral friction, Coriolis force, and non-linear convective acceleration. A finite-element model that operates in the frequency domain was also successfully adopted in modelling tidal circulations [16]. However, extra programming effort and the high run-time costs associated with a finite element model of pollutant transport in natural water bodies, such as bays and estuaries make the implicit finite difference scheme still very attractive in engineering practice. Two-dimensional depth-averaged models are also suitable to model pollutant transport in shallow and broad water regions near pollutant discharging outlets. In these regions, terms that are usually included in large scale simulations, such as atmospheric pressure gradient, tidal force, Coriolis force, evaporation, and sensible heat exchange can be ignored due to the relatively small scale of the problem and the relatively weak influence of these terms. On the other hand, a velocity field with rather smaller scale, such as small scale eddies in a flow field, now needs to be better modelled in order to obtain satisfactory results. Models with primitive variables, such as velocity and water-depth, being solved directly have become quite popular since the 1980’s [17], and have replaced the models with indirect variables, such as the volume flux densities, which were commonly used in the 1970’s [18]. Body-fitted coordinate (BFC) techniques have also been adopted in conjunction with two-dimensional depth-averaged modelling [19]; the influences of river banks and coastlines on the flow fields and pollution concentration fields have been more realistically simulated [20]. Two-dimensional depth-averaged pollutant transport models can still be used when there exist submerged outlets in a domain of interest. However, a successful modelling exercise depends on correctly modelling the submerged outlets, which are mathematically singular points in a computational domain. A widely adopted technique is to add a generalised source term, which consists of all but the transient, convection, and diffusion terms, into the corresponding governing equation. Both Patankar [21] and Yu (221 found that generally satisfactory results were obtained when the additional source term method was used to treat these singular points at the internal nodes. One of the greatest advantages of adopting such generalised source terms is of course to make the extension of the existing algorithms very easy and to render programs currently being developed more versatile. In this paper, a two-dimensional depth-averaged model, with primitive variables being solved directly using the finite volume approach is utilised to simulate the transport-dispersion of waste heat and contaminants in the south estuary of the Yangtze river. The numerical discretization is based on the combination of the finite volume approach and a pressure-velocity correction algorithm. The model has been applied to an engineering project, in which a new electric power plant was planned to be built on the south bank of the estuary of the Yangtze river. The fields of velocity, temperature, concentrations COD (Chemical Oxygen Demand), and BOD (Biochemical Oxygen Demand) in a domain (19.6 km long and 4.2 km wide) near the plant needed to be determined in order to assess whether the industrial waste water and the waste heat discharged from the new plant would be likely to degrade the regional water environment. In this region, there were already six waste heat and pollutant sources (bank-discharging outlets and submergeddischarging outlets) before the new plant was built. Our preliminary numerical results showed that the discharge of industrial waste water from the new plant had negligible effect on the regional water quality, but the discharge from the existing domestic sewage source, Sq, may have certain influences on the warm water bodies discharged from the existing electric power plant Ss and the new plant Sr.

Discharged Waste Heat and Contaminants

2. MODEL

109

FORMULATION

The governing equations for two-dimensional depth-averaged tidal flows in estuaries can be written as: dH x+-

dHii dX

dHe

+dy=

o

7

(1)

(2) 8Hti + dHiz -at ax

aHu2 s---z-gH + dy

(3) where Re is the Reynolds number, bars stand for depth-averaged values through vertical integration; I-~ and 76 are the shear stresses on the water surface and the river bed, respectively; us and us represent the surface velocities in the x and y directions; ceff = v + i& is a depthaveraged effective eddy viscosity coefficient, with i& being depth-averaged turbulent kinematic viscosity and v being the kinematic viscosity of fluids. Correspondingly, the temperature and concentration of pollutants are governed by: 8Hp -& i3Hc

+ aHii?+ + dX

aH=

+

dt+dz where ump, o’ms are the molecular Prandtl number and Schmidt number; FT, Fc are depthaveraged turbulent eddy diffusion coefficients for temperature and concentration, respectively. Equations (l)-(5) are obtained by simply integrating the fundamental three-dimensional Navier-Stokes equations of fluid flows from the bottom to water surface. One should note that the local water depth H cannot be treated as a constant in the process of integration, since the surface elevation & + H is a function of time and position for unsteady tidal flows (see Figure 1). It is because of the variation of the bottom topography and the surface elevation that the last two terms on the right hand sides in equations (2) and (3) are present, which are referred to as the additional source terms.

Figure 1. Illustration of a submerged outlet and a basic control volume.

L.

110

YU

AND

S.-P.

ZHU

Apart from the additional source terms resulting from the vertical integration involved in the depth-averaging process, there are source terms resulting from submerged outlets, such as s,,, %y, ST, and ,!?c of the momentum in the x and y directions, temperature and concentration in the corresponding equations. The source term 3 is linearized as 3 = SC + S, a, where @ stands for any general dependent variable, SC is the constant part of 3, and S, is the coefficient of the first-order term of 3 [21]. All these source terms can be dealt with as if they were part of the boundaries of a computational domain. However, due to their singular behaviour, the boundary conditions associated with these discharging outlets are somehow quite subtle. In our model, the submerged discharging outlets (see Figure 2) were arranged as vertical line sources at the nodal points, with some additional sources in the momentum, temperature, concentration, and pressure-correction equations (the last of which is obtained from the continuity equation as shown in [21]). The discretized forms of these source terms associated with the adoption of control volumes are given in the next section. N

-c

C

c.

, I j wtGL-‘i%,~

W

PY).

Bi,j 'i

I

E

I

S Figure 2. Configuration

of a submerged outlet and nodal grids arrangement

Since the variation of bottom topography and the tidal ranges are generally quite large in the south estuary of the Yangtze river, the variation of the bottom topography and the variation of the surface elevation must be included in the model and be dealt with efficiently. This can be achieved by the introduction of water-depth correction coefficients, which were defined in (231 as:

f&j

W = _h

=

hii + Ahi, ii

’

(6)

where Hij is the local water depth, hij is the local tranquil water depth, Ahij is the variation of hii, and A is a characteristic water depth, which can be chosen to be the averaged tranquil water depth of the main channel of a river reach. These water-depth correction coefficients are generally functions of time. The variation of oiJ provides a consistent means to deal with the variation of water depth and the geometrical changes of boundaries; oij = 0 corresponds to land or islands and non-zero “ii represents a region occupied by water. Although the shallowness of a local water column depends on the horizontal length scale such as tidal wavelength, the definition of water-depth correction coefficient CY%~ in (6) indirectly reflects the shallowness of a local water column, since the horizontal length scale is involved in the nondimensionalized characteristic water depth h. The use of water-depth correction coefficients in conjunction with the finite volume approach will be given in the next section. After the surface velocities in equations (2) and (3) are approximated by the corresponding depth-averaged velocities, ti and C, as shown in [24], and the effective eddy viscosity and diffusion coefficients, i&, I’T, and F, are introduced in the traditional way [25], the system of equations (l)-(5) is closed. However, as the unknown variables, ti, 6, H, and the effective eddy viscosity/diffusion coefficients fiie~, I?*, and I=‘, are all coupled in a highly non-linear fashion, a numerical solution must be sought for these equations.

Discharged Waste Heat and Contaminants

3. NUMERICAL

111

PROCEDURES

With the power-law scheme [21] being used for the convective and diffusion terms, equations (I)-(5) are discretized with the finite volume approach. Although a precise order analysis of the power law for such a highly non-linear system is not available at the moment, the results of a preliminary numerical investigation of a one-dimensional unsteady convection-diffusion equation with a continuous source term [22] indicated that this scheme has lower numerical dissipation and higher numerical accuracy than those of the implicit single step and splitting implicit upwind schemes, in which the diffusion terms are discretized with a second-order central difference scheme and the transient terms are discretized with a first-order implicit difference scheme. All the unknown variables were discretized with spatially staggered grids; the depth-averaged velocities U.and v are computed at the faces of each control volume, whereas the water depth H, the temperature T, and the pollutant concentration 6 are computed at the centre of each cell. The SIMPLEC algorithm [26] was adopted to couple the vectorial quantities U. and V and the scalar quantities H, T, c at each time step. Due to the length of the discretized equations, they have not been included. Because of its unconditional stability, a fully implicit temporal discretization of eqs. (l)-(5) was used; iterations were performed at each time step. The iteration at each time step starts with an initial guess of the pressure (water-depth) field. The velocity components ii and fi corresponding to the guessed pressure field are solved from the momentum equations (2) and (3) with the socalled line-by-line method [21]. Then the pressure (water-depth) values are corrected, with the line-by-line method again, from a pressure correction equation which has been derived in [21] from the continuity equation and velocity correction equations. The pressure-correction equation and the two velocity-correction equations, which were used to calculate the corrected velocities, used in our model are similar to those in [21] and therefore have not been included. Now, the corrected pressure (water depth) is taken as a new guessed pressure (water-depth) and the iteration repeats itself until a pre-specified convergence is achieved. The water-depth correction coefficient technique has been adopted to treat the variations of bottom topography and water surface. With the finite volume approach, the water-depth correction coefficients can be classified, according to the two possible locations of a grid point within a control volume; those associated with the centre point of a control volume, the values of which approximately represent the fraction of volume occupied by fluid in the control volume, and those associated with the centre points of the control faces of a control volume, the values of which approximately represent the void fractions of the corresponding control faces. As mentioned before, care must be taken for the submerged discharging outlets when the discretization is carried out. Since these outlets were treated as part of the boundary for the computational domain, the boundary conditions for a submerged outlet can be derived according to the conservation principle of various dependent variables at the control volume centred at a grid point B [22]. For example, the pressure correction equation (referred as p/-equation later) now contains an additional source term b, which is called mass-source in [21]. For all the cells, where no pollutant outlets are located, this source term is set to zero. On the other hand, if there is an outlet with flow rate Qn in the control volume B, this source term is then set to be equal to Qn. According to the local linearized assumption [21], the additional source term b can be expressed in the linearized form b = S, + S, @, where @ is a general transport variable in the governing equations, and S, and S, are coefficients. In the discretized form of the depth-averaged mathematical model (151, the additional source term appeared as S, Ax Ay H + S, Ax Ay H a’, with Ax, Ay being the x and y direction widths of the control volume, respectively (see Figure 2). In the $-equation, the steady discharge of the vertical linear source at nodal point B should not produce the pressure correction p’. We therefore set Cp= p’ = 0; the additional source term at nodal point B is only Qn = S, Ax Ay H, and the linearized coefficient S, of the constant term

112

L. YU

AND

S.-P.

ZHU

should be equal to &n/Ax Ay H. The derivations of the source terms in the other equations are similar (see [21] for a detailed derivation); they are listed below:

1 QB 2 H Ay(Sz),’

S.

,=_E&! 2 H

l Ay (6~)~’

ii-equation

SPZ,j = -1

g-equation

SC+ =_1% % 2 H Ax (6~)~’

S,. =-Z&Z C%,lfl 2 H

% Ax (by),’

(8)

s. =-Sk 1 Pw 2 H Ax((6y)s’

s..

l Ax (6~)~’

(9)

p’-equation

s,+

T-equation

SCij, z&B*, H AxAy

Plfl,l

pa>,+1

=_s!!? 2 H

= QB --LH AxAy’

(11)

C-equation

(12)

where Qn stands for the discharge strength of the vertical line source at point B, (&r),, (6x), stand for the x-direction distances between B and the east and west control volume faces, respectively, and (by),, (6~)~ are the y-direction distances between B and the north and south control volume faces, respectively (see Figure 2).

4. MODEL 4.1.

General

APPLICATION

Description

In 1990, an electric power plant, with an outlet located at 5’7 in Figure 3, was planned to be built near another existing electric power plant with the existing outlet Ss on the south estuary of the Yangtze river. The outlet ST, from which the cooling water and waste water from an industrial washing process (e.g., water used to wash burnt coal dust) are discharged, was planned to be set up at a location 400m from the river bank. In order to ensure that the outlet is submerged for all the tidal elevations, it was to be placed at a depth of 2.0m below the mean sea level. In the neighbourhood of the new outlet, there were six other discharging sources (four industrial waste water sources indicated by Sr, 5’2, S’s, and Se, one cooling water source indicated by Ss, discharging from the existing electric power plant, and one sewage source indicated by 5’4 in Figure 3). The flow rate as well as other physical parameters for the discharged cooling water are tabulated in Table 1. In order to predict the movement of the warm water bodies and the transport of the pollutant plumes resulting from these discharges, as well as to assess the effect of the main contamination sources on the temperature and quality of river water, a rectangular domain 19.6 km in length (with the outlet Sr located 12 km downstream from the upper reach) and 4.2 km in width was selected as the computational domain. A uniform mesh with 100 x 30 uniform grid points were set up with a resolution of 200 m x 150 m in the along-river and cross-river directions, respectively.

L

*

L

2Lm

SOURCE

POSITIONS

Figure 3. Locations of seven discharging sources.

113

Discharged Waste Heat and Contaminants Table 1. Contamination

Q(m3/s)

sources.

BODbw4

COWwm*)

WC)

Sl

0.111

7.0

109.0

s2

0.012

7.0

200.8

64.5

s3

0.414

7.0

10.0

6.0

s4

6.939 55.0

7.0 16.0

300.5 -

138.5 -

7.0

40.0

25.0

s5 s6

.S7**

0.749 55.0/0.39

16.0

10.6

6.39

* ppm: parts per million. ** S7: cooling water 55.0 m3/s; industrial waste water 0.39 m3/s.

The velocity boundary conditions and the variations of local water depth at the entry section of the upper reach were prescribed according to the measured velocities and water elevations of a full tidal cycle in the dry season of 1990, as shown in Figures 4 and 5. The local tranquil water depth at every grid point was calculated by referring to the topographic map of the river in the region to be modelled.

m

Ttme

( t=hour)

Figure 4. Averaged current velocity section of the upper reach.

Ttme at the entry

( t=hourI

Figure 5. Tidal elevation at the entry section of the upper reach.

Our model was used to determine the velocity and temperature fields and the pollutant movements with a new discharging outlet Sr being added to the existing waste water discharging system. 4.2. The Evaluation of Physical Parameters 4.2.1.

Shear stresses

In the numerical simulation, the surface wind stresses and the surface heat exchange were neglected in this environmental evaluation task. The bottom shear stresses were calculated according to the formula given by Rodi et al. [17]:

rbx

=

I-&$)=

qpii&FTT cosqb Cfp&FTT cos$

’

(13)

’

(14)

114

L. YU

AND

S.-P.

ZHU

where 4 is the inclination angle of the river bank (see Figure 6), Cf is the friction factor, which is proportional to the turbulent eddy viscosity coefficient f&: Cf = C,, + (f - f0)1.5,

(15)

f=$

(16)

in which the subscript zero stands for the values in the centre portion of the river.

Cf was

obtained by utilising the following equation, shown in [17]:

s &I

0

s BO

C&i2dy

=

Crop

0

9

--dy+p cos f#J

soBo (f -

fo)1’5

sdy

= TIPS,

(17)

where P, is a wetted perimeter, Bo is the width of the computational domain, and ~0 is the gross frictional stress of the river reach, which is related to the friction velocity V, as:

U*=G.

(18)

Updated every hour during the computation, U, was taken from the measured data, being within the range 0.029-0.0881 m/s.

Figure 6. The cross-section

4.2.2.

Eddy

viscosity/diffusion

of a typical river bank.

coefficients

For the purpose of an engineering computation, the depth-averaged eddy viscosity/diffusion coefficients were chosen to be linearly related to the water depth, which is a function of time and location, as:

ih = PI u* H, l&kG=P2u*H,

II = P&7

(19) (20)

T he c hoice of these coefficients is still where the local friction velocity is u+ = dm. in dispute in the literature. The correct choice of the proper values of these coefficients has always been one of the most difficult aspects in modelling flows in an estuary. In our modelling exercises, the coefficient ,& in equation (19) was chosen to be 500, which was based on other similar engineering computations reported for the lower reach and the estuary of the same river. With this value, the calculated Gt was within the range lo-100 m2/s, which has a median value 45. Such a median is close to the constant eddy viscosity of 50m2/s adopted in [27]. Furthermore, the constant eddy viscosity adopted in [28] also falls into the above range, which again somehow justified our choice of pi. It was more difficult to determine the coefficients in equation (20), as the field data for temperature and concentration were very scarce. Hence, we selected ps = 0.9, which is suggested in [25], and the corresponding values of I?T, I?c to be between 0.018 and 0.18m2/s, which are very close to the constant value of 0.1 m2/s, suggested in [28].

Discharged Waste Heat and Contaminants

115

4.3. Initial and Boundary Conditions The initial conditions of the velocities, water depth, temperature and concentration were determined by calculating the steady-state solution of equations (1) to (5). Then the time-dependent flow field was computed using a primitive variable method combined with the finite volume approach and the pressure-velocity correction algorithm [21]. Since tidal elevations change rather slowly, a simplified approach was adopted; the water-depth at each grid was set according to the measured tidal elevation at the corresponding time. At each time step, the time step-length At was determined in an inverse way, i.e., At was adjusted until a pre-specified value for the relative mass source for the whole field was satisfied. With the relative mass source being set at 5% for the whole field, the time step-length in the computation turned out to be within the range of 0.25-20 seconds. The residual errors of the two adjacent iterations for temperature and concentration transport equations were set to be less than lo-‘. The boundary conditions for the submerged outlets are prescribed in Section 3 already. Other boundary conditions include those at the river inlet and outlet, those on the open boundary near the river centre, and those on the river banks. At the river inlet, we set u = Vo(t), with Uo(t) being the measured depth-averaged velocity at the river inlet, and assumed that there were no flows in the cross-river direction, i.e., U = 0. The temperature and the concentration were set to T = 7°C and c = lOppm, respectively. At the downstream section of the river, we imposed the zero gradient condition for G, T, and (?, i.e., z

= z

= x

= 0, and assumed again that

there were no flows in the cross-river direction, i.e., V = 0. On the open boundary of the selected computational domain, the velocity in the along-river direction, ii, the temperature, ri’, and the du. dT &? concentration, c, were all assumed to satisfy the symmetry conditions, i.e., = = = 0 8Y dY a?/ and the velocity in the cross-river direction is assumed to be zero, i.e., fi = 0. On the river banks, the non-slip boundary conditions were imposed on the tangential velocities, the non-penetrable conditions were imposed for normal velocities and the concentration, c, and the temperature, T;, was assumed to satisfy the adiabatic conditions. All the above boundary conditions are standard for estuary models in the literature.

5. RESULTS

AND DISCUSSIONS

The numerical simulation was carried out for a full tidal period to provide the necessary information for engineering planning, design and environmental evaluation; the velocity and temperature distributions, the concentrations of COD and BOD were output and are plotted in Figures 7 to 12. In all the modelling exercises so far, we have ignored the heat exchange across the water surface, the interaction of wind shear stress applied on the free surface, and the regression of waste heat and pollutants. The variations in COD and BOD due to the re-oxidation from air across the surface were also assumed to be negligible, as the computational period is rather short, and the computational domain is relatively small. In Figures 7a-7d, the velocity fields, with the new power plant being operated, were plotted at four instants, i.e., at higher high water (HHW), higher low water (HLW), lower low water (LLW), and lower high water (LHW). As can be clearly seen from these figures, the upstream velocities change their directions four times within a tidal cycle, which matches with the boundary conditions being imposed. Near the river bank, there exist regions where eddies were generated; their sizes and strengths were varying with time. Even when the background tidal currents are relatively small, the model seems to perform very well. This is illustrated by the reasonable flow pattern near outlets Ss and 5’7, as shown in Figures 7b and 7d. Since the effect of water depth has been incorporated in the model with water-depth correction coefficients, the computed velocity fields seem to be quite reasonable; especially in the shallow water regions, where the velocities are generally smaller than those in the deep water regions.

L.

116

YU

AND

S.-P.

ZHU

__-~-__--__________~~~~~~~~~~~

_-_-_____-_-_-I___CL

--------_L-_L

--__-__-______--____~~_~~~~~~~

-__-

I_______L__--.~_-__-_~~_

_____

..-

L_LI-_C--_.L-_--_---

-_-_.I_-^LI-WC.---_--_F--___-_-_-__--____^

1=L__

-_-___ _--Pm

_______-_-__ ~____cm.-v_

!_w_________--. ____I_____s__-c

____________-.--~ ___ ________

______

__

____________________-_

__I

_____________________c____________ ~__c____________~-~_~____~~~~~~~-~-~--------

___________-_____-__---

Ti-ee’ee

n .

_m-_-______e __.-_________

L

--cx-~;~~;; A

0.4

n/s

Entry

Veloci!y

HHW; velocity vectors,

t = 7 hours.

(b) HLW; velocity vectors,

t = 9 hours.

(a)

0.04

-

0.7

m/s

m/s

Entry

Entry

Velocity

Velocity (c)

0.02

m/s

Entry

LLW; velocity vectors,

t = 14 hours.

Velocity

(d) LHW; velocity vectors, Figure 7. Velocity

distribution

t = 20 hours.

at different time of a full tidal period.

To demonstrate the effectiveness of the adopted water-depth correction coefficients, a special numerical experiment was conducted in which the velocity field at t = 0 without the utilisation of the water-depth correction coefficients was computed. Figures 8a and 8b were drawn with the same computational conditions. Figure 8a shows the results when the irregular river bed (see Figure 8c) was replaced by a flat bottom. It can be clearly seen, in Figure 8a, that the main stream of the centre portion of the river extends unreasonably from the deep water regions into the

Discharged

shallow water regions; waste heat, and pollutants

water regions the adoption

117

Waste Heat and Contaminants

tend to excessively

move only along the shallow

near the banks of both sides. On the other hand, the main stream of water-depth

correction

the variation of water depth, mainly (see Figure 8b), which is apparently

coefficients

differs from that obtained

in the deep-water more reasonable.

calculated

without

with

considering

region near the centre portion of the river Moreover, in Figure Sa, the variation of

velocity direction is only caused by the boundary shape variation of the river bank; the model with a flat bottom cannot accurately si’mula’ce the variation of velocity directions resulting from the change of local water depth. In Figure Sb, when the water-depth correction coefficients were used, not only the variation of velocity magnitudes but also the variation of velocity directions are influenced by the variation of bottom topography.

(a) Ignoring water-depth

(b) Considering

water-depth

variation;

velocity vectors,

variation;

(c) River sectional Figure 8. Comparison

velocity vectors,

t = 0 hours.

t = 0 hours.

view.

of velocity vectors.

The modelled temperature fields at four different moments are shown in Figures 9a-9d, where the isotherms are drawn from 7°C (background temperature) with 1°C interval. A hot plume is shown to move up (in Figures 9a and 9b) and down (in Figures 9c and 9d) with the main currents of the river. The temperature within the warm water region increases significantly (e.g., 3.6’C at t = 7 hours) and the warm water regions from 5’s and 5’7 being all connected, form a big hot plume along the river bank. This is mainly due to the discharged capacity of the waste heat from both the existing power plant (5’5) and the new plant (ST), that are quite large (both with a capability of I.2 billion watt), and the cooling water discharged from the submerged outlets that is quickly mixed vertically with the environmental water bodies in the neighbourhood of the outlets, and at the same time diffused and dispersed in the horizontal directions. Furthermore, the fact that two outlets were too close to each other and both were located fairly close to the river bank reinforces this effect. The solution to this problem might be to relocate the outlet for the new power plant if the total design discharge cannot be reduced. From our numerical experiments, we also observed that the warm water regions were relatively small during HHW

118

L.

YU

AND

S.-P.

SCALE 900

2900

4000-

l

6900

4900

1

1

1

1

6900

l

I

ZHU

1

inch

18900

I

I

= 2486

I

14900

12900

I

I

m

4

I

I

I 16900

I

I

1-4000

I

-3000 -2000 900

2900

6988

4986

I

0900

10980

(a) HHW; temperature

SCALE 900

2900

4000

I

6988

4900

l

I

I

I

8980

l

I

I

I

I

14900

12900

I

1088

I_

0

16900

t = 7 hours.

1

inch

10900 l

=

,

2486

l

I

m

14900

12900

l

I

I

1 16900

I

I

I4000

l

3000

3000

2000

2000

1000

1000

0 900

2908

4900

69'00

6900

10900

(b) HLW; temperature

SCALE 900 4000c

2900 ,

3000

-

2000

-

,

6900

4900 l

,

0900

I

I

I

I

12900

0

16900

t = 9 hours.

1

ch

in

10900 I

14900

I

= 2486

I

l

I

m

14900

12900 I

4

I 16900

1

1

IJ4688

l

-3000

10000 900

III

1

11

2900

l

6988

4900

’

8900

1

(c) LLW; temperature

SCALE 900 4000

2900

I

I

6900

4900 I

I

I

I

0900 I

lmm

”

l

I 149m

1 h9hCI

t = 14 hours.

1 10900

Ii

l

17900

inch

+

=

12900

ii.11

2486

4

m

14900 I

I 16900

I

114000

3000

3000

2000

2880

1000

1000

0 900

25-88

4900

6900 (d) LHW;

Figure 9. Temperature

8900

18588

temperature

distributions

12900

14900'

16900

0

t = 20 hours.

at different time of a full tidal period.

and LLW periods because of strong convection. During periods of slack water, the discharged waste heat results in the larger warm water regions, as shown in Figures 9b and 9d. In order to test the model-predicted results with the experimental data, a comparison of the modelled temperature distribution and that from a physical model experiment near the moments of HLW and LHW is shown in the Figures 10a and lob, respectively. Unlike the results shown in

Discharged Waste Heat and Contaminants

119

60 (a) Near the moment of HLW

60

70 (b) Near the moment of LHW.

Figure 10. Comparison of surface 3’C temperature rise region (experimental) with depth-averaged 3’C temperature rise region (computational). Shading: Surface 3OC temperature rise region; -: Depth-averaged 3OC temperature rise isotherm.

other figures, the numerical results shown in these two figures were generated after the distance between the outlets 5’s and 5’7 was adjusted to 200m instead of 1000 m. The reason for the adjustment was because 1OOOm between the outlets 5’s and ST was the final version of the design, and 200m distance was used in the experiment, when the physical model was designed prior to the last version of the design blue print being issued. The measured data from the physical model experiment for the surface warm water region with 3°C temperature rise was not complete; the shaded area in Figure 10a did not extend beyond 1200m upstream from Ss. However, even with the limited experimental data set, the simulated depth-averaged warm water region with 3°C temperature rise matches quite well, in general, with that from the experiment. Figure lob shows another comparison near the moment of LHW. The measured surface warm water region in the physical experiment covered about 1.5 km2, which is generally close to the area of depth-averaged region with 3°C temperature rise computed by the model, but the shape of the warm water region has a rather large difference. The computed warm water region does not extend widely into the main channel of the river, and is distributed mainly along the bank.

L. YU

120 2900

908

4900

6988

AND

S.-P.

8900

%W

12900

10900

14900

16900

4000

4000

3000

3000

2000

2000

1000

1000

0 900

2900

4900

6900

0900

(a) HHW; COD concentration 900

2900

4000

4900 11

6988 I

0900

II

I

12900

18900

0

16960

t = 7 hours. 12900

10900 II

14900

I

I

14900 I

I

16980 I

I4000

I

3000 2000 1000

1-0

0 900

2900

4900

6900

8900

113908

(b) HLW; COD concentration 900 4000

2900 I

4900

I

I

2900

6988

I

I

4900

0900

I

I

I

900 4000F

2900 I

I

4900 I

I

6988 I

I

I

I

I

I

16900

I

I

lL900

l

I4000

16900

t = 14 hours. 12900

109e0 I

0

16900

14900

I

12900

10900

090i I

12900

I

(c) LLW; COD concentration

14900

t = 9 hours.

18908

I

9908

6900

12900

I

I

14900 I

I

16900 I

I

l-14000

3000

3000

2000

2000

1000

1000

0 -900

2900

4900

6966

8988

10900

(d) LHW; COD concentration Figure 11. COD concentration

distributions

12900

14900

16900

0 -

t = 20 hours.

at different time of a full tidal period.

This difference may result from the experiment limitations, in which the complete geometrical similarity could not be satisfied due to the fact that the physical model was too small in size. Therefore, the experimental data from the physical model presented deeply stratified flows at this particular moment, whereas it is common knowledge that flows in the estuary of Yangtze River are very much turbulent and, as a consequence, it is probably not reasonable for the stratification being shown in the experimental data. However, owing to the fact that the total volume flow rate of the Yangtze River is quite large even during the winter lower water seasons (from January to March), the temperature rise due to the discharges from both outlets is still lower than the limit for the environmental standard in China. The pollutant concentration distributions are shown in the Figures lla-lld (COD) and Figures 12a-12d (BOD), where the contours are drawn from 1Oppm (background concentration) with a 20ppm interval. From Table 1, it can be seen that the strength of the discharged COD

Discharged Waste Heat and Contaminants 900 4000j-

1868

4900

2900

I

I

6988

I

I

I

I

’

LJJ+Y-h

0900

I

10900

I

I

121 14900

12900

I

I

I

I

16900

I

I

I

I-14088

-

0-

’

900

2980

4900

6900

3900

10900

(a) HHW; BOD concentration 2900

900 4000k

I

I

4900

I

I

h..

I

0900

I

I

10900

I

L

16900

12900

I

I

14900

I

16900

1

1

r-4000

1

-3000

3008 2000

6900

I

14900

12900

1 = 7 hours

‘U.

- 2000

t 1000 0’ 900

’

’

2900

’

4900

6900

B900

10900

(b) HLW; BOD concentration 900

2900

i

4000

I

4900

I

I

6900

I

3900

II

I

12900 t =

I

16900

14900

16900

9 hours.

12900

10900

I

14900

I.1

i

II

14000

i

3000 2000 1000

-

0900

,

2900 I

I

4900 1

900 4000-

2000

-

1000

-

2900

1

I

4900

1

I

-

I

6900

I

I

0900 I

I

t = 14 hours.

10900 I

I

14900

12900 I

I

I

I

16900 I

I

l-4@@@

K o/l0

3000 - 2000

. -

0 900

I

I

2900

I

I

4900

,

6900

8900

18988

(d) LHW; BOD concentration Figure 12. BOD concentration

distributions

1000

0

”

(c) LLW; BOD concentration

3000

*

poo~oa

12900

14488

16900

l

1000 0

t = 20 hours.

at different time of a full tidal period.

that of BOD, which is well-reflected in the computational results; the number of contours in Figures 12a-12d is less than that in Figures lla-lid. Our numerical experiments showed that the discharge of the industrial waste water from the new plant would not significantly change the water quality in the estuary. The numerical experiments also confirmed the conjecture that Sr, S’s, and Sd are the main pollutant sources, as there were three peaks in COD concentration in Figures lla-lld near S I, Ss, and SJ, and two peaks in BOD concentration in Figures 12a-12d near Sz and S*. According to the how rate and concentration of the discharged pollutants at each outlet shown in Table 1, Sd is the main contaminant source. Moreover, due to the fact that outlet S4 is very close to the outlets of the cooling water Ss and ST, the attention should be focused on analysing and evaluating the overlapping region of the pollutant zone and the hot plume. However, because of the strong convection in the estuary of the Yangtze river is greater

than

122

L.

YU

AND

S.-P.

ZHU

tnd highly turbulent flows, the concentrations of COD and BOD near 5’5 and ST are about of 15 ;o 30 ppm, but less than 60 ppm for most of the time. Therefore, although the discharge from Sb vi11 have some influence, the discharged pollutants frorn other sources hardly have any influence In the cooling water discharged from the new electric power plant.

6. CONCLUSIONS In this paper, a two-dimensional depth-averaged model with the variation of water depth being incorporated by means of water-depth correction coefficients is used to simulate and predict the flow, temperature, COD, and BOD concentration fields over a full tidal cycle in the south estuary of the Yangtze River, resulting from the demands of an engineering project. A distinguishing feature of the model is the adoption of the water-depth correction coefficients. Incorporating the water-depth correction coefficients into a depth-averaged model is a very simple and effective computational skill with clear physical background; the simulation and prediction ability and the accuracy of the model can be quite satisfactory in engineering applications. The boundary conditions for the submerged discharging outlets were formed by using the conservation principles of mass, momentum, temperature, and concentration together with the assumption of vertical linear sources. These appear to be very effective in handling boundaries which would otherwise be singular. The discharging boundary conditions used here especially resulted in reasonably good computational results, even when the grid size was rather coarse. The modelled temperature distribution was first compared with that measured from a physical model experiment; reasonable agreement was observed. Then the model was run to simulate the flow field, the temperature distribution, and pollutant concentration, assuming that a new discharging outlet was to be built together with some existing discharging sources in the estuary of the Yangtze River. Through our model simulation, it is found that the discharge of the water used to wash burnt coal dust from the new power plant would not significantly change the present water quality of the estuary. It was found, however, that the discharged warm waste water body would merge with another one into a larger warm water body, which would have a distinct overlapping region with the domestic sewage discharged from 574. Based on this prediction, it has been suggested that further research on the interaction between the discharge from domest,ic sewage and the discharge from the cooling water should be performed to determine possible unfavourable impacts.

REFERENCES 1. G.Z. Forristall, Three-dimensional structure of storm-generated currents, J. Geophys. Res. 79, 2721-2729 (1974). 2. C.J. Hearn and P.E. Holloway, A three-dimensional barotropic model of the response of the Australian northwest shelf to tropical cyclones, J. of Physical Oceanography 20 (l), 60-80 (1990). 3. N.S. Heaps, A two-dimensional numerical sea model, Phil. Truns., Roy. Sot. A265, 93-137 (1969). 4. R.D. Pilgree and D.K. Griffiths, Tidal fronts on the shelf seas around the British Isles, J. Geophys. Res. 83, 4615-4622 (1978). 5. K.P. Holz and G. Nitsche, Tidal wave analysis for estuaries with internal flats, Adv. Water Resow. 5, 142-148 (1982). 6. C.L. Chen and K.K. Lee, Great Lakes river-estuary hydrodynamic finite element model, J. Hydraulic Eng. 117 (ll), 1531-1551 (1991). 7. J.J. Leendertse and E.C. Gritton, A water-quality simulation model for well mixed estuaries and coastal seas, Computation Procedures, Vol. II, R-708-NYC, The Rand Corporation, (1971). 8. J. Kuipers and C.B. Vreugdenhil, Calculations of two-dimensional horizontal flow, Delft Hydraulics Laboratory, Rep. S163, Part I, (1973). 9. A.K. Easton, A tidal model of Corio Bay, Victoria, In NzLmerical Simulation of Fluid Motion, (Edited by J. Noye), 358-369, (1978). 10. D. Fritsch, Ch. Teisson and B. Manoha, Long term simulation of suspended sediment transport application to the Loire Estuary, In Proceedings of the IAHR XX!11 Congress, Technical Session C, pp. 277-285, National Research Council of Canada, Ottawa, Canada, (1989).

Discharged Waste Heat and Contaminants

123

11. S. Onishi, I. Sakai and M. Taga, Study on river behaviour under effects of inertia and Coriolis force, In Proceedings of the IAHR XXIII Congress, Technical Session D, pp. 91-99, National Research Council of Canada, Ottawa, Canada,, (1989). 12. Y.J. Tsai and Y.C. Chang, Two dimensional transient hydrothermal mathematical model, In Proceedings of a lst World Congress on Water Resources, Chicago, Sept. 4, (1973). 13. M. Watanabe, R.F. Harleman and J. Connor, Finite element model for transient two-layer cooling pond behaviour, M.I.T., No. 222, (1975). 14. J.R. Hunter, The user’s manual for two-dimensional numerical hydrodynamic model, Report USO-5, Unit for Coastal and Estuarine Studies, University College of North Wales, (1981). 15. J.Y. Maa, An efficient horizontal two-dimensional hydrodynamic model, Coastal Eng., l-18 (1990). 16. J.J. Westerink, K.D. Stolzenbach and J.J. Conner, A frequency domain finite element model for tidal circulation, Energy Lab., M.I.T., Rep. MIT-EL-85-006, (1985). 17. W.A. Rodi et al., Prediction of flow and pollutant spreading in rivers, In Transport Models for Inland and Coastal Waters, Proceedings of a Symposium on Predictive Ability, (Edited by H.B. Fischer), pp. 63-111, Academic Press, New York, (1981). 18. M.B. Abbott, A. Damsgaard and G.S. Rodenhuis, System 21, Jupiter, A design system for two-dimensional nearly-horizontal flows, J. Hyd., Res. 1 (l), l-28 (1973). 19. Y. Iwssa, Recent development of numerical hydraulics, In Proceedings of 7th Asian and Pacific Regional Division of ZAHR, (November 20-26, 1990), Vol. 4, pp. 15-30, Beijing, China, (1990). 20. P. Yu and H. Liu, Application of boundary fitted coordinate technique on 2-D steady flow of Tail Race River, J. of Hydrodynamics Ser. B 4 (l), 16-23 (1992). 21. S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation and McGraw-Hill Book Company, New York, (1980). 22. L. Yu, The turbulent model and numerical research for the turbulent transport of the contaminants in water environment (in Chinese), Ph.D. Thesis, Hohai University, Nanjing, P.R.C., (1988). 23. L. Yu, Two-dimensional variable water-depth mathematical model for tidal flows and numerical computation for discharges of waste heat and contaminate in south estuary of the Yangtze River, In Proceedings of International Symposium on Environmental Hydraulics, (October 1991), Hong Kong, (1991). 24. L. Yu, A new depth-averaged two-equation (i - 8) turbulent closure model and its application to numerical simulation for a river, J. of Hydrodynamics Ser. B 3 (2), 21-28 (1991). 25. H.B. Fisher, E.J. List, R.C.Y. Koh, J. Imberger and N.H. Brooks, Mtiing in Inland and Coastal Waters, Academic Press, New York, (1978). 26. J.P. Van Doormaal and G.D. Raithby, Enhancement of the SIMPLE method for predicting incompressible fluid flows, Numerical Heat Transfer 7 (2), 147-152 (1984). 27. P. Gu, Experimental Research on Domestic Sewage Discharge at Wai-Guao-Qiao of the Mouth of the Yang&e River (in Chinese), Nanjing Hydraulic Research Institute, Nanjing, (1985). 28. B. Qian and Y. Chen, Two- and Three-Dimensional Numerical Simulation Research on the Effects of Waste Heat Discharge from Ran-Bi Electric Power Plant on Environmental Water Tempemture (in Chinese), Nanjing Hydraulic Research Institute, Nanjing, and Department of Mechanics of Beijing University, Beijing, (1987).