- Email: [email protected]
Realistic numerical simulation of chlorine decay in pipes
PII: S0043-1354(98)00107-9
Wat. Res. Vol. 32, No. 11, pp. 3307±3312, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0043-1354/98 $19.00 + 0.00
REALISTIC NUMERICAL SIMULATION OF CHLORINE DECAY IN PIPES OSMAN N. OZDEMIR* and A. METIN GER Department of Civil Engineering, Middle East Technical University (METU), 06531 Ankara, Turkey (First received April 1997; accepted in revised form February 1998) AbstractÐThe objective of the present study is to develop a numerical mass transport model that can simulate chlorine disappearance mechanisms under low ¯ow conditions. The model accounts for axial convection and radial diusion in a pipe under turbulent, transitional and laminar ¯ow regimes and uses realistic velocity pro®les. The bulk decay coecient assumed to be constant in previous models is taken to be a function of residence time of water in pipes based on the ®ndings of the experiments by laboratory tests. Furthermore, a methodology is introduced to convert simple one dimensional model solutions into re®ned results as if they were calculated by a much more precise two dimensional model. # 1998 Elsevier Science Ltd. All rights reserved Key wordsÐchlorine, decay, diusion, numerical simulation, water supply network
y=distance from the wall (=r0ÿr); Z=dimensionless parameter (=WdL/n); n=kinematic viscosity of water (L2Tÿ1).
NOMENCLATURE A0=dimensionless parameter (=LDr/r20U); A1=dimensionless parameter (=kL/U); A2=dimensionless parameter (=Wdr0/Dr); c=point value of chlorine concentration in the cross-section (MLÿ3); C0=steady state average inlet chlorine concentration (MLÿ3); Cout=steady state average outlet chlorine concentration (MLÿ3); Ct=chlorine concentration at any time t (MLÿ3); C=dimensionless point value of chlorine concentration in the cross-section (=c/C0); Ccup=dimensionless average chlorine concentration of the pipe at any cross-section; Cki=dimensionless average end concentrations for the use of an inlet bulk decay; Ckx=dimensionless average end concentrations for the use of a variable bulk decay; C1D=dimensionless average concentration for 1D solution; C2D=approximate dimensionless average concentration for 2D solution; CT=exact (numerical) dimensionless average concentration for 2D solution; d=pipe diameter (L); Dr=apparent (eective) diusivity of chlorine in the radial direction (L2Tÿ1); f=friction factor; Fv=nondimensional ¯ow parameter depending on the velocity pro®le; FC=correction factor for 1D solution; k=decay constant for bulk ¯ow reaction (Tÿ1); ks=equivalent sand grain roughness (L); L=length of pipe (L); R=dimensionless radial distance from pipe center (=r/r0); r0=pipe radius (L); r=radial distance from pipe center (L); Re=Reynolds number (=Ud/n); U=pipe average velocity (LTÿ1); Wd=chlorine consumption rate at the pipe wall (LTÿ1); V=velocity at any radial location in a pipe section (LTÿ1); X=dimensionless distance along the pipe (=x/L); x=distance along the pipe (L);
INTRODUCTION
*Author to whom all correspondence should be addressed. [Tel: +90-312-4269040; Fax: +90-312-2101262, E-mail: [email protected]].
The emphasis of the Safe Drinking Water Act and its amendments eective in the U.S. is the attainment of standards at the point of consumption. This has led to greater attention to an understanding of the processes eecting the deterioration of water quality in a distribution system. Drinking water is traditionally disinfected at the treatment plant. Because the disinfectant decays as water ¯ows through the pipes, water may receive additional disinfection at some critical locations. Chlorine, which is the most commonly used type of disinfectant, is a highly unstable oxidizing agent that may react with organics, inorganics, as well as the pipe material in the network. The interaction of chlorine with natural organics may lead to formation of trihalomethanes and other by-products that are suspected to be carcinogens. Chlorine is known to be consumed in the bulk liquid phase and at the pipe wall. As indicated by LeChevallier et al. (1988), transport of the disinfectant from the bulk liquid phase into the bio®lm at the walls is an important factor in understanding chlorine decay rates. Higher disinfectant concentrations are required to inactivate bacterial populations in the bio®lm as compared to those suspended in the bulk liquid phase. Several models have been proposed to simulate decay behavior of chlorine in pipes and distribution networks. Hunt and Kroon (1991) developed a network model for
3307
3308
Osman N. Ozdemir and A. Metin Ger Table 1. The ¯ow parameter Fv classi®ed according to Reynolds number Reynolds number Re < 2300 2300 < Re < 30000 Re>30000
Fv 2
2(1 ÿ R ) f 1/2(2.15 log(y/r0) + 1.43) + 1 1
chlorine residual that used a ®rst-order decay reaction with a unique rate constant for each pipe. Biswas et al. (1993) model for chlorine decay within single pipes under steady-state ¯ow conditions included both bulk-¯ow reaction and radial diusion, and subsequent pipe-wall reaction of chlorine. Though Biswas et al. (1993) accounted for radial transport of chlorine, due to restrictions to obtain an analytical solution for the problem, their work was limited to fully turbulent ¯ows only. Rossman et al. (1994) used a lumped mass-transfer coecient to account for the radial transport. While such an approach derived from heat transfer analogies (Skelland, 1974; Edwards et al., 1976) is valid under dierent ¯ow regimes (turbulent and laminar ¯ow), radial velocity variation and bulk decay changes have not been examined. During the operation of water supply networks, low ¯ow (i.e. transitional or laminar) conditions are frequently encountered; prevalent in dead end sections (Buchberger and Wu, 1995). Such low ¯ow conditions corresponding to longer residence times cause loss of chlorine more than that of for turbulent ¯ow conditions. Therefore, it is imperative that accurate estimation of chlorine decay is obtained for both transitional and laminar ¯ow conditions. The objective of the present study is to develop a numerical mass transport model that can simulate chlorine disappearance mechanisms under low ¯ow conditions. The model accounts for axial convection and radial diusion in a pipe under turbulent, transitional and laminar ¯ow regimes and uses realistic velocity pro®les. The bulk decay coecient assumed to be constant in previous models is taken to be a function of residence time of water in pipes based on the ®ndings of the experiments by laboratory tests. Furthermore, a methodology is introduced to convert simple one dimensional model solutions into re®ned results as if they were calculated by a much more precise two dimensional model. MODEL AND THE SOLUTION METHOD
Considering previous works by Hunt et al. (1988); Wable et al. (1991); Sharp et al. (1991); Biswas et al. (1993); Rossman et al. (1994) it is possible to assume that the disappearance of chlor-
Reference
Equation
Rouse, 1949 Rouse, 1949 Hinze, 1975
4 5 6
ine ¯owing through a pipe is governed by ®rstorder kinetics. There are two types of reactions occurring in the pipe causing this disappearance. The ®rst one is the reaction of chlorine with dissolved organic and inorganic matters within the bulk ¯ow. The second type occurs on the pipe wall. The rates of these two reactions can be dierent depending on the characteristics of the treated water, age and type of pipe and ¯ow regime in the pipes. With the above assumptions, the steady state two dimensional conservation of mass equation for dilute concentration of total free chlorine in water ¯owing through a pipe, in the absence of axial dispersion, is V
@ c Dr @ @c r ÿ kc @x r @r @r
1
subject to boundary conditions: at x = 0, c = C0; at r = 0, @c/@r = 0 (symmetry); at r = r0, Dr @c/ @r = ÿ Wd c (wall condition). Where V is the actual velocity pro®le as a function of radial position in the pipe; c is the chlorine concentration in the bulk ¯ow; x is the distance along the pipe; Dr is the apparent (eective) diusivity of chlorine in water in the r direction; r is the radial distance from the center of a pipe; r0: radius of the pipe; k is the decay rate constant in the bulk ¯ow; C0 is the concentration of chlorine at inlet of the pipe; Wd is the constant describing consumption rate of chlorine on the pipe wall. The term on the left hand side accounts for the convective ¯ux of chlorine through the section. The ®rst term on the right hand side accounts for the radial diusion of chlorine towards the walls of the pipe due to wall consumption. The second term on the right hand side represents decay in the bulk ¯ow. The nondimensionalized form of this equation is given below: @ C A0 @ @C Fv 2 R ÿ A1 C R @R @X @R subject to boundary conditions: at X = 0, C = 1; at R = 0, @C/@R = 0; at R = 1, @C/@R = ÿ A2C, where nondimensional parameters are R = r/r0, C = c/C0, X = x/L, A0=L Dr/r20U, A1=k L/U, A2=Wd r0/Dr and the nondimensional velocity Fv
Table 2. Friction factor (f) for any pipe wall type Wall Type
f
Reference
Equation
Smooth Transition Rough
0.316/Re0.25 1/{(1.14 ÿ 2 log(ks/d) + (21.25/Re0.9)}2 1/{(2 log(r0/ks) + 1.75}
Shames, 1992 Shames, 1992 Rouse, 1949
7 8 9
Realistic numerical simulation of chlorine decay in pipes
in equation 2 de®ned as: Fv V=U,
3
where V = velocity at any radial position in the pipe; U = average velocity. Fv assumes the following forms as tabulated in Table 1 for dierent ¯ow regimes. In Table 1, y is the distance from the pipe wall and equal to r0ÿr. The friction factor f appearing in Eq. (5) in Table 1 assumes the forms for dierent pipe wall types as tabulated in Table 2. In this Table 1 ks is the equivalent sand grain roughness; d is the pipe diameter and Re is Reynolds number. A Crank±Nicolson implicit ®nite dierence scheme is applied to integrate the dimensionless form of the governing equation 2; the ®nite dierence form is given as equation 10 below. In equation 10 subscript i denotes the position of the node in the radial direction; i = 1 and i = imax are corresponding to center and wall nodes, respectively. Similarly, k is used to denote the position of the node in the axial direction; k = 1 and k = kmax are for the inlet and outlet sections, respectively. Fv
Ci,k1 ÿ Ci,k DX
Ci,k1 ÿ Ci,k DX
A0
0
0
EXPERIMENTAL PROCEDURE AND FINDINGS
In all the models developed to date the bulk decay rate coecients have been assumed to be a constant. Yet, controlled experiments performed at the USEPA laboratories in Cincinnati to measure possible bulk decay parameter changes along the axis of a pipe proved otherwise. These experiments revealed the fact that chlorine is consumed more under low ¯ow conditions (Table 3), and that the chlorine content of water drops as the residence time is
Ci max,k ÿ Ci maxÿ1,k 11 ÿA2 Ci max,k DR While assessing the value of A2, for Reynolds number less than 10,000, Dr was calculated using the graph prepared by Bischo and Levenspiel (1963) and for larger Reynolds numbers the expression Dr=0.02 Udf 0.5 suggested by Hinze (1975) was used. On the pipe center where R = 0, there is an indeterminacy of type 0/0. L'Hospital's rule is used to eliminate it; the implicit dierence equation becomes:
equations were solved by means of the Thomas algorithm (Chapra and Canale 1990). The run-time of the Fortran program developed for the implementation of the algorithm on a Pentium computer, is in seconds. The solution of the dierence equations lead to nodal concentrations. The crosssectional cup mixing average concentrations were also obtained by applying Simpson's rule to equation 13 below: 1 1 2pRC dR 2pR dR 13 Ccup
Ci1,k1 ÿ 2Ci,k1 Ciÿ1,k1 Ci1,k ÿ 2Ci,k Ciÿ1,k DR2 DR2 A0 Ci1,k1 ÿ Ciÿ1,k1 Ci1,k ÿ Ciÿ1,k A1 ÿ Ci,k1 Ci,k 2R 2DR 2DR 2
A0 2
Equation 10 is valid for all nodes except for those on the pipe wall and those on the pipe axis. On the pipe wall using the boundary condition @C/ @R = ÿ A2C the dierence equation becomes:
Fv
3309
increased in the pipe. The test pipe used in the experiments is 24.3 m. (80 ft) long, and 15 cm. (6 inch) diameter unlined cast iron pipe. In order to measure continuous inlet and outlet chlorine concentrations, two Hach Cl 17 analytical devices measuring free chlorine were installed at the inlet and the outlet sections. The procedure described below was used to determine bulk decay coecients. For each run corresponding to a dierent ¯ow regime, 15 samples of 120 ml were taken both from the inlet and the outlet sections. Samples were stored in amber colored bottles and sealed with Te¯on. The amber color kept the interference of light away. The Te¯on cover was selected to avoid reaction of chlorine with other possible sealing materials. Furthermore, Te¯on is also known as hydrophobic and does not react with chlorine. The bottles were kept at constant temperature of 148C in a refrigerator. Then one bottle a day for 15 consecutive days was sampled for its free chlorine content using EPA approved DPD±FAS method. The same procedure was repeated for
Ci1,k1 ÿ 2Ci,k1 Ciÿ1,k1 Ci1,k ÿ 2Ci,k Ciÿ1,k DR2 DR2
The recursive equation 10 for intermediate points coupled with equations 11 and 12 constitute a set of simultaneous dierence equations. These set of equations were solved for various combinations of dimensionless coecients covering a wide range of geometry and ¯ow regimes. A mesh size of 100 nodes in the radial and 2000 nodes in the axial direction was used. The resulting system of
10
ÿ
A1 Ci,k1 Ci,k 2
12
dierent ¯ow regimes in the pipe. Figure 1. and Fig. 2. include data points obtained from the repeated application of this procedure to dierent ¯ow regimes for the bulk decay at the inlet and outlet sections, respectively. The following ®rst order decay rate equation (equation 14) is used to determine the bulk decay coecient k using the data given in Fig. 1. and Fig. 2. Ct C0 exp ÿkt
14
3310
Osman N. Ozdemir and A. Metin Ger Table 3. Test pipe chlorine decay for dierent ¯ows
Flow GPM (103 m3/s) 0.3 (0.02) 0.5 (0.03) 1.0 (0.06) 1.8 (0.11) 3.0 (0.19) 3.6 (0.23) 4.1 (0.26) 6.8 (0.43) 13.0 (0.82)
Reynolds number 175 292 583 1049 1749 2099 2390 3964 7579
Residence time (h) 6.54 3.92 1.96 1.09 0.65 0.55 0.48 0.29 0.15
Comparison of these ®gures revealed that while there is a single value for bulk decay constant at the inlet section, the value of the constant at the outlet section depends on the ¯ow rate. For the outlet, bulk decay constants are higher for lower ¯ow rates. The results of the bottle decay tests for the samples collected from the inlet yielded a bulk decay constant of 1.0 10ÿ3/h. On the other hand, the same tests for the samples taken from the outlet revealed that there are dierent bulk decay constants ranging from 0.8 10ÿ3/h to 5.0 10ÿ3/h depending on the rate of ¯ow in the pipe. These test results were also con®rmed with heterotrophic plate count (HPC) tests performed at the USEPA in Cincinnati. For the experiments run to obtain the bulk decay constants, extra samples were collected from the inlet and outlet sections. These samples were analyzed by HPC tests; the results are reported in Table 4. The number of colony forming units (CFU) for the samples from the outlet section were found to be very high compared to those for the inlet as the ¯ow rate decreases. These results are in agreement with the results obtained from the bulk decay tests (Fig. 1 and Fig. 2).
Inlet concentration (C0) (mg/l) 0.66 0.71 0.60 0.73 0.79 0.82 0.81 0.90 0.94
Outlet concentration (Cout) (mg/l) 0.46 0.51 0.45 0.66 0.73 0.75 0.77 0.87 0.93
Cout/C0 0.70 0.72 0.75 0.90 0.92 0.92 0.95 0.97 0.99
COMPARISON OF THE PRESENT MODEL WITH THE AVAILABLE ONES
As depicted in Figs 1 and 2 bulk decay constant k is not a ®xed value along the axis of the pipe. The assumption that the k is constant commonly employed in the previous studies is therefore questionable. To demonstrate this eect, a series of numerical runs were made. The numerical simulation depicted in Fig. 3. accentuated the signi®cance of the observation made. Dierences in pipe-end chlorine concentrations for an inlet k and variable k increases drastically with increase in L/d ratio. In Fig. 3, on y axis, Cki denotes for the dimensionless average end concentrations obtained from the application of a single inlet bulk decay constant for the whole length of the pipe whereas, Ckx denotes the one obtained from a linearly changing bulk decay from an inlet to an outlet value. Furthermore, bulk decay alone did not account for the disappearance of total chlorine in the test pipe. To be able to determine the wall decay conTable 4. Heterotrophic plate count (HPC) test results Flow GPM (lt/s)
Number of organisms inlet (CFU/ml)
0.5 (0.03) 1.0 (0.06) 2.5 (0.16) 15.0 (0.95)
6 7 4 3
outlet (CFU/ml) 240 74 23 12
Fig. 1. Bulk decay for the test pipe inlet.
Fig. 2. Bulk decay for the test pipe outlet.
Fig. 3. The eect of bulk decay change on pipe-end concentrations.
Realistic numerical simulation of chlorine decay in pipes
3311
rection factor FC, is introduced such that the 1D solution can be corrected to yield precise prediction of chlorine concentration decay along the pipe axis. The 1D chlorine decay can be represented by C1D
Ct expfÿ k Wd =r0 gt C0
15
The true decay CT is approximated by the present 2D solution, C2D. CT 1C2D Fig. 4. Comparison of the alternative models with the experimental results.
stant representing the decay of chlorine on the pipe wall a dierent series of numerical runs were made of the model. In these numerical simulations the variation of the bulk decay is assumed to be linear; for the bulk decay, the average of the inlet and outlet section was implemented. The wall decay constant value thus obtained is Wd=1.0 10ÿ6 m/s. It corresponds to minimum deviation between the experimental observations and numerical results. The models of Biswas et al. (1993) and of Rossman et al. (1994) were also compared with the experimental results for the outlet cup mixing concentrations. Figure 4 is provided to compare these results with the numerical model output presented in this study. The superiority of the proposed model over the previous ones is evident. The models considered, namely Biswas', Rossman's and the present one, are all in good agreement with the experiments for Reynolds number greater than 10,000. This is mainly because, no detectable chlorine loss can be observed for a 24 m pipe length for Reynolds number above 10,000. As the Reynolds number gets smaller, the Biswas' analytical method and Rossman's plug ¯ow approach yield large deviations from the experimental observations. In general the disagreement between Biswas' model predictions and the experiments are due to the prevailing assumption such as the one about the velocity distribution. Furthermore, the mass-transfer coecient used in Rossman's model was obtained from an experiment designed for dissolution of benzoic acid, cinnamic acid and B-naphtanol coated tube wall through a thin concentration boundary layer. This derivation requires the assumption of in®nite concentration at the wall and zero concentration at the outer edge of a very thin concentration boundary layer.
The aforementioned correction factor FC relates C1D and C2D such that FC
C1D ÿ C2D C2D
17
In other words, if FC and C1D are known, without going through the cumbersome analysis CT can be estimated by CT 1
C1D 1 FC
18
The numerical outputs of equation 15 and the present model were used to determine the functional dependence of FC on the pipe, ¯uid and ¯ow parameters. It can be shown by dimensional analysis that FC is the following form: FC fnc Z, Re
19
where the newly established dimensionless parameter Z is de®ned as Z = Wd L/n. The variation of FC thus obtained is shown in Fig. 5. Using the best ®t lines shown in Fig. 5, the data can be condensed such that FC/Z is related to Reynolds number as shown in Fig. 6. In Fig. 6 the best ®t line (equation 20) corresponding to various ranges of Reynolds numbers are also included. log FC 100=Z ÿ0:9538 log Re=1000 ÿ 1 20 Either Fig. 5 or equation 20 may be used to predict chlorine decay with ease.
THE CORRECTION FACTOR
The present model, even though produces the best ®t to the experimental observations so far, requires a signi®cant computational eort. Therefore, to ease the utility of the method a cor-
16
Fig. 5. The correction factor (FC).
3312
Osman N. Ozdemir and A. Metin Ger
improved chlorine decay simulation signi®cantly. To ease the utility of the model a correction factor is introduced to enhance the use of the simple 1D solution. AcknowledgementsÐSupport given for the preparation of this paper by the Turkish National Science Foundation is acknowledged. The Director of the Drinking Water Research Division Dr Robert M. Clark is also acknowledged for his permission to use the laboratories at EPA in Cincinnati. We also like to thank Jack Tueschler and Chuck Guion of the U.S. Environmental Protection Agency for their sincere helps during the process of collecting the data in the laboratory. Fig. 6. Compact form of equation 19. CONCLUSIONS
It is known that in the design of distribution networks the secondary pipe sizes are chosen according to the emergency ®re demands rather than the domestic consumption. This leads to oversized pipe diameters. Under normal operating conditions very slow velocities can be observed in many parts of the system. This is more severe during night time operations. Furthermore, in either hydraulic or quality solutions of distribution networks continuity has to be assumed at the junction nodes in the system. Therefore, it can readily be accepted that, pipes are connected segments of an entire system rather than individual elements. Thus, the simple one dimensional solution in some of the pipes in the network may cause considerable errors in the system. These errors accumulate at each junction leading to large errors at the end points. The present study shows that as the ¯ow regime in a pipe becomes transitional and then laminar the bulk decay rate constant along the pipe increases. The bulk decay rate increase is observed to be due to increase in number of living organisms consuming chlorine in water. It is therefore concluded that the use of a single bulk decay value for an entire network, especially the value obtained at the water treatment plant euent may not be realistic. It is suggested that measuring the bulk decay at several critical locations be made in order to de®ne site speci®c decay coecients. Furthermore, the water quality models should also be capable of incorporating the changes in bulk decay constants. The present model, incorporating variable bulk decay and using realistic velocity distribution,
REFERENCES
Bischo K. B. and Levenspiel O. (1963) Advances in Chemical Engineering, Vol. 4. Academic Press, New York. Biswas P., Lu C. S. and Clark R. M. (1993) Chlorine concentration decay in pipes. Water Res. 27(12), 1715± 1724. Buchberger S. G. and Wu L. (1995) Model for instantaneous residential water demands. J. Hydrol. Eng. 121(3), 232±246. Chapra S. C. and Canale R. P. (1990) Numerical Methods for Engineers, 2nd ed. McGraw-Hill, International Editions. Edwards D. K., Denny V. E. and Mills A. F. (1976) Transport processes. McGraw-Hill, New York. Hinze J. O. (1975) Turbulence. McGraw-Hill, New York. Hunt W. A. and Kroon J. R. (1991) Model calibration for chlorine residuals in distribution systems. Proc. Water Quality Modeling in Distribution Systems. AWWA Research Foundation, Denver, CO. LeChevallier M. W., Cawthon C. and Lee R. G. (1988) Inactivation of bio®lm bacteria. Appl. Environ. Microbiol. 54, 2492±2499. Rossman L. A., Clark R. M. and Grayman W. M. (1994) Modeling chlorine residuals in drinking-water distribution systems. J. Environ. Eng. 120(4), 803±820. Rouse H. (1949) Elementary Fluid Mechanics. John Wiley and Sons, New York. Sharp W. W., Pfeer J. and Morgan M. (1991) In situ chlorine decay rate testing. Proc. Water Quality Modeling in Distribution Systems. AWWA Research Foundation, Denver, CO. Shames I. H. (1992) Mechanics of Fluids, third edition. McGraw Hill, New York. Skelland A. H. P. (1974) Diusional Mass Transfer. John Wiley and Sons, New York. Wable O., Dumoutier N., Duguet J. P., Jarrige P. A., Gelas G. and Depierre J. F. (1991) Modeling chlorine concentration in a network and applications to Paris distribution network. Proc. Water Quality Modeling in Distribution Systems. AWWA Research Foundation, Denver, CO.
Wat. Res. Vol. 32, No. 11, pp. 3307±3312, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0043-1354/98 $19.00 + 0.00
REALISTIC NUMERICAL SIMULATION OF CHLORINE DECAY IN PIPES OSMAN N. OZDEMIR* and A. METIN GER Department of Civil Engineering, Middle East Technical University (METU), 06531 Ankara, Turkey (First received April 1997; accepted in revised form February 1998) AbstractÐThe objective of the present study is to develop a numerical mass transport model that can simulate chlorine disappearance mechanisms under low ¯ow conditions. The model accounts for axial convection and radial diusion in a pipe under turbulent, transitional and laminar ¯ow regimes and uses realistic velocity pro®les. The bulk decay coecient assumed to be constant in previous models is taken to be a function of residence time of water in pipes based on the ®ndings of the experiments by laboratory tests. Furthermore, a methodology is introduced to convert simple one dimensional model solutions into re®ned results as if they were calculated by a much more precise two dimensional model. # 1998 Elsevier Science Ltd. All rights reserved Key wordsÐchlorine, decay, diusion, numerical simulation, water supply network
y=distance from the wall (=r0ÿr); Z=dimensionless parameter (=WdL/n); n=kinematic viscosity of water (L2Tÿ1).
NOMENCLATURE A0=dimensionless parameter (=LDr/r20U); A1=dimensionless parameter (=kL/U); A2=dimensionless parameter (=Wdr0/Dr); c=point value of chlorine concentration in the cross-section (MLÿ3); C0=steady state average inlet chlorine concentration (MLÿ3); Cout=steady state average outlet chlorine concentration (MLÿ3); Ct=chlorine concentration at any time t (MLÿ3); C=dimensionless point value of chlorine concentration in the cross-section (=c/C0); Ccup=dimensionless average chlorine concentration of the pipe at any cross-section; Cki=dimensionless average end concentrations for the use of an inlet bulk decay; Ckx=dimensionless average end concentrations for the use of a variable bulk decay; C1D=dimensionless average concentration for 1D solution; C2D=approximate dimensionless average concentration for 2D solution; CT=exact (numerical) dimensionless average concentration for 2D solution; d=pipe diameter (L); Dr=apparent (eective) diusivity of chlorine in the radial direction (L2Tÿ1); f=friction factor; Fv=nondimensional ¯ow parameter depending on the velocity pro®le; FC=correction factor for 1D solution; k=decay constant for bulk ¯ow reaction (Tÿ1); ks=equivalent sand grain roughness (L); L=length of pipe (L); R=dimensionless radial distance from pipe center (=r/r0); r0=pipe radius (L); r=radial distance from pipe center (L); Re=Reynolds number (=Ud/n); U=pipe average velocity (LTÿ1); Wd=chlorine consumption rate at the pipe wall (LTÿ1); V=velocity at any radial location in a pipe section (LTÿ1); X=dimensionless distance along the pipe (=x/L); x=distance along the pipe (L);
INTRODUCTION
*Author to whom all correspondence should be addressed. [Tel: +90-312-4269040; Fax: +90-312-2101262, E-mail: [email protected]].
The emphasis of the Safe Drinking Water Act and its amendments eective in the U.S. is the attainment of standards at the point of consumption. This has led to greater attention to an understanding of the processes eecting the deterioration of water quality in a distribution system. Drinking water is traditionally disinfected at the treatment plant. Because the disinfectant decays as water ¯ows through the pipes, water may receive additional disinfection at some critical locations. Chlorine, which is the most commonly used type of disinfectant, is a highly unstable oxidizing agent that may react with organics, inorganics, as well as the pipe material in the network. The interaction of chlorine with natural organics may lead to formation of trihalomethanes and other by-products that are suspected to be carcinogens. Chlorine is known to be consumed in the bulk liquid phase and at the pipe wall. As indicated by LeChevallier et al. (1988), transport of the disinfectant from the bulk liquid phase into the bio®lm at the walls is an important factor in understanding chlorine decay rates. Higher disinfectant concentrations are required to inactivate bacterial populations in the bio®lm as compared to those suspended in the bulk liquid phase. Several models have been proposed to simulate decay behavior of chlorine in pipes and distribution networks. Hunt and Kroon (1991) developed a network model for
3307
3308
Osman N. Ozdemir and A. Metin Ger Table 1. The ¯ow parameter Fv classi®ed according to Reynolds number Reynolds number Re < 2300 2300 < Re < 30000 Re>30000
Fv 2
2(1 ÿ R ) f 1/2(2.15 log(y/r0) + 1.43) + 1 1
chlorine residual that used a ®rst-order decay reaction with a unique rate constant for each pipe. Biswas et al. (1993) model for chlorine decay within single pipes under steady-state ¯ow conditions included both bulk-¯ow reaction and radial diusion, and subsequent pipe-wall reaction of chlorine. Though Biswas et al. (1993) accounted for radial transport of chlorine, due to restrictions to obtain an analytical solution for the problem, their work was limited to fully turbulent ¯ows only. Rossman et al. (1994) used a lumped mass-transfer coecient to account for the radial transport. While such an approach derived from heat transfer analogies (Skelland, 1974; Edwards et al., 1976) is valid under dierent ¯ow regimes (turbulent and laminar ¯ow), radial velocity variation and bulk decay changes have not been examined. During the operation of water supply networks, low ¯ow (i.e. transitional or laminar) conditions are frequently encountered; prevalent in dead end sections (Buchberger and Wu, 1995). Such low ¯ow conditions corresponding to longer residence times cause loss of chlorine more than that of for turbulent ¯ow conditions. Therefore, it is imperative that accurate estimation of chlorine decay is obtained for both transitional and laminar ¯ow conditions. The objective of the present study is to develop a numerical mass transport model that can simulate chlorine disappearance mechanisms under low ¯ow conditions. The model accounts for axial convection and radial diusion in a pipe under turbulent, transitional and laminar ¯ow regimes and uses realistic velocity pro®les. The bulk decay coecient assumed to be constant in previous models is taken to be a function of residence time of water in pipes based on the ®ndings of the experiments by laboratory tests. Furthermore, a methodology is introduced to convert simple one dimensional model solutions into re®ned results as if they were calculated by a much more precise two dimensional model. MODEL AND THE SOLUTION METHOD
Considering previous works by Hunt et al. (1988); Wable et al. (1991); Sharp et al. (1991); Biswas et al. (1993); Rossman et al. (1994) it is possible to assume that the disappearance of chlor-
Reference
Equation
Rouse, 1949 Rouse, 1949 Hinze, 1975
4 5 6
ine ¯owing through a pipe is governed by ®rstorder kinetics. There are two types of reactions occurring in the pipe causing this disappearance. The ®rst one is the reaction of chlorine with dissolved organic and inorganic matters within the bulk ¯ow. The second type occurs on the pipe wall. The rates of these two reactions can be dierent depending on the characteristics of the treated water, age and type of pipe and ¯ow regime in the pipes. With the above assumptions, the steady state two dimensional conservation of mass equation for dilute concentration of total free chlorine in water ¯owing through a pipe, in the absence of axial dispersion, is V
@ c Dr @ @c r ÿ kc @x r @r @r
1
subject to boundary conditions: at x = 0, c = C0; at r = 0, @c/@r = 0 (symmetry); at r = r0, Dr @c/ @r = ÿ Wd c (wall condition). Where V is the actual velocity pro®le as a function of radial position in the pipe; c is the chlorine concentration in the bulk ¯ow; x is the distance along the pipe; Dr is the apparent (eective) diusivity of chlorine in water in the r direction; r is the radial distance from the center of a pipe; r0: radius of the pipe; k is the decay rate constant in the bulk ¯ow; C0 is the concentration of chlorine at inlet of the pipe; Wd is the constant describing consumption rate of chlorine on the pipe wall. The term on the left hand side accounts for the convective ¯ux of chlorine through the section. The ®rst term on the right hand side accounts for the radial diusion of chlorine towards the walls of the pipe due to wall consumption. The second term on the right hand side represents decay in the bulk ¯ow. The nondimensionalized form of this equation is given below: @ C A0 @ @C Fv 2 R ÿ A1 C R @R @X @R subject to boundary conditions: at X = 0, C = 1; at R = 0, @C/@R = 0; at R = 1, @C/@R = ÿ A2C, where nondimensional parameters are R = r/r0, C = c/C0, X = x/L, A0=L Dr/r20U, A1=k L/U, A2=Wd r0/Dr and the nondimensional velocity Fv
Table 2. Friction factor (f) for any pipe wall type Wall Type
f
Reference
Equation
Smooth Transition Rough
0.316/Re0.25 1/{(1.14 ÿ 2 log(ks/d) + (21.25/Re0.9)}2 1/{(2 log(r0/ks) + 1.75}
Shames, 1992 Shames, 1992 Rouse, 1949
7 8 9
Realistic numerical simulation of chlorine decay in pipes
in equation 2 de®ned as: Fv V=U,
3
where V = velocity at any radial position in the pipe; U = average velocity. Fv assumes the following forms as tabulated in Table 1 for dierent ¯ow regimes. In Table 1, y is the distance from the pipe wall and equal to r0ÿr. The friction factor f appearing in Eq. (5) in Table 1 assumes the forms for dierent pipe wall types as tabulated in Table 2. In this Table 1 ks is the equivalent sand grain roughness; d is the pipe diameter and Re is Reynolds number. A Crank±Nicolson implicit ®nite dierence scheme is applied to integrate the dimensionless form of the governing equation 2; the ®nite dierence form is given as equation 10 below. In equation 10 subscript i denotes the position of the node in the radial direction; i = 1 and i = imax are corresponding to center and wall nodes, respectively. Similarly, k is used to denote the position of the node in the axial direction; k = 1 and k = kmax are for the inlet and outlet sections, respectively. Fv
Ci,k1 ÿ Ci,k DX
Ci,k1 ÿ Ci,k DX
A0
0
0
EXPERIMENTAL PROCEDURE AND FINDINGS
In all the models developed to date the bulk decay rate coecients have been assumed to be a constant. Yet, controlled experiments performed at the USEPA laboratories in Cincinnati to measure possible bulk decay parameter changes along the axis of a pipe proved otherwise. These experiments revealed the fact that chlorine is consumed more under low ¯ow conditions (Table 3), and that the chlorine content of water drops as the residence time is
Ci max,k ÿ Ci maxÿ1,k 11 ÿA2 Ci max,k DR While assessing the value of A2, for Reynolds number less than 10,000, Dr was calculated using the graph prepared by Bischo and Levenspiel (1963) and for larger Reynolds numbers the expression Dr=0.02 Udf 0.5 suggested by Hinze (1975) was used. On the pipe center where R = 0, there is an indeterminacy of type 0/0. L'Hospital's rule is used to eliminate it; the implicit dierence equation becomes:
equations were solved by means of the Thomas algorithm (Chapra and Canale 1990). The run-time of the Fortran program developed for the implementation of the algorithm on a Pentium computer, is in seconds. The solution of the dierence equations lead to nodal concentrations. The crosssectional cup mixing average concentrations were also obtained by applying Simpson's rule to equation 13 below: 1 1 2pRC dR 2pR dR 13 Ccup
Ci1,k1 ÿ 2Ci,k1 Ciÿ1,k1 Ci1,k ÿ 2Ci,k Ciÿ1,k DR2 DR2 A0 Ci1,k1 ÿ Ciÿ1,k1 Ci1,k ÿ Ciÿ1,k A1 ÿ Ci,k1 Ci,k 2R 2DR 2DR 2
A0 2
Equation 10 is valid for all nodes except for those on the pipe wall and those on the pipe axis. On the pipe wall using the boundary condition @C/ @R = ÿ A2C the dierence equation becomes:
Fv
3309
increased in the pipe. The test pipe used in the experiments is 24.3 m. (80 ft) long, and 15 cm. (6 inch) diameter unlined cast iron pipe. In order to measure continuous inlet and outlet chlorine concentrations, two Hach Cl 17 analytical devices measuring free chlorine were installed at the inlet and the outlet sections. The procedure described below was used to determine bulk decay coecients. For each run corresponding to a dierent ¯ow regime, 15 samples of 120 ml were taken both from the inlet and the outlet sections. Samples were stored in amber colored bottles and sealed with Te¯on. The amber color kept the interference of light away. The Te¯on cover was selected to avoid reaction of chlorine with other possible sealing materials. Furthermore, Te¯on is also known as hydrophobic and does not react with chlorine. The bottles were kept at constant temperature of 148C in a refrigerator. Then one bottle a day for 15 consecutive days was sampled for its free chlorine content using EPA approved DPD±FAS method. The same procedure was repeated for
Ci1,k1 ÿ 2Ci,k1 Ciÿ1,k1 Ci1,k ÿ 2Ci,k Ciÿ1,k DR2 DR2
The recursive equation 10 for intermediate points coupled with equations 11 and 12 constitute a set of simultaneous dierence equations. These set of equations were solved for various combinations of dimensionless coecients covering a wide range of geometry and ¯ow regimes. A mesh size of 100 nodes in the radial and 2000 nodes in the axial direction was used. The resulting system of
10
ÿ
A1 Ci,k1 Ci,k 2
12
dierent ¯ow regimes in the pipe. Figure 1. and Fig. 2. include data points obtained from the repeated application of this procedure to dierent ¯ow regimes for the bulk decay at the inlet and outlet sections, respectively. The following ®rst order decay rate equation (equation 14) is used to determine the bulk decay coecient k using the data given in Fig. 1. and Fig. 2. Ct C0 exp ÿkt
14
3310
Osman N. Ozdemir and A. Metin Ger Table 3. Test pipe chlorine decay for dierent ¯ows
Flow GPM (103 m3/s) 0.3 (0.02) 0.5 (0.03) 1.0 (0.06) 1.8 (0.11) 3.0 (0.19) 3.6 (0.23) 4.1 (0.26) 6.8 (0.43) 13.0 (0.82)
Reynolds number 175 292 583 1049 1749 2099 2390 3964 7579
Residence time (h) 6.54 3.92 1.96 1.09 0.65 0.55 0.48 0.29 0.15
Comparison of these ®gures revealed that while there is a single value for bulk decay constant at the inlet section, the value of the constant at the outlet section depends on the ¯ow rate. For the outlet, bulk decay constants are higher for lower ¯ow rates. The results of the bottle decay tests for the samples collected from the inlet yielded a bulk decay constant of 1.0 10ÿ3/h. On the other hand, the same tests for the samples taken from the outlet revealed that there are dierent bulk decay constants ranging from 0.8 10ÿ3/h to 5.0 10ÿ3/h depending on the rate of ¯ow in the pipe. These test results were also con®rmed with heterotrophic plate count (HPC) tests performed at the USEPA in Cincinnati. For the experiments run to obtain the bulk decay constants, extra samples were collected from the inlet and outlet sections. These samples were analyzed by HPC tests; the results are reported in Table 4. The number of colony forming units (CFU) for the samples from the outlet section were found to be very high compared to those for the inlet as the ¯ow rate decreases. These results are in agreement with the results obtained from the bulk decay tests (Fig. 1 and Fig. 2).
Inlet concentration (C0) (mg/l) 0.66 0.71 0.60 0.73 0.79 0.82 0.81 0.90 0.94
Outlet concentration (Cout) (mg/l) 0.46 0.51 0.45 0.66 0.73 0.75 0.77 0.87 0.93
Cout/C0 0.70 0.72 0.75 0.90 0.92 0.92 0.95 0.97 0.99
COMPARISON OF THE PRESENT MODEL WITH THE AVAILABLE ONES
As depicted in Figs 1 and 2 bulk decay constant k is not a ®xed value along the axis of the pipe. The assumption that the k is constant commonly employed in the previous studies is therefore questionable. To demonstrate this eect, a series of numerical runs were made. The numerical simulation depicted in Fig. 3. accentuated the signi®cance of the observation made. Dierences in pipe-end chlorine concentrations for an inlet k and variable k increases drastically with increase in L/d ratio. In Fig. 3, on y axis, Cki denotes for the dimensionless average end concentrations obtained from the application of a single inlet bulk decay constant for the whole length of the pipe whereas, Ckx denotes the one obtained from a linearly changing bulk decay from an inlet to an outlet value. Furthermore, bulk decay alone did not account for the disappearance of total chlorine in the test pipe. To be able to determine the wall decay conTable 4. Heterotrophic plate count (HPC) test results Flow GPM (lt/s)
Number of organisms inlet (CFU/ml)
0.5 (0.03) 1.0 (0.06) 2.5 (0.16) 15.0 (0.95)
6 7 4 3
outlet (CFU/ml) 240 74 23 12
Fig. 1. Bulk decay for the test pipe inlet.
Fig. 2. Bulk decay for the test pipe outlet.
Fig. 3. The eect of bulk decay change on pipe-end concentrations.
Realistic numerical simulation of chlorine decay in pipes
3311
rection factor FC, is introduced such that the 1D solution can be corrected to yield precise prediction of chlorine concentration decay along the pipe axis. The 1D chlorine decay can be represented by C1D
Ct expfÿ k Wd =r0 gt C0
15
The true decay CT is approximated by the present 2D solution, C2D. CT 1C2D Fig. 4. Comparison of the alternative models with the experimental results.
stant representing the decay of chlorine on the pipe wall a dierent series of numerical runs were made of the model. In these numerical simulations the variation of the bulk decay is assumed to be linear; for the bulk decay, the average of the inlet and outlet section was implemented. The wall decay constant value thus obtained is Wd=1.0 10ÿ6 m/s. It corresponds to minimum deviation between the experimental observations and numerical results. The models of Biswas et al. (1993) and of Rossman et al. (1994) were also compared with the experimental results for the outlet cup mixing concentrations. Figure 4 is provided to compare these results with the numerical model output presented in this study. The superiority of the proposed model over the previous ones is evident. The models considered, namely Biswas', Rossman's and the present one, are all in good agreement with the experiments for Reynolds number greater than 10,000. This is mainly because, no detectable chlorine loss can be observed for a 24 m pipe length for Reynolds number above 10,000. As the Reynolds number gets smaller, the Biswas' analytical method and Rossman's plug ¯ow approach yield large deviations from the experimental observations. In general the disagreement between Biswas' model predictions and the experiments are due to the prevailing assumption such as the one about the velocity distribution. Furthermore, the mass-transfer coecient used in Rossman's model was obtained from an experiment designed for dissolution of benzoic acid, cinnamic acid and B-naphtanol coated tube wall through a thin concentration boundary layer. This derivation requires the assumption of in®nite concentration at the wall and zero concentration at the outer edge of a very thin concentration boundary layer.
The aforementioned correction factor FC relates C1D and C2D such that FC
C1D ÿ C2D C2D
17
In other words, if FC and C1D are known, without going through the cumbersome analysis CT can be estimated by CT 1
C1D 1 FC
18
The numerical outputs of equation 15 and the present model were used to determine the functional dependence of FC on the pipe, ¯uid and ¯ow parameters. It can be shown by dimensional analysis that FC is the following form: FC fnc Z, Re
19
where the newly established dimensionless parameter Z is de®ned as Z = Wd L/n. The variation of FC thus obtained is shown in Fig. 5. Using the best ®t lines shown in Fig. 5, the data can be condensed such that FC/Z is related to Reynolds number as shown in Fig. 6. In Fig. 6 the best ®t line (equation 20) corresponding to various ranges of Reynolds numbers are also included. log FC 100=Z ÿ0:9538 log Re=1000 ÿ 1 20 Either Fig. 5 or equation 20 may be used to predict chlorine decay with ease.
THE CORRECTION FACTOR
The present model, even though produces the best ®t to the experimental observations so far, requires a signi®cant computational eort. Therefore, to ease the utility of the method a cor-
16
Fig. 5. The correction factor (FC).
3312
Osman N. Ozdemir and A. Metin Ger
improved chlorine decay simulation signi®cantly. To ease the utility of the model a correction factor is introduced to enhance the use of the simple 1D solution. AcknowledgementsÐSupport given for the preparation of this paper by the Turkish National Science Foundation is acknowledged. The Director of the Drinking Water Research Division Dr Robert M. Clark is also acknowledged for his permission to use the laboratories at EPA in Cincinnati. We also like to thank Jack Tueschler and Chuck Guion of the U.S. Environmental Protection Agency for their sincere helps during the process of collecting the data in the laboratory. Fig. 6. Compact form of equation 19. CONCLUSIONS
It is known that in the design of distribution networks the secondary pipe sizes are chosen according to the emergency ®re demands rather than the domestic consumption. This leads to oversized pipe diameters. Under normal operating conditions very slow velocities can be observed in many parts of the system. This is more severe during night time operations. Furthermore, in either hydraulic or quality solutions of distribution networks continuity has to be assumed at the junction nodes in the system. Therefore, it can readily be accepted that, pipes are connected segments of an entire system rather than individual elements. Thus, the simple one dimensional solution in some of the pipes in the network may cause considerable errors in the system. These errors accumulate at each junction leading to large errors at the end points. The present study shows that as the ¯ow regime in a pipe becomes transitional and then laminar the bulk decay rate constant along the pipe increases. The bulk decay rate increase is observed to be due to increase in number of living organisms consuming chlorine in water. It is therefore concluded that the use of a single bulk decay value for an entire network, especially the value obtained at the water treatment plant euent may not be realistic. It is suggested that measuring the bulk decay at several critical locations be made in order to de®ne site speci®c decay coecients. Furthermore, the water quality models should also be capable of incorporating the changes in bulk decay constants. The present model, incorporating variable bulk decay and using realistic velocity distribution,
REFERENCES
Bischo K. B. and Levenspiel O. (1963) Advances in Chemical Engineering, Vol. 4. Academic Press, New York. Biswas P., Lu C. S. and Clark R. M. (1993) Chlorine concentration decay in pipes. Water Res. 27(12), 1715± 1724. Buchberger S. G. and Wu L. (1995) Model for instantaneous residential water demands. J. Hydrol. Eng. 121(3), 232±246. Chapra S. C. and Canale R. P. (1990) Numerical Methods for Engineers, 2nd ed. McGraw-Hill, International Editions. Edwards D. K., Denny V. E. and Mills A. F. (1976) Transport processes. McGraw-Hill, New York. Hinze J. O. (1975) Turbulence. McGraw-Hill, New York. Hunt W. A. and Kroon J. R. (1991) Model calibration for chlorine residuals in distribution systems. Proc. Water Quality Modeling in Distribution Systems. AWWA Research Foundation, Denver, CO. LeChevallier M. W., Cawthon C. and Lee R. G. (1988) Inactivation of bio®lm bacteria. Appl. Environ. Microbiol. 54, 2492±2499. Rossman L. A., Clark R. M. and Grayman W. M. (1994) Modeling chlorine residuals in drinking-water distribution systems. J. Environ. Eng. 120(4), 803±820. Rouse H. (1949) Elementary Fluid Mechanics. John Wiley and Sons, New York. Sharp W. W., Pfeer J. and Morgan M. (1991) In situ chlorine decay rate testing. Proc. Water Quality Modeling in Distribution Systems. AWWA Research Foundation, Denver, CO. Shames I. H. (1992) Mechanics of Fluids, third edition. McGraw Hill, New York. Skelland A. H. P. (1974) Diusional Mass Transfer. John Wiley and Sons, New York. Wable O., Dumoutier N., Duguet J. P., Jarrige P. A., Gelas G. and Depierre J. F. (1991) Modeling chlorine concentration in a network and applications to Paris distribution network. Proc. Water Quality Modeling in Distribution Systems. AWWA Research Foundation, Denver, CO.