- Email: [email protected]
Modelling traffic air pollution in road tunnels
Pergamon
Atmospheric
Pm: S1352-2310(96)00296-g
MODELLING
TRAFFIC
AIR POLLUTION
ROBERTO
Environment
Vol. 31, No. 10, pp. 1539-1551, 1997 0 1997 Else&r Science Ltd AU rights reserved. Printed in Great Britain 1352-2310/97 $17.00 + 0.00
IN ROAD TUNNELS
BELLASIO
Scientific Consultant, Via Dolomiti 8, 20017 Rho (MI), Italy (First received 20 December 1995 and in final form 13 July 1996. Published March 1997)
Abstract-This paper presents two models for the description of air pollutant concentrations in road tunnels due ts> traffic. Turbulence is assumed to depend on both atmospheric turbulence and vehicle’s motion. Emissions are calculated with the COPERT90 methodology (Eggleston et al., 1991). Up to 34 different vehicle categories can be considered at the same time, with an arbitrary number of vehicles travelling inside the tunnel. Emissions are calculated as a function of the position inside the tunnel and of the time. Three pollutants can be simulated with the current version of the models, CO, NO, and VOC. It is also possible to consider vehicles with null emissions. The models are able to consider the effects of an arbitrary number of sinks. Flow rates and outdoor concentrations are a function of the time for each sink. The equation of conservation of mass has been solved with the control volumes method. Particular attention has been given to the formulation of stability conditions. Sensitivity analysis was conducted to verify the model answer to different input parameters such as initial concentration, boundary concentrations and vehicle-induced turbulence. Examples of application are given for a tunnel with urban traffic regime, including passenger cars with different fuels, light duty trucks, heavy duty trucks and motorbikes, and for an underground railway. Simulations have been carried out for the five working days of the week. 0 1997 Elsevier Science Ltd. All rights reserved. Key word index: Tunnel, finite volumes, piston effect.
INTRODUCTION
The interest about air pollution has grown enormously in the last decades due to a better comprehension of the dangerous effects that some pollutants have on human health (Kagawa, 1984; WHO, 1987). Air pollution in urban areas is due to different sources, but the major responsible source is vehicle traffic. Road-traffic-related pollution should be measured at locations where maximum concentrations are expected (EU directive, 85/203/EEC). These locations can be determined by means of appropriate models capable of describing the atmospheric dispersion in particular geometries where transport processes are not well established. An example of these geometries is the street canyon 1:hat has been extensively studied in the last years botlh from an experimental (Rotach, 1995) and a modelling point of view (Yamartino and Wiegand, 1986; Lanzani and Tamponi, 1995). Other locations where high concentrations can be expected are the tunnels in which there are high emissions of pollutants .in relatively small volumes and possible scarce venlilation. Even if people do not spend very long time,3 in tunnels these high concentrations can be extremely dangerous. For example, carbon monoxide has an attitude about 200 times greater
than oxygen to link with haemoglobin (Kagawa, 1984), forming COHb in a short time. This paper describes the modelling of air pollution due to traffic inside tunnels. The author does not know any other work on the simulation of pollutant dispersion inside tunnels. For this reason, some ad hoc working hypotheses are adopted; they can be verified with validation trials and, if necessary, they can be improved in future works. The development of two Eulerian models is illustrated. The first one is a box model that gives the concentration only as function of time; the second one has been developed with the finite-volume method and is able to give concentrations as a function of time and of position inside the tunnel. In both the models, turbulence is treated considering both atmospheric turbulence and turbulence due to vehicle’s movement. Emissions are calculated with the COPERT90 methodology (Eggleston et al., 1991) as functions of vehicle speed for 34 different vehicle types. The numerical model is able to consider also vehicles with null emissions such as, for example, a train inside the underground tunnel. The next section describes the theoretical formulation which is the same for both the models. Then the numerical formulation for the box model and for the one using the finite-volume
1539
R. BELLASIO
1540
methodology is illustrated. A section is dedicated to the description of emission calculation. Finally, the last two sections present the sensitivity analysis and examples of application.
THEORY
The basic equation solved by the models is the equation of conservation of mass inside the tunnel. It is assumed that the mixing along the directions normal to the tunnel axis (y and z) is effective enough due to the vehicle’s motion so that the pollutants are well mixed. This hypothesis is supported by the work of Staehelin et al. (1995). Hence, the concentration is assumed to be a function of time t and of the distance x inside the tunnel from the left entrance, C = C(x, t). The one-dimensional equation of conservation of mass for a turbulent flow obtained starting from the usual equation of conservation of mass and applying the Reynolds’ decomposition of the variables, is the following:
aed -z+ax=-ax ac
ach
+s-P
(11
where C represents the (mean) concentration (pgm-3), u is the wind speed inside the tunnel (m s- ‘), S = S(x, t) is the source of pollutant (pgmm3 s-i) and P = P(x, t) is the sink term (pgrne3 s-i). The source term is a function of position and time because vehicles move in the tunnel causing emissions in different locations and different times. The source formulation will be discussed later. The model is able to consider two types of sinks, fans and aspirators. Fans catch polluted air inside the tunnel and push it out. Aspirators catch “clear” air outside the tunnel and push it into the tunnel; this last type of sink can simulate openings along the tunnel. The wind speed u represents the mean wind component in the tunnel. It is considered positive when the wind blows from left to right and negative when the wind blows in the opposite direction. The term C’U’ represents the turbulent flux of mass inside the tunnel. Often in the lower layer of the atmosphere this term is satisfactorily described by means of the gradient theory. Turbulence inside the tunnels is not only due to the atmospheric turbulence but also due to the vehicle’s motion; in this work an ad hoc formulation of the turbulent term is adopted. It is assumed that the turbulent mass flux is the sum of two terms, one related to the atmospheric turbulence and the other to the vehicle’s motion: -7-t
C u = T~tm+ TV,,,
where T indicates a turbulent term. The atmospheric turbulent flux T,,, is described by means of the gradient theory (Zannetti, 1990). It is assumed to be proportional to the concentration gradient via the turbulent dispersion coefficient K,.
The turbulent dispersion coefficient can be calculated as a function of the friction velocity u* and of the Monin-Obukhov length I. (Zannetti, 1990). The horizontal dispersion coefficient K, is proportional to the vertical dispersion coefficient K,, K, = /lKL. Since no radiation hits the ground at the centre of the tunnel, neutral stability is assumed, the resulting K,, averaged over the tunnel cross-section, is K = k/W& X 2
(3)
where H, is the tunnel height (m), k is the von Karman constant (k = 0.4) and u* is the friction velocity (m s- ‘). The turbulent dispersion coefficients outside the tunnel are calculated with the similarity theory. In the rest of the tunnel the coefficients are calculated interpolating from the external values to the internal one given by equation (3). The turbulent flux due to the vehicle’s motion, TVeh, is described assuming that each vehicle determines, with its movement, the movement of a quantity of air with a speed that is proportional to the ratio between the vehicle cross-section and the tunnel cross-section (piston effect). The turbulent flux of mass at a point will be the product between this piston velocity and the concentration at that point. Hence, the turbulent transport of mass is described with C’u’= - K, 2
+ aC F
ceil
uVeh
(4)
where CIis an empirical constant, A,,,, is the crosssection of a particular vehicle type (m’), Ace,, is the tunnel cross-section (m2) and uVehis the vehicle speed (ms-‘). Substituting equation (4) into equation (l), we have
ac x
acu
+dx
a =-z
qCuve,-K~~
>
+S-P
(5) where rl is the following non-dimensional &$ cell
quantity: (6)
The vehicle-induced turbulence is then treated as an advection term; in other words, the vehicle’s motion contributes to increase or decrease the wind speed inside the tunnel. Equation (5) can be solved provided the initial concentration and boundary conditions are specified. The external concentration near the tunnel entrances is influenced by the high value of concentration inside the tunnel. For this reason, the boundary concentration is assumed to be constant only farther from the tunnel. Near the tunnel the boundary concentration is calculated as a function of the time solving a three-dimensional conservation equation in a box. Equation (5) has been solved in two different ways, with a box-model formulation assuming that the pollutant is well mixed inside the whole tunnel (C = C(t))
Modelling traffic air pollution and with a one-dimensional formulation (C = C(x, t)).
THE BOX-MODEL
finite-volume numerical
concentration C exits from it. It must be noted that when C,,, is greater than C then P* is a source term. The general numerical solution of equation (8) is Max( -w,O)
FORMULATION
The relatively simple formulation of the box-model helps to understand the physical phenomena responsible for the temporal variation of concentration inside the tunnel. Let us consider a tunnel characterised by a length L (m) and a cross-section A (m’), the concentration at the left boundary of the tunnel is called CBL@gmm3) and the one at the right boundary is called CBR(pg m - ‘). The integration of equation (5) over the tunnel volume gives the following equation stating the integral mass conservation:
AL; =A
[
C(u + rpV,h) -
-A w
Kg
1 Left
1
ac + V%eh) - Kxz Righ,
[
+ s* - p*.
(7)
Equation (7) states that the temporal variation of mass is due to the difference between the mass entering or exiting in the tunnel from the left and the mass exiting or entering in the tunnel from the right, plus the difference between the mass “produced” inside the tunnel and the mass exiting from the sinks. The source term S* and the sink term P* have the dimensions of mass per unit of time (pgs-‘); they are, respectively, equal to the source term S and the sink term P times the tunnel volume AL. The numerical formulation of advection and dispersion in equation. (7) is the following: A,!,:
1541
= A[CBL lMax(w, 0) - C Max( - w, 0)]
+ Max(w, 0) + T )) +-
+
CaLAt
L C&t L
+ -$(S*
- R*(C” - co,,))
(11)
where At is the discrete temporal step, C” is the concentration at time t” = n At and C”+ ’is the concentration at time t”+ ’ = (n + 1)At. Equation (8) can be solved iteratively provided that an initial condition Ci is furnished. In order to achieve numerical stability, equation (8) is solved with a temporal step At such that
L
At <
(12) Max( - w, 0) + Max(w, 0) + %
+ T
It is obtained forcing the coefficient of C” in equation (11) to be positive. It can be seen that this condition is not very critical; as an example, for a tunnel 1 km long, a speed w of 5 m s- ‘, a turbulent atmospheric coefficient of 2 m2 s-l, a sink rate R* of 10 m3 SK’and a tunnel cross-section A of 50 m2 it is obtained a At of about 200 s.
THE ANALYTICAL
SOLUTION
OF THE BOX-MODEL
FORMULATION
- A [C Max(w, 0) - CBRMax( - w, 0)] Under special conditions it is possible to give an analytical solution of equation (8). Some of these L conditions are not very realistic but they are useful +S* --P* (8) because the comparison between the analytical solution and the numerical solution gives an idea of the where the speed w (rn s- ‘) is given by goodness of the numerical scheme adopted. If w > 0 and CaL = CaR = C,,, = Ca then equation w = u + fluve,,. (9) (8) assumes the form The temporal derivative is discretised by means of the finite-difference method. The sink term P* can be AL$f=ACB W+%i-T described by the following relation: ( > -
C-
AK,rp
CBL
Cm-C +AK,----L
Pl:= R*(C- C,",)
(10)
where R* (m3 s-‘) is the rate at which air is exchanged between the tunnel and the external environment. C,,, (pgme3) is the external pollutant concentration near the sinks. Equation (10) indicates that when a volume of air with external concentration Goutenters into the tunnel, the same volume of air with internal
+S* >
(13)
whose solution, assuming that w, K,,S*,R* are constants and Ci = 0 is c(t)=C(m)(l-exp(
-4))
(14)
R. BELLASIO
1542
where the asymptotic concentration
C(c0) is given by
s*
C(m) = CB + Aw+~
2AK,
(15) $R*
and where the characteristic time r is given by L
Z= w+y
2K,
R*’ +A
(16)
It can be seen that when t $ r, equation (14) tends to the value C(a) , after t x 57 the concentration C(t) is greater than the 99.3% of the final concentration C(a). The characteristic time r represents the time after which the concentration is equal to the 63.2% of C(a). The characteristic time r increases with the tunnel length L and the tunnel cross-section A, while it decreases with the advection speed w, the dispersion coefficient K, and the sink rate R*. The final concentration C(a) increases with the source term S, the external concentration Ca and the tunnel length L, while it decreases with the sink rate R*, the advection speed w, the atmospheric dispersion coefficient K, and the tunnel cross-section A.
speed due to the vehicle’s motion, and is described by equation (9). It can be seen that the advection term is described by the up-wind or down-wind formulation as a function of the wind direction, and the dispersion term is described with the finite-difference method. Applying the finite-difference method to the temporal derivative of the concentration in equation (17), the concentration at the time step (n + l)At can be written in terms of the concentration at time step n At. Hence, starting from an initial condition Ci it is possible to calculate the concentrations at each time step. The stability condition for equation (17) has been found to be Ax
At< Kx, + KX+, + y Ax
+ Max(w, 0) + Max( - w, 0) (19)
The value of At is calculated using the values that maximise the denominator of equation (19). Another stability condition arises from the three-dimensional equation of conservation of mass solved outside the tunnel to obtain the external concentration near its entrances. This condition is practically the sum of three terms similar to the one in equation (19), one for each direction.
THE FINITE-VOLUME FORMULATION
Equation (5) is solved with the finite-volume method (Patankar, 1980) which guarantees a perfect mass conservation. The tunnel is divided into grid boxes along the x direction and the equation is integrated over each of these boxes. The following integral equation of conservation of mass is obtained:
THE EMISSION SUBMODEL
Emission factors for traffic were calculated as functions of vehicle speed with the COPERT90 methodology (Eggleston et al., 1991). Only hot engines emissions were considered. Vehicles were grouped in categories on the basis of fuel type, the age of the vehicle, the engine cylinder capacity and the weight. AAX: = A(@i-Qi+l) + AAXSi A total of 34 vehicle categories has been considered, each identified by a number as illustrated in Table 1. - A AX Ri(Ci - CpU’) (17) The term PC indicates passenger cars, LDT indicates light duty tracks, HDT indicates heavy duty tracks where Ax is the grid size (m), @i(pgm-‘s- ‘) represand 2W indicates two wheelers. Subscripts indicate ents the mass flux at the position (i - 1)Ax and the fuel type, g stands for gasoline, d stands for diesel cDi+i(~grn-~ s- ‘) represents the mass flux at the posiand lpg stands for liquid petroleum gas. Passenger tion iAx. These fluxes describe the turbulent dispersion and the advection of mass. As an example the flux cars are also split according to the different steps of international legal conformity (ECE classes). Buses Qi at position “i” is presented: are described as heavy-duty trucks. @i = Cj_iMax(w,O) - CjMax(-w,O) Vehicle’s fluxes in the tunnel must be specified at each hour of simulation. For each vehicle type _K cj-cj-l W-9 the average number of vehicles per unit time cp x, Ax ’ (vehicles h-i) and the average speed must be specified. Index “i” represents the positions of the bounding Vehicles entering from left and from right are treated surface of the grid box while index j represents the independently. The temporal distance between two calculation position; each calculation position is in vehicles of the same type in the same hour remains the middle of a grid box. Then, for example, fluxes are constant. It is calculated as r, = 3600/p and is exevaluated at the grid box faces and concentrations are pressed in seconds. Vehicles are treated as mobile evaluated at the middle of each grid box. The velocity sources, this means that, considering for simplicity term w is the sum of the component of the atmoonly one vehicle in the tunnel, the emissions are alspheric wind speed inside the tunnel and the turbulent ways zero but not in the positions occupied by the
Modelling traffic air pollution
1543
partial differential equations, a comparison between the numerical solution of the box model and the analytical solution has been carried out. If this comCylinder parison gives good results, it is reasonable to think ID Vehicle Category capacity that even the numerical scheme adopted for the soluPRE ECE cc< 1.4 1 tion with the finite-volume methodology is correct. In peg 2 PRE ECE 1.42.0 3 pa other. 4 ECE 15-00 cc< 1.4 F-3 Only vehicles ECElS-04 with cylinder capacity ECE 15-00 1.42.0 peg ECE 15-01 cc< 1.4 I peg were examined. A tunnel 1 km long, 8 m wide and 5 m 8 ECE 15-01 1.42.0 9 peg with a wind speed of 2 m s- ‘, a turbulent dispersion 10 ECE 15-02 cc < 1.4 PC% coefficient of 2.3 mz s-l, a boundary concentration of ECE 15-02 1.42.0 15 peg observe their effect over the pollutant concentration ECE 15-04 16 cc< 1.4 PC&? ECE 15-04 1.42.0 18 pa lent dispersion coefficient were changed together beImpr. Conv. cc< 1.4 19 pa cause the dispersion coefficient depends on the Impr. Conv. 1.42.0 26 practically equal. The asymptotic concentration was All 21 PClpg LDTg All 28 determined as the final concentration after a 6 h simuLDTd All 29 lation. The characteristic time was determined as the _ HDTd 3.5t < wght < 16t 30 time at which the concentration over the asymptotic _ HDTd wght<16t 31 concentration is the nearest to 63.2%. The maximum _ 2w cc<50 32 percent difference between the analytical concentra2w 2 stroke cc>50 33 2w 4 stroke 34 cc>50 tion and the numerical one is about 1.5%; the remaining are few decimals or hundreds of percent. The maximum percent difference between the charactervehicle during the temporal step considered. Foristic times is about 1.8. It is useful to note that this characteristic time could not be calculated exactly mulated mathematically, if P1 is the position of the vehicle at time tl, the vehicle’s position P2 at time because the calculation time step used was 10 s. t2 will be According to equation (16) a variation of the number of vehicles per hour or a variation of their speed Pz =:Pi + UV& - tJ (20) (i.e. a variation of S*), does not affect the characteristic where uVehis the vehicle’s speed. During the temporal time z. The same is true of variation in the tunnel cross-section A because in this particular case the sink step At = t2 - tl emissions will be different from zero rate is null. The modification of the concentration and only at the points between PI and P2.For this reason an emission submodel has been constructed in order of the characteristic time as a function of the input variables as calculated by the numerical model agrees to calculate the position of all the vehicles inside the with the one given by the analytical solution. tunnel. The number of vehicles inside the tunnel increases every time a new vehicle enters and decreases every time a vehicle exits. For the box model formulation the emissions are SENSITIVITY ANALYSIS calculated only as a function of the number of vehicles inside the tunnel, since the concentration is not a funcSensitivity analysis is a fundamental part of the tion of the position inside the tunnel. evaluation of a model because it identifies its critical inputs and allows the evaluator to determine if nature displays the same sensitivity to these inputs as the model (Ermak and Merry, 1988). COMPARISON BETWEEN THE ANALYTICAL In order to carry out the sensitivity analysis of the AND THE NUMERICAL SOLUTION finite-volume model, a scenario has been created with In order to obtain an evaluation of the goodness of a tunnel 1 km long, a width of 8 m and a height of 5 m. the numerical scheme adopted for the solution of the The tunnel length has been divided in 50 grids 20 m
Table 1. Vehicle categories for which emission factors are calculated
1544
R. BELLASIO Table 2. Comparison between the analytical solution and the numerical solution of the box model Input
w (ms-‘)
K, (m’s_‘)
R* (m3s-r)
2.3 1.2 3.4 4.5 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3
0 0 0 0 2 4 6 8 0 0 0 0 0 0 0 0 0
2 1 3 4 2 2 2 2 2 2 2 2 2 2 2 2 2
1600
Analytical
A (m2) F (vehh-‘) 40 40 40 40 40 40 40 40 50 60 70 40 40 40 40 40 40
Uveh(kmh-i)
100 100 100 100 100 100 100 100 100 100 100 100 100 100 50 150 200
100 100 100 100 100 100 100 100 100 100 100 50 75 125 100 100 100
Numerical
r (s)
C, (pgme3)
7 (s)
C, (pgmm3)
498.9 997.6 332.6 249.4 486.7 475.1 464.1 453.6 498.9 498.9 498.9 498.9 498.9 498.9 498.9 498.9 498.9
2483.7 3967.2 1989.2 1741.9 2447.6 2413.2 2380.4 2349.1 2187.0 1989.2 1847.9 3569.4 2626.4 2844.7 1741.9 3225.6 3967.5
490.0 990.0 330.0 250.0 480.0 470.0 460.0 450.0 500.0 500.0 500.0 500.0 500.0 500.0 500.0 500.0 500.0
2484.0 3967.0 1989.0 1742.0 2448.0 2413.0 2380.0 2349.0 2187.0 1989.0 1848.0 3569.0 2617.0 2831.0 1769.0 3235.0 3967.0
1
0
20
40
60
00
100
120
140
160
180
200
Time (I)
Fig. 1. Effect of different initial concentrations (Ci) on the temporal variation of concentration calculated with the finite-volume model. Concentrations are calculated at the middle of the tunnel, x = 500 m, n = 1: (~)C~=1mgm~3;n=2:(~)C~=2mgm~3;n=3:(O)C~=3mgm~3;n=4:(~)C~=4mgm~3; n=5:(x)C:=5mgm-3;n=6:(-)Cp=6mgm-3.
long. Only CO emissions from one vehicle type (ECE
15-04, CC > 2000 cm3) entering from the left side of the tunnel have been considered. Vehicles travel at a speed of 100 km h-’ and with a number of 1000 vehicles h- ‘. The component of the wind speed along x is 2ms-‘, the external ones along y and z are 0 m s-l. Neutral conditions are assumed outside the tunnel (l/L = 0 m- ‘) and a roughness length of 0.1 m is used. The boundary condition farther from the tunnel is 2 mgm-3 and the initial condition is 0 mgmm3. Starting from this reference scenario, only one input variable per time has been changed in order
to understand its influence over the calculated concentrations. It is important to understand how much the initial concentration affects the calculated concentration inside the tunnel. Figure 1 shows the influence of different initial concentrations over the concentration calculated in the middle of the tunnel. The vertical axis reports the relative percent variation of concentration with respect to the one calculated with null initial concentration. The relative percent variation of concentration was calculated as lOO*(Cl - Co)/Co where Cf refers the nth initial concentration and
Modelling traffic air pollution Co refers the null initial concentration. Seven different initial concentrations have been used: 0, 1, 2, 3, 4, 5 and 6mgmm3, respectively. It is possible to note that there are big differences in the initial phases of the simulation but they tend to disappear as the simulation evolves. After about 2 min all the concentrations are practically equal, then it is possible to state that, after a reasonable time, the model results become independent of the initial concentration. However, the relative importance of the initial concentration will increase for a small number of vehicles inside the tunnel. Boundary conditions play an important role in the internal concentrations. Greater the concentration upwind to the tunnel, the greater is the mass that flows into the tunnel. It can be shown that if the boundary concentrations farther from the tunnel are increased, then the mean concentration inside the tunnel increases. The mean concentration is calculated by averaging the values of concentration over all the grid cells of the tunnel at a fixed time. In particular, the final value of the mean concentration is about linearly dependent on the boundary conditions. It is noteworthy that the downwind concentration outside the tunnel and near to it assumes very large values; this can indicate serious air pollution problems in places where tunnels are situated. In order to describe correctly the vehicle-induced turbulence an appropriate value of the empirical parameter c[ must be used. This value can be determined by means of a comparisons between observed and calculated concentrations; its determination could be the aim cd a next paper. At present it is useful to understand the role of CIover the calculated
0
a00
600
1545
concentration. Figure 2 shows that if the value of c(is increased, then the mean concentration inside the tunnel decreases since the total speed of the flow inside the tunnel grows. In fact, this result has been obtained with the wind blowing in the same direction in which vehicles are travelling. Figure 3 shows the mean concentration inside the tunnel when the vehicles go in the direction opposite to the wind direction. The situation is completely different from the one depicted in Fig. 2. When the value of tl increases the mean concentration initially increases due to the fact that the total speed of the flow inside the tunnel decreases. The concentration increases until TVreaches the value of 0.8, then it decreases as the total flow speed changes its direction and increases in absolute value. For tl equal to 1 the flow direction coincides with the vehicle’s direction. For c1equal to 0.6 and 0.8 the flow speed is so small that no steady state is reached during the time shown in Fig. 3. From these figures it is clear that the vehicle-induced turbulence plays a very important role over the dispersion of pollutants inside tunnels. Figure 4 finally illustrates the amount of mass exchanged between the two grid cells in the middle of the tunnel during a 3 h simulation; it has been obtained with the vehicles directed in the same direction of the wind. The total mass exchanged due to atmospheric turbulence is of the order of grams, it is small with respect to the one exchanged due to advection and vehicles-induced turbulence which is of the order of kilograms. Increasing the value of parameter CIthe piston effect becomes even more important than the advection. The mass exchanged for atmospheric turbulence and for advection decreases when c(increases only because the mass
am
1200
15M)
Time (8)
Fig. 2. Effect of the empirical parameter c(over the mean concentration as a function of time. The wind blows in the same direction in which vehiclesare travelling: (0) tl = 0, ( x ) a = 0.2,(0) c(= 0.4,( + ) c(= 0.6, (A) c( = 0.8, ( - ) tl = 1.
18W
1546
R. BELLASIO
I
0
300
600
900
1200
1500
1800
Time (8)
Fig. 3. Effect of the empirical parameter c( over the mean concentration as a function of time. The wind blows in the direction opposite to the one in which vehicles are travelling: (0) c( = 0, (x) ct= 0.2, (0) ct = 0.4, ( + ) c1= 0.6, (A) a = 0.8, ( - ) c(= 1.
Fig. 4. Total mass exchanged between the two grid cells at the middle of the tunnel due to turbulent atmospheric dispersion (0) , advection (0) and vehicle-induced turbulence (A) during a 3 h simulation. available to be transported by these two effects decreases. It is clear that c( does not have any direct influence over atmospheric turbulence and advection.
EXAMPLES
OF APPLICATION
This section describes two examples of application of the finite-volume model. The first one is the application of the model to a fictitious tunnel with different vehicles during five working days and the second is the application of the model inside a hypothetical underground railway during the same days.
A road tunnel 1 km long, 10 m wide and 5 m high has been considered; the tunnel has been divided into 100 grids 10 m long. An urban traffic regime with nine vehicle’s categories has been assumed: cars with different engines and different fuels, motorbikes, light-duty vehicles and heavy-duty vehicles. The number of vehicles per hour and per vehicle category, entering from each side of the tunnel is illustrated in Table 3; the first row of this table contains the identification number for the vehicle category. The vehicle’s speed has been put equal to 80 km h-’ from hour 1 to hour 7, from hour 10 to hour 12 and from hour 14 to hour 17; during the rest of the day it has been put to
Modelling traffic air pollution 60 km h- ’ due to the intense traffic of the rush hours.
The wind along x has been maintained constant at a speed of 2 m s-l, but blowing from left to right from 7 a.m. to 6 p.m. and from right to left during the rest of the day. The external wind speed along the y direction has been maintainled constant at 1 m s- ’ and the vertical component at 0 m s- ‘. The roughness length
Table 3. Number of vehicles entering in each side of the tunnel at each hour. R’owscontain the number of vehicles for each hour and columns contain different vehicle categories Vehicle ID Hour 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
16
17
18
25
27
29
30
32
34
48 29 14 8 8 19 53 87 89 69 76 82 85 68 78 77 85 99 106 97 91 74 67 67
38 23 12 6 6 15 43 70 71 55 61 66 68 54 62 62 68 79 84 78 73 59 54 54
10 6 3 2 2 4 11 17 18 14 15 16 17 14 16 15 17 20 21 19 18 15 13 13
13 8 4 2 2 5 15 24 25 19 21 23 24 19 22 22 24 27 29 27 25 21 19 19
6 3 2 1 1 2 6 10 11 8 9 10 10 8 9 9 10 12 13 12 11 9 8 8
7 4 2 1 1 3 9 31 38 46 38 39 16 24 40 47 46 34 19 13 12 8 7 7
14 9 4 3 2 5 14 25 28 36 23 22 9 15 32 33 25 18 17 23 22 19 17 17
6 4 2 1 1 3 7 21 19 7 22 6 17 12 7 8 5 12 17 13 12 9 8 8
1 1 1 0 0 2 5 12 13 4 2 6 10 10 3 3 2 4 6 4 4 2 2 2
9
10
1547
has been put at 0.1 m and the inverse of the Monin-Obukhov length has been maintained constant at 0.1 m- 1 for the whole day; this assumption is reasonable for a wintry day. Two sinks have been considered placed at x = 330 m and x = 670 m with a constant rate of aspiration of 5 m3 s- ’ for the whole day. The boundary concentrations farther from the tunnel and the external concentrations for the sinks are illustrated in Fig. 5 both for carbon monoxide and nitrogen oxides. Due to the different values of concentration, two axes are used in this figure, the left one for CO and the right one for NO,. Figures 6 and 7 show, respectively, the time series of CO and NO, concentration during the third day of simulation at x = 330 m and x = 670 m. For both the pollutants the concentrations are very high with respect to the ones usually measured in the atmosphere. They reach a maximum peak during the morning rush hourofabout19mgm-3forCOand3.5mgm-3for NO, at x = 670 m. During the evening rush hour maximum concentrations are reached at x = 330 m; they are about 19 mgmm3 for CO and 3 mgme3 for NO,. These different positions for the maximum concentrations between the morning and the evening are due to the different wind directions. It can be noted that when the wind is blowing from left to right (from 7 a.m. to 6 p.m.) the concentrations increase moving from left to right inside the tunnel. The opposite is true when the wind blows in the other direction. The minimum values are reached during the night at about 5 p.m., they are equal to about 1 mgmm3 for CO and 0.3 mgmW3 for NO,. Small concentrations during the night are due to a small number of vehicles travelling inside the tunnel and to small concentrations at the boundary. The strong oscillations of concentration
IWO
0 12
3
4
5
6
7
8
11
12
13
14
15
16
17
18
19
20
21
22
23
Time (h)
Fig. 5. Boundary concentrations far from the tunnel for CO ( x ) and NO,(O) used for the example of applications. Units for CO are reported on the left vertical axis and those for NO, are reported on the right vertical axis.
24
1548
R. BELLASIO
60 Time (h)
Fig. 6. Carbon monoxide concentrations
during the third day of simulation at x = 330m (A) and at x = 670 m (0).
70
72
Time (h)
Fig. 7. Nitrogen oxides concentrations
during the third day of simulation at x = 330 m (A) and at x = 670 m (0).
the night are due to the small number of vehicles travelling inside the tunnel. When the number of vehicles is small the piston effect is practically a discrete phenomenon: pollutant mass is pushed and pulled at a point when a vehicle passes at that point. When the number of vehicles increases the piston effect becomes more continuous, it becomes practically an extra wind to be added everywhere to the atmospheric one. Figure 8 illustrates the external concentrations near the tunnel entrances for CO and NO, during the third day of simulation. It can be seen that the external concentration next to the tunnel can be very high. For CO this concentration is about two times the during
boundary concentration farther from the tunnel and for NO, it is even high. These high values of concentrations are principally due to the pollutant mass exiting from the tunnel and, in part, to the mass emitted there. The finite-volume numerical model has been used to simulate the CO pollution due to external emissions in the underground railway. The underground trains do not emit any pollutant because they are electric powered, but their motion is responsible for the piston effect inside the tunnel. An underground tunnel 10 km long, 20 m wide and 5 m high has been considered. The tunnel length has been divided into 100 grids 100 m long. The initial concentration for
Modelling traffic air pollution
1549
14000
12000
1 -
pL P
8M)o
E :: $0
a,
-
2ooo
2z
0 4.6
&Ad-A-A 50
52
54
58
58
60
62
64
66
,
3
i2
Time (h)
Fig. 8. External concentrations of CO and NO, near to the tunnel entrances during the third day of simulation: (a) CO concentration to the left of the tunnel; ( + ) CO concentration to the right of the tunnel; ( - ) NO, concentration to the left of the tunnel; (A) NO, concentration to the right of the tunnel. 6wn
T
01: 46
1
50
52
54
50
58
60
62
64
86
68
70
72
Time (h)
Fig. 9. Carbon monoxide concentration during the third day of simulation inside the underground railway. Concentrations are calculated at x = 2 km (0) x = 6 km ( + ) and x = 8 km ( - ).
CO has been put at O~grn-~. A constant wind of 0.5 m s- ’ blowing from left to right along the x direction, an external wind speed of 0.5 ms-’ along the y direction and a null vertical wind have been considered. The inverse of the Monin-Obukhov length was of 0.1 m-‘. The train speed has been considered constant to 80 km h-’ and the number of trains per hour has been varied during the day. In particular, it has been assumed that no trains are present from hour 1 to hour 6,20 trains h- ’ at hour 7 and from hour 22 to hour 24 and 30 trains h-’ during the remaining hours. Four passenger entrances have been considered at 2,4,6 and 8 km; external air has been
assumed to enter with a constant flux of 15 m3 s- ’ in the first, the second and the fourth entrance and of 20 m3 s-l in the third one which was assumed to be the biggest one. External concentration relative to the entrances and boundary concentrations have been put at those shown in Fig. 5. It must be noticed that, if the external concentration is greater than the internal one, the entrances are sources of pollutants for the subway. A period of 5 days has been simulated. Figure 9 shows the temporal variation of CO concentration corresponding to the underground entrances for the third day of simulation at x = 2 km, x = 6 km and x = 8 km. The temporal variation, a part for the
1550
R. BELLASIO
IWO
0
IIr”
0
mm
2oca
3Qoo
4ooo
Em0
7mu
x(m)
Fig. 10. Carbon monoxide concentration inside the underground railway at three times of the third day of simulation: 03 : 00 (O), 09 : 00 ( + ) and 19 : 00 ( - ).
due to the piston effect, resembles the temporal variation of the external concentration with two peaks corresponding to the rush hours during the morning and the evening. It is worth noting that oscillations start after 6 a.m., when the trains start, which means that they are really due to the piston effect. The minimum concentrations are reached from 5 to 6 a.m., when the external concentration is also minimum. Figure 10 shows the CO concentration inside the underground as a function of the position at three fixed times corresponding to hours 03:00,09:00 and 19:00 of the third day of simulation. Concentrations are very different due to the different values of the external concentration. At 03:00, when no trains are running, there are no oscillations around the mean value of concentration since the piston effect is not present. Moreover, since the external CO concentration is about 0.8 mg m -3 and the internal one, before the first entrance, is about 1 mgmW3, the internal concentration at 3 a.m. corresponding to each entrance of the subway decreases. At 09:OO and 19:00 oscillations due to the piston effect are well evident. The external concentration at these hours is about 7 mgmm3 while the internal one is smaller; hence the internal concentration corresponding to each entrance increases. Concentrations at 19 are roughly constant while those at 9 are constant in the first 2 km inside the tunnel and then start to decrease. This fact could indicate that after some time something similar to a steady state is reached inside the whole tunnel. oscillations
CONCLUSION
Two models are developed for the simulation of pollutants due to traffic in road tunnels: a box-model
and a finite-volume model. An ad hoc formulation is used for the turbulence due to the vehicle’s motion inside the tunnel. The boundary conditions are maintained constant only farther from the tunnel. Close to the tunnel a three-dimensional equation of conservation of mass is solved in order to consider the perturbation induced by the high values of the internal concentration. Both the formulation of the piston effect and the determination of the external concentration close to the tunnel are first working hypothesis and need to be improved. The numerical scheme of solution of the partial differential equations is tested with the box model. A very good agreement between the numerical solution and the analytical solution is obtained, both in terms of value of concentration and of times at which they are reached. Sensitivity analysis, carried out using the finite-volume model, has demonstrated that it displays a correct physical behaviour to the variation of the input parameters. Examples of application are illustrated for a road tunnel and for a subway. Very high concentrations are found as a result of the scarce dispersion and, for the road tunnel, of the big quantity of pollutants emitted. If the working hypothesis used in this work for the description of the vehicle-induced turbulence is correct, measures of concentration at a point inside the tunnel should give very different values even for relatively close times. The comparison of the model’s results against experimental data will be the subject of future work.
Acknowledgements-The author wishes to thank Dr M. Tamponi of PMIP-USSL 16 (Lecco) and Dr R. Bianconi of IDEA snc (Milan) for their useful advice, and Dr A. N. Skouloudis of the Joint Research Centre of Ispra (Varese) for the useful discussions about modular programming.
Modelling traffic air pollution REFERENCES
Eggleston H. S., Gaudioso D., Gorissen N., Joumard R., Rijkeboer R. C., Samaras Z. and Zierock K. H. (1991) CORINAIR working group on emission factors for calculating 1990 emissions from road traffic. Vol. 1: methodology and emission factors. Final Report. Ermak L. E. and Merry M. H. (1988) A methodology for evaluating heavy-gas dispersion models. LLNL, contract W-7405-ENG-48. Kagawa J. (1984) Health effects of air pollutants and their management. Atmospheric Environment 18, 613-620. Lanzani G. and Tamponi M. (1995) Microscale Lagrangian particle model for 1he dispersion of primary pollutants in a street canyon. Sensitivity analysis and first validation trials. Atmospheric Environment 29, 346553475. Patankar S. V. (1980) Numerical Heat Transfer and Fluid Flow. McGraw-Hill, New York.
1551
Rotach M. W. (1995) Profiles of turbulence statistics in and above an urban street canyon. Atmospheric Environment 29, 1473-1486.
Staehelin J., Schlapfer K., Biirgin T., Steinemann U., Schneider S., Brunner D., Baumle M., Meier M., Zahner C., Keiser S., Stahel W. and Keller C. (1995) Emission factors from road traffic from a tunnel study (Gubrist tunnel, Switzerland). Part 1: concept and first result. Sci Total Envir. (169) 141-147.
WHO (1987) Air quality guidelines for Europe. World Health Organisation Regional Publications, European Ser. No. 23. Yamartino R. J. and Wiegand G. (1986) Development and evaluation of simple models for the flow, turbulence and pollutant concentration fields within an urban street canyon. Atmospheric Environment 20, 2137-2156. Zannetti P. (1990) Air Pollution Modeling. Theories, Computational Methods and Available Software. Computational mechanics publications. Southampton, Boston.
Atmospheric
Pm: S1352-2310(96)00296-g
MODELLING
TRAFFIC
AIR POLLUTION
ROBERTO
Environment
Vol. 31, No. 10, pp. 1539-1551, 1997 0 1997 Else&r Science Ltd AU rights reserved. Printed in Great Britain 1352-2310/97 $17.00 + 0.00
IN ROAD TUNNELS
BELLASIO
Scientific Consultant, Via Dolomiti 8, 20017 Rho (MI), Italy (First received 20 December 1995 and in final form 13 July 1996. Published March 1997)
Abstract-This paper presents two models for the description of air pollutant concentrations in road tunnels due ts> traffic. Turbulence is assumed to depend on both atmospheric turbulence and vehicle’s motion. Emissions are calculated with the COPERT90 methodology (Eggleston et al., 1991). Up to 34 different vehicle categories can be considered at the same time, with an arbitrary number of vehicles travelling inside the tunnel. Emissions are calculated as a function of the position inside the tunnel and of the time. Three pollutants can be simulated with the current version of the models, CO, NO, and VOC. It is also possible to consider vehicles with null emissions. The models are able to consider the effects of an arbitrary number of sinks. Flow rates and outdoor concentrations are a function of the time for each sink. The equation of conservation of mass has been solved with the control volumes method. Particular attention has been given to the formulation of stability conditions. Sensitivity analysis was conducted to verify the model answer to different input parameters such as initial concentration, boundary concentrations and vehicle-induced turbulence. Examples of application are given for a tunnel with urban traffic regime, including passenger cars with different fuels, light duty trucks, heavy duty trucks and motorbikes, and for an underground railway. Simulations have been carried out for the five working days of the week. 0 1997 Elsevier Science Ltd. All rights reserved. Key word index: Tunnel, finite volumes, piston effect.
INTRODUCTION
The interest about air pollution has grown enormously in the last decades due to a better comprehension of the dangerous effects that some pollutants have on human health (Kagawa, 1984; WHO, 1987). Air pollution in urban areas is due to different sources, but the major responsible source is vehicle traffic. Road-traffic-related pollution should be measured at locations where maximum concentrations are expected (EU directive, 85/203/EEC). These locations can be determined by means of appropriate models capable of describing the atmospheric dispersion in particular geometries where transport processes are not well established. An example of these geometries is the street canyon 1:hat has been extensively studied in the last years botlh from an experimental (Rotach, 1995) and a modelling point of view (Yamartino and Wiegand, 1986; Lanzani and Tamponi, 1995). Other locations where high concentrations can be expected are the tunnels in which there are high emissions of pollutants .in relatively small volumes and possible scarce venlilation. Even if people do not spend very long time,3 in tunnels these high concentrations can be extremely dangerous. For example, carbon monoxide has an attitude about 200 times greater
than oxygen to link with haemoglobin (Kagawa, 1984), forming COHb in a short time. This paper describes the modelling of air pollution due to traffic inside tunnels. The author does not know any other work on the simulation of pollutant dispersion inside tunnels. For this reason, some ad hoc working hypotheses are adopted; they can be verified with validation trials and, if necessary, they can be improved in future works. The development of two Eulerian models is illustrated. The first one is a box model that gives the concentration only as function of time; the second one has been developed with the finite-volume method and is able to give concentrations as a function of time and of position inside the tunnel. In both the models, turbulence is treated considering both atmospheric turbulence and turbulence due to vehicle’s movement. Emissions are calculated with the COPERT90 methodology (Eggleston et al., 1991) as functions of vehicle speed for 34 different vehicle types. The numerical model is able to consider also vehicles with null emissions such as, for example, a train inside the underground tunnel. The next section describes the theoretical formulation which is the same for both the models. Then the numerical formulation for the box model and for the one using the finite-volume
1539
R. BELLASIO
1540
methodology is illustrated. A section is dedicated to the description of emission calculation. Finally, the last two sections present the sensitivity analysis and examples of application.
THEORY
The basic equation solved by the models is the equation of conservation of mass inside the tunnel. It is assumed that the mixing along the directions normal to the tunnel axis (y and z) is effective enough due to the vehicle’s motion so that the pollutants are well mixed. This hypothesis is supported by the work of Staehelin et al. (1995). Hence, the concentration is assumed to be a function of time t and of the distance x inside the tunnel from the left entrance, C = C(x, t). The one-dimensional equation of conservation of mass for a turbulent flow obtained starting from the usual equation of conservation of mass and applying the Reynolds’ decomposition of the variables, is the following:
aed -z+ax=-ax ac
ach
+s-P
(11
where C represents the (mean) concentration (pgm-3), u is the wind speed inside the tunnel (m s- ‘), S = S(x, t) is the source of pollutant (pgmm3 s-i) and P = P(x, t) is the sink term (pgrne3 s-i). The source term is a function of position and time because vehicles move in the tunnel causing emissions in different locations and different times. The source formulation will be discussed later. The model is able to consider two types of sinks, fans and aspirators. Fans catch polluted air inside the tunnel and push it out. Aspirators catch “clear” air outside the tunnel and push it into the tunnel; this last type of sink can simulate openings along the tunnel. The wind speed u represents the mean wind component in the tunnel. It is considered positive when the wind blows from left to right and negative when the wind blows in the opposite direction. The term C’U’ represents the turbulent flux of mass inside the tunnel. Often in the lower layer of the atmosphere this term is satisfactorily described by means of the gradient theory. Turbulence inside the tunnels is not only due to the atmospheric turbulence but also due to the vehicle’s motion; in this work an ad hoc formulation of the turbulent term is adopted. It is assumed that the turbulent mass flux is the sum of two terms, one related to the atmospheric turbulence and the other to the vehicle’s motion: -7-t
C u = T~tm+ TV,,,
where T indicates a turbulent term. The atmospheric turbulent flux T,,, is described by means of the gradient theory (Zannetti, 1990). It is assumed to be proportional to the concentration gradient via the turbulent dispersion coefficient K,.
The turbulent dispersion coefficient can be calculated as a function of the friction velocity u* and of the Monin-Obukhov length I. (Zannetti, 1990). The horizontal dispersion coefficient K, is proportional to the vertical dispersion coefficient K,, K, = /lKL. Since no radiation hits the ground at the centre of the tunnel, neutral stability is assumed, the resulting K,, averaged over the tunnel cross-section, is K = k/W& X 2
(3)
where H, is the tunnel height (m), k is the von Karman constant (k = 0.4) and u* is the friction velocity (m s- ‘). The turbulent dispersion coefficients outside the tunnel are calculated with the similarity theory. In the rest of the tunnel the coefficients are calculated interpolating from the external values to the internal one given by equation (3). The turbulent flux due to the vehicle’s motion, TVeh, is described assuming that each vehicle determines, with its movement, the movement of a quantity of air with a speed that is proportional to the ratio between the vehicle cross-section and the tunnel cross-section (piston effect). The turbulent flux of mass at a point will be the product between this piston velocity and the concentration at that point. Hence, the turbulent transport of mass is described with C’u’= - K, 2
+ aC F
ceil
uVeh
(4)
where CIis an empirical constant, A,,,, is the crosssection of a particular vehicle type (m’), Ace,, is the tunnel cross-section (m2) and uVehis the vehicle speed (ms-‘). Substituting equation (4) into equation (l), we have
ac x
acu
+dx
a =-z
qCuve,-K~~
>
+S-P
(5) where rl is the following non-dimensional &$ cell
quantity: (6)
The vehicle-induced turbulence is then treated as an advection term; in other words, the vehicle’s motion contributes to increase or decrease the wind speed inside the tunnel. Equation (5) can be solved provided the initial concentration and boundary conditions are specified. The external concentration near the tunnel entrances is influenced by the high value of concentration inside the tunnel. For this reason, the boundary concentration is assumed to be constant only farther from the tunnel. Near the tunnel the boundary concentration is calculated as a function of the time solving a three-dimensional conservation equation in a box. Equation (5) has been solved in two different ways, with a box-model formulation assuming that the pollutant is well mixed inside the whole tunnel (C = C(t))
Modelling traffic air pollution and with a one-dimensional formulation (C = C(x, t)).
THE BOX-MODEL
finite-volume numerical
concentration C exits from it. It must be noted that when C,,, is greater than C then P* is a source term. The general numerical solution of equation (8) is Max( -w,O)
FORMULATION
The relatively simple formulation of the box-model helps to understand the physical phenomena responsible for the temporal variation of concentration inside the tunnel. Let us consider a tunnel characterised by a length L (m) and a cross-section A (m’), the concentration at the left boundary of the tunnel is called CBL@gmm3) and the one at the right boundary is called CBR(pg m - ‘). The integration of equation (5) over the tunnel volume gives the following equation stating the integral mass conservation:
AL; =A
[
C(u + rpV,h) -
-A w
Kg
1 Left
1
ac + V%eh) - Kxz Righ,
[
+ s* - p*.
(7)
Equation (7) states that the temporal variation of mass is due to the difference between the mass entering or exiting in the tunnel from the left and the mass exiting or entering in the tunnel from the right, plus the difference between the mass “produced” inside the tunnel and the mass exiting from the sinks. The source term S* and the sink term P* have the dimensions of mass per unit of time (pgs-‘); they are, respectively, equal to the source term S and the sink term P times the tunnel volume AL. The numerical formulation of advection and dispersion in equation. (7) is the following: A,!,:
1541
= A[CBL lMax(w, 0) - C Max( - w, 0)]
+ Max(w, 0) + T )) +-
+
CaLAt
L C&t L
+ -$(S*
- R*(C” - co,,))
(11)
where At is the discrete temporal step, C” is the concentration at time t” = n At and C”+ ’is the concentration at time t”+ ’ = (n + 1)At. Equation (8) can be solved iteratively provided that an initial condition Ci is furnished. In order to achieve numerical stability, equation (8) is solved with a temporal step At such that
L
At <
(12) Max( - w, 0) + Max(w, 0) + %
+ T
It is obtained forcing the coefficient of C” in equation (11) to be positive. It can be seen that this condition is not very critical; as an example, for a tunnel 1 km long, a speed w of 5 m s- ‘, a turbulent atmospheric coefficient of 2 m2 s-l, a sink rate R* of 10 m3 SK’and a tunnel cross-section A of 50 m2 it is obtained a At of about 200 s.
THE ANALYTICAL
SOLUTION
OF THE BOX-MODEL
FORMULATION
- A [C Max(w, 0) - CBRMax( - w, 0)] Under special conditions it is possible to give an analytical solution of equation (8). Some of these L conditions are not very realistic but they are useful +S* --P* (8) because the comparison between the analytical solution and the numerical solution gives an idea of the where the speed w (rn s- ‘) is given by goodness of the numerical scheme adopted. If w > 0 and CaL = CaR = C,,, = Ca then equation w = u + fluve,,. (9) (8) assumes the form The temporal derivative is discretised by means of the finite-difference method. The sink term P* can be AL$f=ACB W+%i-T described by the following relation: ( > -
C-
AK,rp
CBL
Cm-C +AK,----L
Pl:= R*(C- C,",)
(10)
where R* (m3 s-‘) is the rate at which air is exchanged between the tunnel and the external environment. C,,, (pgme3) is the external pollutant concentration near the sinks. Equation (10) indicates that when a volume of air with external concentration Goutenters into the tunnel, the same volume of air with internal
+S* >
(13)
whose solution, assuming that w, K,,S*,R* are constants and Ci = 0 is c(t)=C(m)(l-exp(
-4))
(14)
R. BELLASIO
1542
where the asymptotic concentration
C(c0) is given by
s*
C(m) = CB + Aw+~
2AK,
(15) $R*
and where the characteristic time r is given by L
Z= w+y
2K,
R*’ +A
(16)
It can be seen that when t $ r, equation (14) tends to the value C(a) , after t x 57 the concentration C(t) is greater than the 99.3% of the final concentration C(a). The characteristic time r represents the time after which the concentration is equal to the 63.2% of C(a). The characteristic time r increases with the tunnel length L and the tunnel cross-section A, while it decreases with the advection speed w, the dispersion coefficient K, and the sink rate R*. The final concentration C(a) increases with the source term S, the external concentration Ca and the tunnel length L, while it decreases with the sink rate R*, the advection speed w, the atmospheric dispersion coefficient K, and the tunnel cross-section A.
speed due to the vehicle’s motion, and is described by equation (9). It can be seen that the advection term is described by the up-wind or down-wind formulation as a function of the wind direction, and the dispersion term is described with the finite-difference method. Applying the finite-difference method to the temporal derivative of the concentration in equation (17), the concentration at the time step (n + l)At can be written in terms of the concentration at time step n At. Hence, starting from an initial condition Ci it is possible to calculate the concentrations at each time step. The stability condition for equation (17) has been found to be Ax
At< Kx, + KX+, + y Ax
+ Max(w, 0) + Max( - w, 0) (19)
The value of At is calculated using the values that maximise the denominator of equation (19). Another stability condition arises from the three-dimensional equation of conservation of mass solved outside the tunnel to obtain the external concentration near its entrances. This condition is practically the sum of three terms similar to the one in equation (19), one for each direction.
THE FINITE-VOLUME FORMULATION
Equation (5) is solved with the finite-volume method (Patankar, 1980) which guarantees a perfect mass conservation. The tunnel is divided into grid boxes along the x direction and the equation is integrated over each of these boxes. The following integral equation of conservation of mass is obtained:
THE EMISSION SUBMODEL
Emission factors for traffic were calculated as functions of vehicle speed with the COPERT90 methodology (Eggleston et al., 1991). Only hot engines emissions were considered. Vehicles were grouped in categories on the basis of fuel type, the age of the vehicle, the engine cylinder capacity and the weight. AAX: = A(@i-Qi+l) + AAXSi A total of 34 vehicle categories has been considered, each identified by a number as illustrated in Table 1. - A AX Ri(Ci - CpU’) (17) The term PC indicates passenger cars, LDT indicates light duty tracks, HDT indicates heavy duty tracks where Ax is the grid size (m), @i(pgm-‘s- ‘) represand 2W indicates two wheelers. Subscripts indicate ents the mass flux at the position (i - 1)Ax and the fuel type, g stands for gasoline, d stands for diesel cDi+i(~grn-~ s- ‘) represents the mass flux at the posiand lpg stands for liquid petroleum gas. Passenger tion iAx. These fluxes describe the turbulent dispersion and the advection of mass. As an example the flux cars are also split according to the different steps of international legal conformity (ECE classes). Buses Qi at position “i” is presented: are described as heavy-duty trucks. @i = Cj_iMax(w,O) - CjMax(-w,O) Vehicle’s fluxes in the tunnel must be specified at each hour of simulation. For each vehicle type _K cj-cj-l W-9 the average number of vehicles per unit time cp x, Ax ’ (vehicles h-i) and the average speed must be specified. Index “i” represents the positions of the bounding Vehicles entering from left and from right are treated surface of the grid box while index j represents the independently. The temporal distance between two calculation position; each calculation position is in vehicles of the same type in the same hour remains the middle of a grid box. Then, for example, fluxes are constant. It is calculated as r, = 3600/p and is exevaluated at the grid box faces and concentrations are pressed in seconds. Vehicles are treated as mobile evaluated at the middle of each grid box. The velocity sources, this means that, considering for simplicity term w is the sum of the component of the atmoonly one vehicle in the tunnel, the emissions are alspheric wind speed inside the tunnel and the turbulent ways zero but not in the positions occupied by the
Modelling traffic air pollution
1543
partial differential equations, a comparison between the numerical solution of the box model and the analytical solution has been carried out. If this comCylinder parison gives good results, it is reasonable to think ID Vehicle Category capacity that even the numerical scheme adopted for the soluPRE ECE cc< 1.4 1 tion with the finite-volume methodology is correct. In peg 2 PRE ECE 1.4
Table 1. Vehicle categories for which emission factors are calculated
1544
R. BELLASIO Table 2. Comparison between the analytical solution and the numerical solution of the box model Input
w (ms-‘)
K, (m’s_‘)
R* (m3s-r)
2.3 1.2 3.4 4.5 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3
0 0 0 0 2 4 6 8 0 0 0 0 0 0 0 0 0
2 1 3 4 2 2 2 2 2 2 2 2 2 2 2 2 2
1600
Analytical
A (m2) F (vehh-‘) 40 40 40 40 40 40 40 40 50 60 70 40 40 40 40 40 40
Uveh(kmh-i)
100 100 100 100 100 100 100 100 100 100 100 100 100 100 50 150 200
100 100 100 100 100 100 100 100 100 100 100 50 75 125 100 100 100
Numerical
r (s)
C, (pgme3)
7 (s)
C, (pgmm3)
498.9 997.6 332.6 249.4 486.7 475.1 464.1 453.6 498.9 498.9 498.9 498.9 498.9 498.9 498.9 498.9 498.9
2483.7 3967.2 1989.2 1741.9 2447.6 2413.2 2380.4 2349.1 2187.0 1989.2 1847.9 3569.4 2626.4 2844.7 1741.9 3225.6 3967.5
490.0 990.0 330.0 250.0 480.0 470.0 460.0 450.0 500.0 500.0 500.0 500.0 500.0 500.0 500.0 500.0 500.0
2484.0 3967.0 1989.0 1742.0 2448.0 2413.0 2380.0 2349.0 2187.0 1989.0 1848.0 3569.0 2617.0 2831.0 1769.0 3235.0 3967.0
1
0
20
40
60
00
100
120
140
160
180
200
Time (I)
Fig. 1. Effect of different initial concentrations (Ci) on the temporal variation of concentration calculated with the finite-volume model. Concentrations are calculated at the middle of the tunnel, x = 500 m, n = 1: (~)C~=1mgm~3;n=2:(~)C~=2mgm~3;n=3:(O)C~=3mgm~3;n=4:(~)C~=4mgm~3; n=5:(x)C:=5mgm-3;n=6:(-)Cp=6mgm-3.
long. Only CO emissions from one vehicle type (ECE
15-04, CC > 2000 cm3) entering from the left side of the tunnel have been considered. Vehicles travel at a speed of 100 km h-’ and with a number of 1000 vehicles h- ‘. The component of the wind speed along x is 2ms-‘, the external ones along y and z are 0 m s-l. Neutral conditions are assumed outside the tunnel (l/L = 0 m- ‘) and a roughness length of 0.1 m is used. The boundary condition farther from the tunnel is 2 mgm-3 and the initial condition is 0 mgmm3. Starting from this reference scenario, only one input variable per time has been changed in order
to understand its influence over the calculated concentrations. It is important to understand how much the initial concentration affects the calculated concentration inside the tunnel. Figure 1 shows the influence of different initial concentrations over the concentration calculated in the middle of the tunnel. The vertical axis reports the relative percent variation of concentration with respect to the one calculated with null initial concentration. The relative percent variation of concentration was calculated as lOO*(Cl - Co)/Co where Cf refers the nth initial concentration and
Modelling traffic air pollution Co refers the null initial concentration. Seven different initial concentrations have been used: 0, 1, 2, 3, 4, 5 and 6mgmm3, respectively. It is possible to note that there are big differences in the initial phases of the simulation but they tend to disappear as the simulation evolves. After about 2 min all the concentrations are practically equal, then it is possible to state that, after a reasonable time, the model results become independent of the initial concentration. However, the relative importance of the initial concentration will increase for a small number of vehicles inside the tunnel. Boundary conditions play an important role in the internal concentrations. Greater the concentration upwind to the tunnel, the greater is the mass that flows into the tunnel. It can be shown that if the boundary concentrations farther from the tunnel are increased, then the mean concentration inside the tunnel increases. The mean concentration is calculated by averaging the values of concentration over all the grid cells of the tunnel at a fixed time. In particular, the final value of the mean concentration is about linearly dependent on the boundary conditions. It is noteworthy that the downwind concentration outside the tunnel and near to it assumes very large values; this can indicate serious air pollution problems in places where tunnels are situated. In order to describe correctly the vehicle-induced turbulence an appropriate value of the empirical parameter c[ must be used. This value can be determined by means of a comparisons between observed and calculated concentrations; its determination could be the aim cd a next paper. At present it is useful to understand the role of CIover the calculated
0
a00
600
1545
concentration. Figure 2 shows that if the value of c(is increased, then the mean concentration inside the tunnel decreases since the total speed of the flow inside the tunnel grows. In fact, this result has been obtained with the wind blowing in the same direction in which vehicles are travelling. Figure 3 shows the mean concentration inside the tunnel when the vehicles go in the direction opposite to the wind direction. The situation is completely different from the one depicted in Fig. 2. When the value of tl increases the mean concentration initially increases due to the fact that the total speed of the flow inside the tunnel decreases. The concentration increases until TVreaches the value of 0.8, then it decreases as the total flow speed changes its direction and increases in absolute value. For tl equal to 1 the flow direction coincides with the vehicle’s direction. For c1equal to 0.6 and 0.8 the flow speed is so small that no steady state is reached during the time shown in Fig. 3. From these figures it is clear that the vehicle-induced turbulence plays a very important role over the dispersion of pollutants inside tunnels. Figure 4 finally illustrates the amount of mass exchanged between the two grid cells in the middle of the tunnel during a 3 h simulation; it has been obtained with the vehicles directed in the same direction of the wind. The total mass exchanged due to atmospheric turbulence is of the order of grams, it is small with respect to the one exchanged due to advection and vehicles-induced turbulence which is of the order of kilograms. Increasing the value of parameter CIthe piston effect becomes even more important than the advection. The mass exchanged for atmospheric turbulence and for advection decreases when c(increases only because the mass
am
1200
15M)
Time (8)
Fig. 2. Effect of the empirical parameter c(over the mean concentration as a function of time. The wind blows in the same direction in which vehiclesare travelling: (0) tl = 0, ( x ) a = 0.2,(0) c(= 0.4,( + ) c(= 0.6, (A) c( = 0.8, ( - ) tl = 1.
18W
1546
R. BELLASIO
I
0
300
600
900
1200
1500
1800
Time (8)
Fig. 3. Effect of the empirical parameter c( over the mean concentration as a function of time. The wind blows in the direction opposite to the one in which vehicles are travelling: (0) c( = 0, (x) ct= 0.2, (0) ct = 0.4, ( + ) c1= 0.6, (A) a = 0.8, ( - ) c(= 1.
Fig. 4. Total mass exchanged between the two grid cells at the middle of the tunnel due to turbulent atmospheric dispersion (0) , advection (0) and vehicle-induced turbulence (A) during a 3 h simulation. available to be transported by these two effects decreases. It is clear that c( does not have any direct influence over atmospheric turbulence and advection.
EXAMPLES
OF APPLICATION
This section describes two examples of application of the finite-volume model. The first one is the application of the model to a fictitious tunnel with different vehicles during five working days and the second is the application of the model inside a hypothetical underground railway during the same days.
A road tunnel 1 km long, 10 m wide and 5 m high has been considered; the tunnel has been divided into 100 grids 10 m long. An urban traffic regime with nine vehicle’s categories has been assumed: cars with different engines and different fuels, motorbikes, light-duty vehicles and heavy-duty vehicles. The number of vehicles per hour and per vehicle category, entering from each side of the tunnel is illustrated in Table 3; the first row of this table contains the identification number for the vehicle category. The vehicle’s speed has been put equal to 80 km h-’ from hour 1 to hour 7, from hour 10 to hour 12 and from hour 14 to hour 17; during the rest of the day it has been put to
Modelling traffic air pollution 60 km h- ’ due to the intense traffic of the rush hours.
The wind along x has been maintained constant at a speed of 2 m s-l, but blowing from left to right from 7 a.m. to 6 p.m. and from right to left during the rest of the day. The external wind speed along the y direction has been maintainled constant at 1 m s- ’ and the vertical component at 0 m s- ‘. The roughness length
Table 3. Number of vehicles entering in each side of the tunnel at each hour. R’owscontain the number of vehicles for each hour and columns contain different vehicle categories Vehicle ID Hour 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
16
17
18
25
27
29
30
32
34
48 29 14 8 8 19 53 87 89 69 76 82 85 68 78 77 85 99 106 97 91 74 67 67
38 23 12 6 6 15 43 70 71 55 61 66 68 54 62 62 68 79 84 78 73 59 54 54
10 6 3 2 2 4 11 17 18 14 15 16 17 14 16 15 17 20 21 19 18 15 13 13
13 8 4 2 2 5 15 24 25 19 21 23 24 19 22 22 24 27 29 27 25 21 19 19
6 3 2 1 1 2 6 10 11 8 9 10 10 8 9 9 10 12 13 12 11 9 8 8
7 4 2 1 1 3 9 31 38 46 38 39 16 24 40 47 46 34 19 13 12 8 7 7
14 9 4 3 2 5 14 25 28 36 23 22 9 15 32 33 25 18 17 23 22 19 17 17
6 4 2 1 1 3 7 21 19 7 22 6 17 12 7 8 5 12 17 13 12 9 8 8
1 1 1 0 0 2 5 12 13 4 2 6 10 10 3 3 2 4 6 4 4 2 2 2
9
10
1547
has been put at 0.1 m and the inverse of the Monin-Obukhov length has been maintained constant at 0.1 m- 1 for the whole day; this assumption is reasonable for a wintry day. Two sinks have been considered placed at x = 330 m and x = 670 m with a constant rate of aspiration of 5 m3 s- ’ for the whole day. The boundary concentrations farther from the tunnel and the external concentrations for the sinks are illustrated in Fig. 5 both for carbon monoxide and nitrogen oxides. Due to the different values of concentration, two axes are used in this figure, the left one for CO and the right one for NO,. Figures 6 and 7 show, respectively, the time series of CO and NO, concentration during the third day of simulation at x = 330 m and x = 670 m. For both the pollutants the concentrations are very high with respect to the ones usually measured in the atmosphere. They reach a maximum peak during the morning rush hourofabout19mgm-3forCOand3.5mgm-3for NO, at x = 670 m. During the evening rush hour maximum concentrations are reached at x = 330 m; they are about 19 mgmm3 for CO and 3 mgme3 for NO,. These different positions for the maximum concentrations between the morning and the evening are due to the different wind directions. It can be noted that when the wind is blowing from left to right (from 7 a.m. to 6 p.m.) the concentrations increase moving from left to right inside the tunnel. The opposite is true when the wind blows in the other direction. The minimum values are reached during the night at about 5 p.m., they are equal to about 1 mgmm3 for CO and 0.3 mgmW3 for NO,. Small concentrations during the night are due to a small number of vehicles travelling inside the tunnel and to small concentrations at the boundary. The strong oscillations of concentration
IWO
0 12
3
4
5
6
7
8
11
12
13
14
15
16
17
18
19
20
21
22
23
Time (h)
Fig. 5. Boundary concentrations far from the tunnel for CO ( x ) and NO,(O) used for the example of applications. Units for CO are reported on the left vertical axis and those for NO, are reported on the right vertical axis.
24
1548
R. BELLASIO
60 Time (h)
Fig. 6. Carbon monoxide concentrations
during the third day of simulation at x = 330m (A) and at x = 670 m (0).
70
72
Time (h)
Fig. 7. Nitrogen oxides concentrations
during the third day of simulation at x = 330 m (A) and at x = 670 m (0).
the night are due to the small number of vehicles travelling inside the tunnel. When the number of vehicles is small the piston effect is practically a discrete phenomenon: pollutant mass is pushed and pulled at a point when a vehicle passes at that point. When the number of vehicles increases the piston effect becomes more continuous, it becomes practically an extra wind to be added everywhere to the atmospheric one. Figure 8 illustrates the external concentrations near the tunnel entrances for CO and NO, during the third day of simulation. It can be seen that the external concentration next to the tunnel can be very high. For CO this concentration is about two times the during
boundary concentration farther from the tunnel and for NO, it is even high. These high values of concentrations are principally due to the pollutant mass exiting from the tunnel and, in part, to the mass emitted there. The finite-volume numerical model has been used to simulate the CO pollution due to external emissions in the underground railway. The underground trains do not emit any pollutant because they are electric powered, but their motion is responsible for the piston effect inside the tunnel. An underground tunnel 10 km long, 20 m wide and 5 m high has been considered. The tunnel length has been divided into 100 grids 100 m long. The initial concentration for
Modelling traffic air pollution
1549
14000
12000
1 -
pL P
8M)o
E :: $0
a,
-
2ooo
2z
0 4.6
&Ad-A-A 50
52
54
58
58
60
62
64
66
,
3
i2
Time (h)
Fig. 8. External concentrations of CO and NO, near to the tunnel entrances during the third day of simulation: (a) CO concentration to the left of the tunnel; ( + ) CO concentration to the right of the tunnel; ( - ) NO, concentration to the left of the tunnel; (A) NO, concentration to the right of the tunnel. 6wn
T
01: 46
1
50
52
54
50
58
60
62
64
86
68
70
72
Time (h)
Fig. 9. Carbon monoxide concentration during the third day of simulation inside the underground railway. Concentrations are calculated at x = 2 km (0) x = 6 km ( + ) and x = 8 km ( - ).
CO has been put at O~grn-~. A constant wind of 0.5 m s- ’ blowing from left to right along the x direction, an external wind speed of 0.5 ms-’ along the y direction and a null vertical wind have been considered. The inverse of the Monin-Obukhov length was of 0.1 m-‘. The train speed has been considered constant to 80 km h-’ and the number of trains per hour has been varied during the day. In particular, it has been assumed that no trains are present from hour 1 to hour 6,20 trains h- ’ at hour 7 and from hour 22 to hour 24 and 30 trains h-’ during the remaining hours. Four passenger entrances have been considered at 2,4,6 and 8 km; external air has been
assumed to enter with a constant flux of 15 m3 s- ’ in the first, the second and the fourth entrance and of 20 m3 s-l in the third one which was assumed to be the biggest one. External concentration relative to the entrances and boundary concentrations have been put at those shown in Fig. 5. It must be noticed that, if the external concentration is greater than the internal one, the entrances are sources of pollutants for the subway. A period of 5 days has been simulated. Figure 9 shows the temporal variation of CO concentration corresponding to the underground entrances for the third day of simulation at x = 2 km, x = 6 km and x = 8 km. The temporal variation, a part for the
1550
R. BELLASIO
IWO
0
IIr”
0
mm
2oca
3Qoo
4ooo
Em0
7mu
x(m)
Fig. 10. Carbon monoxide concentration inside the underground railway at three times of the third day of simulation: 03 : 00 (O), 09 : 00 ( + ) and 19 : 00 ( - ).
due to the piston effect, resembles the temporal variation of the external concentration with two peaks corresponding to the rush hours during the morning and the evening. It is worth noting that oscillations start after 6 a.m., when the trains start, which means that they are really due to the piston effect. The minimum concentrations are reached from 5 to 6 a.m., when the external concentration is also minimum. Figure 10 shows the CO concentration inside the underground as a function of the position at three fixed times corresponding to hours 03:00,09:00 and 19:00 of the third day of simulation. Concentrations are very different due to the different values of the external concentration. At 03:00, when no trains are running, there are no oscillations around the mean value of concentration since the piston effect is not present. Moreover, since the external CO concentration is about 0.8 mg m -3 and the internal one, before the first entrance, is about 1 mgmW3, the internal concentration at 3 a.m. corresponding to each entrance of the subway decreases. At 09:OO and 19:00 oscillations due to the piston effect are well evident. The external concentration at these hours is about 7 mgmm3 while the internal one is smaller; hence the internal concentration corresponding to each entrance increases. Concentrations at 19 are roughly constant while those at 9 are constant in the first 2 km inside the tunnel and then start to decrease. This fact could indicate that after some time something similar to a steady state is reached inside the whole tunnel. oscillations
CONCLUSION
Two models are developed for the simulation of pollutants due to traffic in road tunnels: a box-model
and a finite-volume model. An ad hoc formulation is used for the turbulence due to the vehicle’s motion inside the tunnel. The boundary conditions are maintained constant only farther from the tunnel. Close to the tunnel a three-dimensional equation of conservation of mass is solved in order to consider the perturbation induced by the high values of the internal concentration. Both the formulation of the piston effect and the determination of the external concentration close to the tunnel are first working hypothesis and need to be improved. The numerical scheme of solution of the partial differential equations is tested with the box model. A very good agreement between the numerical solution and the analytical solution is obtained, both in terms of value of concentration and of times at which they are reached. Sensitivity analysis, carried out using the finite-volume model, has demonstrated that it displays a correct physical behaviour to the variation of the input parameters. Examples of application are illustrated for a road tunnel and for a subway. Very high concentrations are found as a result of the scarce dispersion and, for the road tunnel, of the big quantity of pollutants emitted. If the working hypothesis used in this work for the description of the vehicle-induced turbulence is correct, measures of concentration at a point inside the tunnel should give very different values even for relatively close times. The comparison of the model’s results against experimental data will be the subject of future work.
Acknowledgements-The author wishes to thank Dr M. Tamponi of PMIP-USSL 16 (Lecco) and Dr R. Bianconi of IDEA snc (Milan) for their useful advice, and Dr A. N. Skouloudis of the Joint Research Centre of Ispra (Varese) for the useful discussions about modular programming.
Modelling traffic air pollution REFERENCES
Eggleston H. S., Gaudioso D., Gorissen N., Joumard R., Rijkeboer R. C., Samaras Z. and Zierock K. H. (1991) CORINAIR working group on emission factors for calculating 1990 emissions from road traffic. Vol. 1: methodology and emission factors. Final Report. Ermak L. E. and Merry M. H. (1988) A methodology for evaluating heavy-gas dispersion models. LLNL, contract W-7405-ENG-48. Kagawa J. (1984) Health effects of air pollutants and their management. Atmospheric Environment 18, 613-620. Lanzani G. and Tamponi M. (1995) Microscale Lagrangian particle model for 1he dispersion of primary pollutants in a street canyon. Sensitivity analysis and first validation trials. Atmospheric Environment 29, 346553475. Patankar S. V. (1980) Numerical Heat Transfer and Fluid Flow. McGraw-Hill, New York.
1551
Rotach M. W. (1995) Profiles of turbulence statistics in and above an urban street canyon. Atmospheric Environment 29, 1473-1486.
Staehelin J., Schlapfer K., Biirgin T., Steinemann U., Schneider S., Brunner D., Baumle M., Meier M., Zahner C., Keiser S., Stahel W. and Keller C. (1995) Emission factors from road traffic from a tunnel study (Gubrist tunnel, Switzerland). Part 1: concept and first result. Sci Total Envir. (169) 141-147.
WHO (1987) Air quality guidelines for Europe. World Health Organisation Regional Publications, European Ser. No. 23. Yamartino R. J. and Wiegand G. (1986) Development and evaluation of simple models for the flow, turbulence and pollutant concentration fields within an urban street canyon. Atmospheric Environment 20, 2137-2156. Zannetti P. (1990) Air Pollution Modeling. Theories, Computational Methods and Available Software. Computational mechanics publications. Southampton, Boston.