Combined use of LS-PIV and CFD for the characterization of turbulent flow in the contact chamber of ‘COSTA AZUL’ wastewater treatment plant, Carlos Paz

Water is an invaluable resource, essential for all living creatures. The continuous population growth in cities has led to the degradation of many natural ecosystems such as rivers and oceans, partly as a result of the discharge of inadequately treated urban and industrial wastewater. Therefore, from a sanitary and environmental viewpoint, it is of utmost importance to ensure the proper sanitation of urban centers and the depuration of their wastewater through sanitary works.

The requirement that the treated liquid meets quality standards in accordance with the assimilative capacity of the receiving water body encourages the development of new technologies and the optimization of existing treatment systems with the aim of improving performance in the elimination of pollutants and the sustainability of the treatment systems.

The design of this type of infrastructure is usually carried out according to sanitary engineering criteria. The hydraulic design should allow: (1) the ordinary operation of the plant when the flow is significantly less than the design flow, without causing undesirable deposition of solids in the treatment chain; and (2) that the processing units of the plant can carry the maximum flow without causing hydraulic breakage and without the flow turbulence affecting the physical processes that take place in the plant (Mays 1999).

This raises the need of accurately monitoring the operating flow rate of the plant, not only to control the treatment unit but also to characterize the flow within the unit. When operating treatment plants are evaluated, hydrodynamic singularities in the turbulent flow are often observed. These singularities can affect the quantification of the operating flow rate, also its efficiency and in consequence the final quality of treated effluent. This is true for all the treatment processes of the plant (desanders, settling tanks, contact chamber).

In recent years, the LS-PIV technique implementation has experienced an important development. This technique has been successfully implemented for the characterization of turbulent flows in large experimental laboratory facilities (e.g. in physical models, Patalano & García 2016), in river bends (Patalano et al. 2014) and initially in prototype-scale treatment plants (Ragessi et al. 2013; Ragessi 2017).

On the other hand, CFD has become a cost-effective alternative for the study of turbulent flow in water and wastewater treatment process units (Goula et al. 2008; Shahrokhi et al. 2011; Zhang 2014; Patziger & Kiss 2015; Wei et al. 2016; Ragessi et al. 2019).

In the last five years, papers have been published in an effort to present best practices that serve as guidelines (or protocols) for the use of CFD models, applied not only to improve the design but also to the daily monitoring and operation of treatment units (Laurent et al. 2014; Patziger et al. 2016; Wicklein et al. 2016).

However, the state of the art indicates that the characterization of hydrodynamic, turbulence, unstable flow and their effects on the treatment plant components remain an engineering challenge (De Clercq et al. 2002; Samstag et al. 2016). This is also valid for disinfection in contact chambers, where it is necessary to deepen the study on turbulence and/or mixing factor effects on disinfection effectiveness (Angeloudis 2014; Netshidaulu 2015).

The above indicates that the characterization of the flow turbulence in water treatment plant components remains a challenge and, furthermore, the desirability of complementing modern experimental techniques and numerical tools (CFD) for detailed flow characterization. This paper presents the case of the Costa Azul wastewater treatment plant (WWTP), located in Carlos Paz city, Córdoba, Argentina. The effluents of this plant are discharged into the reservoir of the San Roque Dam, used mainly to supply water to the city of Córdoba Capital and to the cities of Sierras Chicas (metropolitan area of Córdoba). In consequence, it is to great importance to monitor the depuration processes, including flow rate measurement, to evaluate the potential impact on the quality of the receiving water body and the preservation of the public health.

Currently, the WWTP lacks an accurate flow monitoring system at the treatment outlet. Near the contact chamber outlet section, an ultrasonic sensor VEGASON 61 is installed for continuous level measurement (see location on Figure 1). Unfortunately, an electrical failure put it out of operation.

This problem led to the development of a methodology to calibrate a discharge-measurement system at the outlet of the WWTP contact chamber, prior to discharge into San Roque Lake.

In a first stage of this work, large-scale particle image velocimetry (LS-PIV) and computational fluid dynamics (CFD) were applied in combination, both at a prototype scale in the contact chamber. The details of LS-PIV and CFD model used are described in the ‘Materials and methods’ section.

Experimental results and numerical simulations made it possible to determinate two operating flow rates, build a head-flow curve, obtain a detail characterization of flow hydrodynamics and estimate characteristic parameters of the contact chamber.

Finally, in a second stage, five design alternatives were modeled to improve the flow hydraulics in the curved area of the contact chamber, with the objective of minimizing stagnation zones, recirculation and shortcuts developed by the flow, in order to achieve the most uniform possible distribution of residence time in the chamber.

All experimental and numerical results achieved are presented in the ‘Results’ section. Finally, the principal conclusion of this work is summarized.

The LS-PIV technique was applied in two zones located upstream of the discharge weir, at the outlet of the contact chamber (see Figure 1). In zone 1 (180° curve), the velocity field was determined to evaluate separation zones, recirculation and low velocity zones. In zone 2 (straight section), the surface velocity field was used to estimate the operating flow rate.

On two occasions, videos were recorded with a resolution of 1,280 × 720 pixels. Wood chips were used as tracer. Among the reasons for this choice is that they have little inertia, are dragged by the flow almost without disturbing it, provide a good contrast on the water surface and are biodegradable. The size of the wood chips varies from the order of 2 or 3 mm to 1 cm.

First, images were extracted from the video, with a constant time interval between them of 33.32 ms (30 fps), for a video length of 1800s. Then, the computational package PIVlab (Thielicke & Stamhuis 2014), which is included in RIVeR (Patalano et al. 2017), was used to process these images. For this purpose, the discrete fast fourier transformation (DFFT) algorithm was used, using double pass, with square interrogation window sizes of 64 px for the first pass and 32 px for the second pass, using Windows deformation between them, with a 50% overlap in both cases.

Second, the mean velocity fields were estimated and the results obtained were rectified. For this purpose, four control points were defined in the channel located in the same plane as the free surface of the flow, and the distances between them were taken: see for example the control points in Figure 1. Knowing these distances, and the pixel coordinates of these points in an image, it is possible to remove the perspective distortion of the results from the resolution of the two-dimensional homographic matrix.

Finally, the flow rate was estimated from the mean surface velocity profile and the surveyed section (width and depth).

The ratio depends on many parameters, such as the geometry section, the bed roughness, the geometry upstream and downstream of the section, secondary flows, and others. It is up to the user to choose the most reliable value of ⁠. Cheng et al. (2004) concluded that it is reliable to use the surface velocity as an index to determine a river flow discharge because the value of always falls in the same range, from 0.80 to 0.93.

In this case, the channel has a uniform rectangular section but the measured flow is downstream of a curve. Therefore, a value for the ratio between the mean flow velocity and the surface velocity was adjusted based on the CFD model results (see Model 1). In this case, was equal to 0.82. This is explained in detail in the ‘Flow velocities’ section.

The estimated flow rates during the two in-situ measurements (made with LS-PIV) did not change significantly. At the moment, the WWTP is operating well below its flow capacity because the sewer network only serves a smaller percentage of the city. For this reason, it was decided to use a numerical model to simulate a wide range of flows and to calibrate the weir height-discharge curve at the outlet of the chamber.

Two numerical schemes were selected using the free and open package OpenFOAM^® (Open Field Operation and Manipulation, OpenCFD Ltd 2007).

First, a section of the flume and the discharge weir were modeled to determine the weir discharge curve (see Figure 2).

Here is the mean i component of velocity vector (i = 1–3), is mean pressure, is water density, is momentum diffusivity (viscosity), k is turbulent kinetic energy, is turbulent viscosity, is rate of dissipation of k, is Dirac delta function, is ratio of momentum diffusivity (viscosity) and k is diffusivity, is the ratio of momentum diffusivity (viscosity) and is diffusivity, ⁠, ⁠, are − model parameters. In this work, ⁠, ⁠, ⁠, and based on Launder & Spalding (1972, 1974) recommendations as it is implemented in OpenFOAM^®. These parameters have not been modified in this work.

The scheme PIMPLE solver (pimpleFoam in OpenFOAM^®) was used for numerical integration. PIMPLE merges the PISO (Issa 1986) and SIMPLE (Patankar 1980).

The mesh was generated with Gmsh^® (Geuzaine & Remacle 2009) and was discretized with cells of variable size, from 0.5 to 20 cm on a side, prioritizing the highest resolution in the area near the walls, in the water-air surface and in the outlet weir. A mesh convergence analysis was performed in order to select a mesh that offered acceptable results with a relatively low computational cost.

The boundary conditions (BCs) imposed on the numerical domain were at the flume inlet and outlet, at the free surface, and the walls (flume and weir). To match the experimental conditions, a variableHeightFlowRateInletVelocity boundary was applied to U at the inlet boundary, which adjusted the flow depth and velocity at the inlet to obtain a constant flow rate (⁠⁠) and a variableHeightFlowRate boundary condition was used to enable the water level to oscillate. The outlet boundary was controlled by a zero gradient (Neuman BC) velocity condition, and no-slip boundary conditions were imposed on all walls. The top surface of the mesh was considered a free surface, and a pressureInletOutletVelocity condition for U was imposed to allow the flow to enter and leave the simulation domain freely. Regarding the variables in the RANS model, the values of k, ⁠, and at the inlet and outlet boundaries were set to arbitrary low values and allowed to develop within the simulation domain. The standard wall function was employed so that the governing equations could be solved correctly.

Second, the inflow and a contact chamber portion were simulated as shown in Figure 3. A flow rate of 100 L/s was modeled, equal to that observed in the second in-situ measurement, to estimate the residence time using tracer particles in the resolved flow.

In this case, also a RANS – model (see Equations (1)–(5)) and PIMPLE scheme were used, but to numerically solve a single-phase flow field. This hypothesis was assumed because in prototype, the water surface does not show level fluctuations for a constant flow.

For this reason, in this model a rigid surface was imposed with the slip boundary condition, which does not alter the tangential velocity component on the top surface. In this way it is possible to use a one-phase (water) model. This simplification made it possible to reduce the computational effort required. The height of the water column, for the flow rate modeled, was defined from the height-flow curve calibrated with the interFoam model results.

The other boundary conditions (BCs) imposed were a constant uniform value for U at the inlet to obtain a constant flow rate of 100 L/s, the outlet boundary was controlled by a zero gradient (Neuman BC) velocity condition and no-slip boundary conditions were imposed on the side and bottom walls. The values of k, ⁠, and at the inlet and outlet boundaries were set to arbitrary low values and allowed to develop within the simulation domain. The same standard wall function as in Model 1 was employed. For p a zero gradient (Neuman BC) was used at the inlet and equal to zero at the outlet.

The mesh was completely generated by Gmsh^® and discretized with cells of variable size, from 0.5 to 10 cm on a side. The maximum resolution was used near the side and bottom wall.

In the post-processing stage, tracer particle propagation analysis was incorporated in order to evaluate more accurately the residence time in a portion of the contact chamber. To do this, a particle tracer filter from ParaView (Ahrens et al. 2005) was used in non-stationary flow. The post-processing was made on 200 seconds of simulation results.

Third and finally, using the same pimpleFoam scheme, a two-dimensional model was used to simulate the flow for five design alternatives. The objective of each design is to improve hydraulics in the curves, with the criteria of reducing dead flow zones, recirculation and short circuits that could generate a wide range of residence times.

The boundary conditions (BCs) were set similar to Model 1. The only difference is that empty conditions were used in the top and bottom surface so the vertical direction is not solved.

The mesh was generated by combining Gmsh^® and the snappyHexMesh of openFOAM^®. The resolution is similar to Model 1, but the mesh has only one cell in the vertical direction.

The design alternatives are presented in Figure 4. In Scheme 0, the apex of the bulkhead was extended 1.0 m into the curve. In Schemes 1 and 2, deflecting bulkheads are incorporated to orient the flow in the curve area. The design of Scheme 3 modifies the shapes of the apex in the curve and incorporates chamfering of the outer area of the curve. Last, the design of Scheme 4 modifies the apex of the curve incorporating an airfoil.

In the 180° curve, the experimental results show a strong flow separation and the development of a recirculation zone with low velocities downstream of the apex (see blue color, in Figure 5(a)). The recirculation zone extends beyond 9 m downstream and uniform flow cannot be developed. Also, in the outer zone of the curve (see white rectangles in Figure 5(a)), very low velocities are observed, forming stagnation zones.

In Figure 5(b), the streamlines show a recirculation zone immediately downstream of the curve as a consequence of flow separation. The effect of recirculation gradually decreases along the reach, although it does not complete disappear and its impact on the transverse velocity profile extends over an extended length. It is noteworthy that when measured 4 m downstream from the curve, it was detected that the velocity profile is still affected by the separation.

The behavior of the surface velocity profile was first analyzed in four sections downstream of the curve. The first section is located 1.0 m downstream of the curve and the distance between sections is 1.0 m. Sections 1–4 are shown in Figure 5(c) and the surface velocity profiles are presented in Figure 6(b).

Negative velocity components (opposite to the main flow direction) can be seen in most of the profile in Sections 1 and 2. The effect of recirculation progressively decreases downstream and in Section 4, the furthest from the septum apex (4.0 m), the velocities are positive throughout the profile.

These results are from the second measurement with LS-PIV. However, the results of the first in-situ measurement show a similar flow behavior.

In Figure 6(b), the two remaining profiles (Sections 5a and 5b) are located 7 m downstream from the curve. From this surface velocity profiles flow rate was estimated. Sections 5a and 5b profiles are the results of the first and second LS-PIV measurements, respectively.

In the first measurement flow depth was 1.48 m and flow rate Q was equal to 0.12 m³/s and in the second measurement flow depth was 1.485 m and Q was equal to 0.11 m³/s.

In both cases was used to calculate Q. The value was calculated as the ratio of the depth-average velocity to surface velocity: u at 0.6 of water depth and at the top surface, both from Model 1 results. This was carried out for two flow rates: 100 and 140 L/s.

Figure 6(b) shows the dimensionless velocity profile. In this case, was equal to 1.223, so was equal to 0.82. This value is among those expected according to the literature.

The weir discharge curve (⁠) was determined from the Model 1 results (see Figure 2). Figure 7 shows the curve for a range of flow rates from 50 to 500 L/s and compares it with the Q estimated from LS-PIV measurements (points labeled LS-PIV in Figure 7). The experimental results include a 10% error bar.

The curve fitted and experimental data show a good agreement, with differences less than 10% for low flow rates. It was possible to fit an within the recommended values (0.8–0.95), and the LS-PIV technique proved to be a cost-effective alternative to estimate two operating flow rates.

Another source of error in estimated Q could be due to the precision to which the water level H must be measured in the flume. It should be noted that the fitted curve is very sensitive to the H values: for example, a 1 cm variation in H produces a variation of 17 L/s. This variation can represent 34% for smaller flows (of the order of 50 L/s). In other words, a great deal of precision in the measurement of H is necessary to minimize the uncertainty in the determination of Q. This is due to the fact that the weir is 2 m long, approximately twice the width of the channel. In addition, the weir is located lateral to the main flow direction.

The experimental results and the numerical model so far presented indicate that to minimize uncertainty in the determination of the flow rate, it is advisable to change the location of the weir (transverse to the main flow direction), and minimize its length or change to a triangular type weir.

In this case, L was determined from the flow lines reported by the LS-PIV results in the curve and a straight section. The results were extrapolated to the entire volume of the chamber to determine the average flow trajectory. L is 319.4 m long.

The values of the contact chamber parameters are presented in Table 2 for the two measurements performed together with results from the CFD model (see ‘Analysis of residence time using CFD model’ section).

First and second LS-PIV measurement results show approximately a 30 minute difference between and ⁠. Future work will evaluate this difference for higher flows, since when measurements were made it was not possible to modify the flow discharge in the WWTP contact chamber.

The WWTP operators indicated that the weir crest level was raised from the original level in order to increase the total volume of the contact chamber and the ⁠, and to ensure the disinfection. The parameters reported in Table 1 indicate that the control structure (weir) at the outlet of the contact chamber configures a flow of low F values (less than 0.1), which facilities the propagation of instabilities along the channel and hinders the development of a uniform flow. As a result, strong recirculation and stagnation zones develop, as well as short circuits (shortcuts) that affect the mixing process. The existence of low and high velocity zones results in a wide range of residence times, with minimum values that can be well below the mean residence time (close to ⁠).

The residence time dispersion could also be observed from the analysis of the propagation of tracer particles in different sections of the contact chamber by a CFD model. Figure 8 presents the result of the propagation of a set of 50 tracer particles for 100, 150 and 200 seconds of simulation. All the particles were seeded at the same time in the input section of the model (see Figure 3).

The trajectory was analyzed for groups of five particles close to each other. The mean trajectory (⁠⁠) of each group was measured over 200 seconds. Then, group mean velocity (⁠⁠) was calculated as ⁠, where seconds. Finally, group residence time (⁠⁠) was calculated as ⁠, where L is the average chamber trajectory.

Table 1 presents the results of ⁠, and for each group of particles.

In Table 2, contact chamber parameters are presented for the two LS-PIV measurements performed, together with results from Model 2. In the Model 2 column, u and result from the average of and ⁠, respectively.

LS-PIV results show a 2.1 times greater than the numerical model. This is because Model 2 allows analyzing the flow behavior (using tracer particles) in a larger volume of the contact chamber than the LS-PIV measurement area.

The results in Table 1 show that most of the tracer particles fulfill the recommended minimum residence time. In Argentina, the standards recommend a minimum time of 20 minutes for peak flow (National Entity of Public Sanitation Works – ENHOSA, Argentina), while some international standards recommend a minimum of 20 minutes (US-EPA 1986).

However, the numerical model results also indicate that the tracer particles describe very different trajectories. It should be noted that some particles left the chamber domain before 20 minutes and the results correspond to a flow rate of 100 L/s, well below the average (250 L/s) and maximum (500 L/s) capacity of WWTP. In future work, the minimum residence time for medium and high flows will be evaluated. These results should be taken into account for chlorine dosing.

Finally, in an attempt to reduce the recirculation zones, stagnation and shortcuts developed by the flow in the chamber, five alternatives were evaluated to improve the hydraulic design. The schemes are presented in Figure 4.

Next, the results of velocity and flow lines in each of the alternatives are presented in Figure 9.

Table 3 presents four parameters to compare flow in five schemes: the recirculation length and width (⁠ and ⁠) downstream of the curve; low velocity zone width (⁠⁠) on the outside of the curve; and the maximum and mean velocity ratio (⁠⁠) in the section.

The extension of the bulkhead (Scheme 0) reduces the volume of recirculation and stagnation of flow on the outside of the curve, but without significant differences with the current situation, while the design that combines extending the bulkhead by 1 m and incorporating two ‘I’ and ‘L’ deflectors (Scheme 1) achieved a significant reduction in the recirculation zone downstream of the curve, in length and width. The low velocity zone on the outside of the curve was significantly reduced. The bulkhead in ‘I’ form (Scheme 2) does not allow reduction of the recirculation zone downstream of the curve and the low-velocity zone. The design of Scheme 3, with a chamfered profile at the apex of the curve and 45° corner trimming, reduces the recirculation zone downstream of the curve, but a low-speed zone continues to develop on the outside of the curve. The curve airfoil (Figure 4) did not prevent the flow from separating and the development of instabilities (eddies, vortices) that travel downstream, preventing the development of a uniform flow in the straight section of the flume.

In summary, Scheme 1, which contemplates extending the bulkhead 1 m inward from the curve and incorporating two baffles (one L-shaped), drastically reduces the recirculation and stagnation zone of the flow: a, b, c and have the lowest values in combination. This design could reduce dispersion in residence time distribution.

In the design of chlorine contact tanks, the principle concern must be to achieve maximum disinfection efficiency with a minimum chlorine residual. Therefore, the numerical simulation results suggest that the hydraulic design of the contact chamber should consider two aspects: first, minimize flow shortcuts to ensure sufficient and uniform contact time to inactive pathogens and, second, avoid disinfectant concentration due to prolonged residence time in dead zones, which may promote the formation of disinfection by-products.

The experimental results and the numerical simulation allowed the obtaining of a detailed characterization of the flow hydrodynamics, the determination of flow rates and the calibration of the head-flow curve at the treatment outlet.

The velocity field obtained with the LS-PIV technique reports that the flow recirculation and the stagnation zones affect a large percentage of the contact chamber volume. Consequently, the average flow trajectory is much smaller than theorical and the effective residence time is much shorter than estimated based on the mean flow rate. It is important to evaluate this for higher flow rates.

The numerical simulation allowed reconstructing the height-flow curve of the weir for a wide range of flow rates. Based on the results, it is recommended to change the location of the weir (transverse to the main flow direction), and minimize its length or change to a triangular type weir, to minimize the uncertainty in the determination of the flow rate.

The results of the numerical model indicate that a set of tracer particles describes very different trajectories. The numerical calculation indicates that most of the tracer particles meet the minimum residence time of 20 minutes (indicated by ENHOSA standards), but some particles left the chamber domain before 20 minutes. This should be taken into account for chlorine dosing.

The LS-PIV technique, with low implementation and calibration costs, proved to be a cost-effective alternative to estimate two operating flow rates and characteristic parameters of the contact chamber. This alternative, complemented with the use of CFD-type models, allowed more detailed analysis of the flow behavior and residence time to be carried out throughout the contact chamber.

LS-PIV results show a 2.1 times greater than the numerical model. This is because Model 2 allows analyzing of the flow behavior (using tracer particles) in a larger volume of the contact chamber than the LS-PIV measurement area. This difference suggests that a full-scale experiment should be performed with LS-PIV to more accurately estimate residence time.

In particular, the results of the numerical simulation suggest that the Scheme 1 allows a major reduction of the recirculation and low-velocity zones. The main benefit of Scheme 1 is that it would result in a decreased dispersion in the residence time distribution.

In general, the numerical simulations performed indicate that the hydraulic design of the contact chamber should consider two aspects: first, minimize flow shortcuts to ensure sufficient and uniform contact time to inactive pathogens and, second, avoid disinfectant concentration due to prolonged residence time in dead zones, which may promote the formation of disinfection by-products.

All relevant data are included in the paper or its Supplementary Information.