*Article* **Validation of a Computational Fluid Dynamics Model for a Novel Residence Time Distribution Analysis in Mixing at Cross-Junctions**

#### **Daniel Hernández-Cervantes 1, Xitlali Delgado-Galván 2, José L. Nava 2, P. Amparo López-Jiménez 3, Mario Rosales <sup>1</sup> and Jesús Mora Rodríguez 2,\***


Received: 30 April 2018; Accepted: 28 May 2018; Published: 5 June 2018

**Abstract:** In Water Distribution Networks, the chlorine control is feasible with the use of water quality simulation codes. EPANET is a broad domain software and several commercial computer software packages base their models on its methodology. However, EPANET assumes that the solute mixing at cross-junctions is "complete and instantaneous". Several authors have questioned this model. In this paper, experimental tests are developed while using Copper Sulphate as tracer at different operating conditions, like those of real water distribution networks, in order to obtain the Residence Time Distribution and its behavior in the mixing as a novel analysis for the cross-junctions. Validation tests are developed in Computational Fluid Dynamics, following the *k-ε* turbulence model. It is verified that the mixing phenomenon is dominated by convection, analyzing variation of Turbulent Schmidt Number vs. experimental tests. Having more accurate mixing models will improve the water quality simulations to have an appropriate control for chlorine and possible contaminants in water distribution networks.

**Keywords:** water distribution networks; EPANET; safe water

#### **1. Introduction**

In recent years, useful life of a several number of Water Distribution Networks (WDN) has been exceeded, so that there are certain points or sectors that do not reach enough disinfectant. This problem could produce the formation of biofilm on pipe walls [1,2], promoting bacterial growth and its transport through water. On the other hand, the system's own useful life, the intermittency in the supply, the leaks level, and the demands behaviour cause a great variety of residence times of the drinking water and variations of reaction of the disinfectant in the network [3,4].

Chlorine disinfection is one of the principal factors in drinking water treatment processes. This chemical agent is used mainly to ensure the destruction of pathogen organisms that could be present in water and that are coming from the sources of supply. The disinfection prevents gastrointestinal diseases due to the consumption of contaminated water. However, it is also true that chemical disinfection has caused unwanted risks due to by-products that are generated in the reactions with contaminants present in water, such as trihalomethanes (THM), among others [1,2,5–9].

Finding solutions to the problem of inadequate use of chlorine is an issue that involves various aspects. When WDN presents a high degree of reaction due to pipe conditions and residence times, concentrations in some areas decrease rapidly [2,3,7]. Due to this, the absence of the necessary amount of disinfectant induces the possible consumption of water that is contaminated by users. A commonly used solution is the increased chlorine concentration from the sources. However, demand nodes near the sources distribute water with the high chlorine contents and due to factors, such as residence times, the formation of harmful by-products increases, which in the same way affects the consumer population [9,10].

Location of strategic points to control the problem of some solute distribution can be reached through the applications, in which it is necessary to improve the incomplete mixing model. Diverse authors [11–14] have developed models based on heuristic technics that use the results of water quality analysis to determine the best location of booster chlorination stations to optimize their usage. Other authors [15–20] developed researches that determine the strategic zones for the sensors location to detect pollutants, with the aim of finding accidental or intentional action of harmful agents and with the idea that the consumer has the risk of consuming contaminated water. Both kinds of applications predict the best zones for the water quality control devices; however, they are dependent on the way, in which EPANET (is a computer program that performs extended period simulation of hydraulic and water quality behavior within pressurized pipe networks [21]. The name comes from the U.S. Environmental Protection Agency and Networks) simulates the distribution of a solution through the network; the predictions are not accurate as many cross-junctions are presented in the networks that are studied.

An incomplete mixing analysis is presented in experimental tests that use injection pulses of copper sulphate (CuSO4) as a tracer. Residence Time Distribution (RTD) curves are good features and a novel proposal on mixing cross-junctions, in order to calibrate simulation processes, with the help of stimulus-response techniques [22–24]. There were generated turbulent regime scenarios in Reynolds numbers from 43,343 to 155,793. Therefore, different operating conditions that are close to real WDN were performed. To analyse the tracer distribution in the outlets of the cross-junction, the curves of the pulses are registered, visualizing by RTD curves the approximation of the Computational Fluid Dynamics (CFD) simulations. According to the cross-junction mixing models that are reviewed in the literature, this is the first time that the hydraulic pipe flow has been simulated with success using RTD analysis, which was addressed by solving the RANS equations with the k-epsilon turbulence model coupled to the diffusion-convection equation.

#### *Complete Mixing Model*

In EPANET and in most hydraulic simulation programs, it is assumed that the solute mixing is complete and instantaneous [21]. Therefore, the concentration of the substance (for example: residual chlorine) that leaves the cross-junction, will be the weighted average of the solute concentrations of the arriving pipes. For a given node *k*, the concentration leaving the node is obtained by Equation (1).

$$\mathbf{C}\_{i} = \frac{\sum\_{j \in I\_{k}} Q\_{j} \mathbf{C}\_{j} + Q\_{k, \text{ext}} \mathbf{C}\_{k, \text{ext}}}{\sum\_{j \in I\_{k}} Q\_{j} + Q\_{k, \text{ext}}} \tag{1}$$

where: *i*, link with flow leaving node *k*; *Ik*, set of links with flow into *k*; *Lj*, length of link; *Qj*, flow (volume/time) in link *j*; *Qk,ext*, external source flow entering the network at node *k*; and, *Ck, ext*, concentration of the external source flow entering at node *k*.

Several authors [25–36] have numerically and experimentally demonstrated that the mixing in these unions is far from being "complete and instantaneous". Most of these researches were based on physical experiments and simulations through Computational Fluid Dynamics (CFD) using tracers for the concentration distribution analysis. Some researchers [26–28,31–36] implement experimental scale models of a cross-junction and tee junctions with different pipe diameters, which are driven by elevated tanks where tracer is injected into one of them. Once the mix of the flows was made, mixing coefficients were formulated based on their experimental measurements. Some others based their approaches on CFD simulations directly, even for small geometric configurations [25,29,30]. The most

commonly used tracers were salts, such as NaCl. Other authors considered colored dyes and high speed photography too [34,35]. Conducting a study directly with chlorine as a solute turns out to be somewhat complicated, although Mompremier [33] used in-line chlorine loggers/dosers to analyze the mixing at different conditions.

The developed studies carried out to date allow for the derivation of mathematical functions for the computation of the incomplete mixing, even two codes have been formulated to be implemented as toolkits for EPANET (one of them, called AZRED II and developed by [36]). However, these applications have not been disseminated globally yet, because most of them were performed by scaled experimental devices reaching low velocities and turbulent flows that were close to 4000. In some cases, they are still subject to specific calibrations of each network and there is a lack of experimental sufficiency in real networks of high interconnection.

The lack of precision in water quality simulations can lead to erroneous projects of monitoring systems of the location of re-chlorination systems necessary for the control of disinfection [28]. The accuracy in these models is also necessary to simulate the spatial-temporal dispersion of chemical and microbial agents during accidental or intentional contamination events.

This paper presents a validation of the phenomenon of mixing in cross-junctions by means of CFD, in which the distribution of the disinfectant is not considered as homogenous mixing. The presented scenarios were performed at velocity and pressure conditions that were similar to those real WDN thanks to the impulsion and conduction equipment used. This avoids making measurement adjusts or extrapolations, that makes it different from many other jobs working at lower turbulence flows. The model is validated by a novel analysis of RTD curves for mixing cross-junction, resulting from the execution of tests of different velocity conditions and its corresponding measurement of copper sulfate injection response. The correct model meshing degree is verified for the study of the phenomenon in the CFD model.

#### **2. Materials and Methods**

#### *2.1. Experimental Model*

The cross-junction used for the mixing analysis is situated on an experimental network of the "Universidad de Guanajuato"; the pipe material is galvanized iron of four inches in diameter (Figure 1). The flows are conducted by a 15HP hydraulic pump, which reaches 30 L per second at pressures of 1.78 bar for the pipes on the cross-junction. Open valves are available at each border for flow velocity control for inlets and outlets. They are located at 2 and 7.5 m upstream water from the cross-junction. It is done to reduce turbulence generated by the obstruction of the valve, and so that the flow could be uniformized along the pipeline.

**Figure 1.** Cross-junction boundaries. The distances correspond to the pressure registration points.

Operating conditions that are close to WDN can be approximated by the physical characteristics of the experimental network. All of the mixing tests were developed contemplating two inlets and two outlets flow. The boundary distances (shown in Figure 1), that were established for the numerical model, are the location points of the pressure measuring devices. The north inlet (N) corresponds to the upper boundary. The West inlet (W) is the left boundary, and in the same way for the outlets boundaries on the East (E) and South (S) (Figure 1). In addition, the numerical model is performed with the objective of maintaining the flow direction with a sense from left to right and from top to bottom.

North and West inlets flows are measured with the support of two flow meters, of them being of electromagnetic type, DOROT brand, model MC608-A (Figure 2a). The second one is a propeller meter, of the company Hidrónica®, model MP-0400 (Figure 2b), both with datalogger capabilities. For measures at the outlets, the outflowing water is conducted up to four flowmeters, two for each outlet. For the East outlet, the flowmeters are of electromagnetic type, of Badger Meter Inc. brand, model M2000 (Figure 2c). For the South outlet, the water is conducted to two propeller volumetric flow meters (Figure 2d).

**Figure 2.** Equipment to measure flow (electromagnetic: (**a**,**c**), propeller volumetric: (**b**,**d**)) and pressure around the cross-junction (oscilloscope and pressure transducer: (**e**)).

A digital storage oscilloscope with four measuring channels, a Tektronix® brand model TDS-2004C (Figure 2e) was available to the pressure registration on the four boundaries of the cross-junction. In this case, the oscilloscope visually represents the electrical signals that are captured by pressure transducers. These transducers are of WIKA® brand model S-10, comprise a measuring range of 0–10 V of electrical signal intensity corresponding to the range −1 to 9 bar of pressure.

The tracer that was used for the tests is copper sulphate pentahydrate (CuSO4·5H2O) of fine grade, Fermont® brand. The 0.25 mol/L solutions are prepared, which return an effective response to the concentration measurement. The tracer injection is done by pumping, with a monophasic equipment, Siemens brand, 1HP power and rotor speed of 3515–3535 rpm. The tracer is stored in a volumetric test tube of four liters size, and then, in each test, the tracer is pumped 150 mL approximately. The injection point is located at 3.50 m before the cross-junction (Figure 3). The pump equipment allows for overcoming the pressures of water flow in the network and the controlled intrusion of tracer that is injected through the opening and closing ball valve of a 0.5 inch.

A potentiostat-galvanostat, brand BIOLOGIC, model SP-150 was used to apply an electric potential, with which copper adhesion will occur on the electrodes (Figure 4). This equipment has software that allows for current intensity registration at fixed time intervals. The electrodes were placed at the outlet boundaries of the pipes, close to the pressure transducers. The electrodes must be made of a metal with suitable properties for electrical conduction. They were developed with the characteristics that are shown in the Table 1.

**Figure 3.** Tracer injection and registration points.

**Figure 4.** Equipment for stimulus-response technique.


The velocity alternatives for the case of two inlets and two outlets are described in Figure 5. To simulate equal velocities is difficult because of working with high magnitudes of velocity and pressure. Therefore, the scenarios were approximated as much as possible to cover the four-velocity combinations in Figure 5. It has been determined that the flow velocities are closer when the variation between both of the inlets or outlets boundaries are by a percentage with a relative difference is less than 25%, otherwise, for the variations that are presented with a relative difference of 25% or greater, it is considered that the velocities are different. The objective is to cover trials that involve the different mixing conditions, likely in a cross joint. The alternatives are classified in these four groups, according to their velocities. Table 2 contains contents values of the four scenarios V1 to V4 proposed.

*Water* **2018**, *10*, 733

$$\begin{array}{ccccc} \Sigma\_{\mathsf{A}} & \mathsf{V}\_{\mathsf{A}} = \mathsf{V}\_{\mathsf{W}} & \mathsf{C}\_{\mathsf{A}} \neq \mathsf{V}\_{\mathsf{W}} \\\\ \Sigma\_{\mathsf{A}} & \mathsf{V}\_{\mathsf{E}} = \mathsf{V}\_{\mathsf{S}} & \mathsf{C} & \mathsf{V}\_{\mathsf{E}} = \mathsf{V}\_{\mathsf{S}} \\\\ \Sigma\_{\mathsf{E}} & \mathsf{V}\_{\mathsf{K}} = \mathsf{V}\_{\mathsf{W}} & \mathsf{C}\_{\mathsf{A}} & \mathsf{V}\_{\mathsf{K}} \neq \mathsf{V}\_{\mathsf{W}} \\\\ \end{array}$$

**Figure 5.** Boundary velocities for experiments.

**Table 2.** Contents values of the four scenarios V1 to V4 proposed.


#### *2.2. Formulation of Numerical Simulation*

#### 2.2.1. Turbulent Flow

The numerical simulations for the dynamic flow are based on the mass conservation Equation (2) and the Reynolds Averaged Navier-Stokes (RANS). These last equations consider the turbulent flow applying Reynolds decomposition by its average over a small-time increment; the flow quantities are decomposed in a temporal mean and a fluctuating component [37]. The decomposition results in the RANS Equation (3) [38,39].

$$\nabla \cdot (\rho u) = 0 \tag{2}$$

$$\rho \left( \mu \cdot \nabla \right) \mu = -\nabla P + \nabla \cdot \left( \left( \mu + \mu\_T \right) \left( \nabla \cdot \mu + \left( \nabla \cdot \boldsymbol{u} \right)^T \right) \right) \tag{3}$$

where *u* is the mean-averaged velocity vector, *P* is the pressure, *μ* is the dynamic viscosity, and *ρ* is the fluid density.

The fluctuating component introduces a closure problem [37], which could be solved by diverse type of models: Eddy Viscosity or First-Order and Second-Moment Closure (SMC). The standard *k-ε* model is one of SMC that works with high Reynolds number turbulence valid only for fully turbulent free shear flows that cannot be integrated all the way to the wall [40]. Modeling flows close to solid walls requires integration of the two equations over a fine grid to correctly capture the turbulent quantities inside the boundary layer [40,41]. Then, the so-called Reynolds stresses can be stated in terms of a turbulent viscosity *μT*, according to the standard *k-ε* turbulence model; Equations (4) to (6).

$$
\mu\_T = \rho \mathbf{C}\_\nu \frac{k^2}{\varepsilon} \tag{4}
$$

$$\rho(\mu \nabla)k = \nabla \cdot \left( \left( \mu + \frac{\mu\_T}{\sigma\_k} \right) \nabla k \right) + P\_k - \rho \varepsilon \tag{5}$$

$$
\rho \mu \cdot \nabla \varepsilon = \nabla \cdot \left( \left( \mu + \frac{\mu\_T}{\sigma\_\varepsilon} \right) \nabla \varepsilon \right) + \mathbb{C}\_{\varepsilon 1} \frac{\varepsilon}{\kappa} P\_k - \mathbb{C}\_{\varepsilon 2} \rho \frac{\varepsilon^2}{k} \tag{6}
$$

where *k* is the turbulent kinetic energy, *ε* is the turbulent energy dissipation rate, *Pk* is the energy production term (*Pk* = *μT*[∇*u* : (∇*u* + (∇*u*) *<sup>T</sup>*]), and *<sup>C</sup><sup>μ</sup>* (0.09), *Ce1* (1.44)*, Ce2* (1.92), *<sup>σ</sup><sup>k</sup>* (1), and *σε* (1.3) are dimensionless constant values that are obtained by data fitting over a wide range of turbulent flows [38]. K. Hanjalik [41] used the *k-ε* model over other types of turbulence models for simulation in pressure pipes with high Reynolds numbers, with the same empirical constants. Furthermore, RANS *k-ε* models have proposed good results also in mixing and concentration analysis, when these models have been validated with measurements, as in [42,43].

#### 2.2.2. Boundary Conditions

Due to the high values of the Reynolds number that were reached in the scenarios, velocities in regions near the walls model decrease rapidly. The fluid velocity becomes dependent on these boundaries. The wall functions used are based on a universal distribution of velocities depending on the proximity of the wall, Equation (7).

$$
\mu^{+} = \frac{1}{\kappa} \ln \left( E y^{+} \right) \tag{7}
$$

where *u*<sup>+</sup> is the normalized velocity component in the logarithmic layer of the wall, *κ* is the Von Karman constant (0.4187), *E* is a constant that depends on the roughness of the wall, and *y*<sup>+</sup> is the dimensionless distance from the wall.

The input parameters that are assigned to the CFD analysis (shown in Table 2) were established them in the following way, at the same distance like in the experimental model:

Inlet velocity at N; Inlet velocity at W; Pressure at E outlet; and, Pressure at S outlet.

The velocity values (Table 2) are multiplied by the unit normal vector n, *u* = −*v*n, then, the input velocity is a vector of equal magnitude in the surface input boundary. The parameters *k*<sup>0</sup> and *ε*<sup>0</sup> must be defined, which are obtained from the Equations (8) and (9).

$$k\_0 = \frac{3}{2} (v\_0 I\_T)^2 \tag{8}$$

$$
\varepsilon\_0 = \frac{C\_\mu^{\frac{3}{4}} k^{\frac{3}{2}}}{L\_T} \tag{9}
$$

Turbulent intensity *IT* (dimensionless) varies between 0.05 and 0.10. For the turbulent scale length, *LT* (long unit) can be obtained according to the pipe radius *r*, obtaining 7% of it. *LT* = 0.07*r* [39]. Therefore, the values were fixed on *IT* = 0.05 and *LT* = 0.07 × 0.0508 m = 0.003556 m.

#### 2.2.3. Residence Time Distribution

To sketch RTD curve in a fluid for transport of diluted species in dynamic regime, an averaged convection-diffusion equation for turbulent flow is used, Equation (10).

$$\frac{\partial \mathbf{C}}{\partial t} = -\mathbf{u} \cdot \nabla \mathbf{C} + \nabla \cdot (D + D\_T) \nabla \mathbf{C} \tag{10}$$

where *u* is the velocity vector determined from Equation (2) of the hydrodynamic analysis, *C* is the tracer concentration, and *D* is the diffusion coefficient of the tracer. Under turbulent flow conditions, the concentration and the fluctuation parameters can be grouped in a turbulent diffusion term, *DT*. This term can be evaluated considering the Schmidt Number (turbulent), *ScT*, for which, the equation proposed by Kays-Crawford [41] was used, Equations (11) and (12).

$$\mathcal{S}c\_T = \frac{\mu\_T}{(\rho D\_T)}\tag{11}$$

$$Sc\_{T} = \frac{1}{\frac{1}{2\text{Sc}\_{T\infty}} + \frac{0.3\mu\_{T}}{\sqrt{\text{Sc}\_{T\infty}\rho D}} - \left(0.3\frac{\mu\_{T}}{\rho D}\right)^{2} \left(1 - \exp\left(-\frac{\rho D}{0.3\mu\_{T}\sqrt{\text{Sc}\_{T\infty}}}\right)\right)}\tag{12}$$

where *ScT*<sup>∞</sup> is the value of *ScT* at a distance far from the wall, and it is fixed with a value *ScT*<sup>∞</sup> = 0.85. It is considered to be a perfect mix at the tracer input boundary due to the distance with which it is injected.

The RTD curve is estimated by Equation (13) and it describes the distribution of the tracer through model at certain fixed time periods, from the injection of the tracer to the point of contact with the electrodes at outlets.

$$E(t) = \frac{\mathcal{C}(t)}{\int\_0^\infty \mathcal{C}(t)dt} \tag{13}$$

where *C (t)* is the concentration that is time-dependent on time. The development of RTD curves was executed taking *C (t)* from the Equation (10).

Both DTR curves obtained, experimental and simulated will have different dimensions, so it is necessary to make a parameterization. A normalized function *E* (*θ*) can be defined as [23]:

$$E(\theta) = \tau E(t) \tag{14}$$

where the dimensionless time *θ* is obtained by the spatial residence time *τ* = *V*/*Q*, *V* is the fluid volume of the model, and *Q* is the volumetric flowrate in each outlet.

The boundary and the initial conditions for solving the Equation (10) can be established as the tracer concentration, that is zero, *C* = 0, before the tracer injection in the reactor (*t* = 0) is performed. To simulate the tracer pulse injection, a mean form of the RTD is employed RTD at the outlets.

The mixing analysis was implemented in CFD, by the COMSOL Multiphysic software, 4.3b version in a computer of 64-bits operating system, Dual Core Xeon 2.30 GHz processor, and 96 GB of RAM. COMSOL has efficient simulation capabilities and complete graphic settings over other simulation commercial software packages [44]. This allows for better describing of different development patterns by function mathematical assignments in the simulation of tracer injection and the calculation of turbulent flow implemented in this study.

The discretization of the cross-junction volume was done using an unstructured mesh with tetrahedral, pyramid, and prism elements for simulation. According to the COMSOL manual [45], the mesh depends on the fluid flow and on the accuracy required. In many cases, the standard k-epsilon model with wall functions, that in this case, may deliver an accurate enough result at a much lower computational cost.

#### **3. Results**

#### *3.1. Mesh Selection and Sensibility Analysis*

The COMSOL software already has predefined mesh grades for simulations in turbulent flow, varying only the mesh density degree. A sensitivity analysis (Figure 6) was carried out among nine different fineness degrees to determine which density is adequate for the numeric representation regarding the measurements that were made of pressure and velocity in the real experiment. Also, it is verified that by means of the value of *y*<sup>+</sup> there is an adequate representation of the velocity profile.

**Figure 6.** Mesh sensibility analysis.

The sensitivity analysis for scenario V1 is shown in Figure 6. For the degrees of coarse mesh (extremely coarse, extra coarser, coarser, coarse, and normal), the speed is having a variation for both outlets. However, from the FINE fineness level, besides the speeds adjusted and the increase of the degree of fineness, there is no other significant change. Therefore, the appropriate degree of mesh for the simulations is the FINE grade, it is conformed of 493,438 to 444,610 mesh elements, with the minimum and maximum edge lengths of 1.92, 0.362 cm, respectively; the number of elements depends on the drop of pressure that is generated from the experimental valves that control the flow that in the CFD was represented with an obstruction on the flow. Although only one case of sensitivity analysis is shown, it is not necessary to perform one for each scenario, since the geometry of the model varies insignificantly.

According to Veerstek & Malalasekera [38] and Moukalled, Margani & Darwish [46], *y*<sup>+</sup> values in a range of (11.06 to 12.00) have important effects to turbulence production near the wall. As shown in the graph of Figure 6, the *u* velocity with respect to the mesh degree (by number of elements) stops varying from the FINE mesh degree. Therefore, this mesh degree was selected for the CFD model, which guarantees that the mesh tetrahedra established in the contours maintains a minimum *y*<sup>+</sup> value of 11.1 (those were verified plotting values of *y*<sup>+</sup> in a mesh graphic). The wall roughness was assumed to have a negligible effect. The solver that was employed was iterative, a generalized minimal residual, and a relative tolerance of accuracy of the CFD simulations considered a convergence criterion below 1 × <sup>10</sup><sup>−</sup>5. The typical solution around these meshes elements was unchanged.

#### *3.2. CFD Simulation*

The tracer injection must be defined in CFD simulations. However, the internal geometry of the sphere valve for the injection is controlled manually, and because of that, it is very irregular and difficult to obtain a chemical input pattern. Therefore, an interpolation was made between the tracer curves obtained at the outputs of the cross-junction. With this, a tracer input average is obtained at the injection point, and with this, the simulation is performed.

In the Table 3, it can be verified that the approximation between the experimental and CFD results has a significant proximity. For velocity values, the average error is 0.76% and the maximum value reached is 1.72%. The errors in terms of pressure are on average 0.80% on the N boundary and 0.92% on the W boundary. In part, pressure differences that are presented will have to do with the decimals with which the readings were recorded. The values in CFD have 16 decimals of floating point, while the oscilloscope used registers of an average value of signal, in a given time, with two decimals.


**Table 3.** Approximation errors between experimental and simulated tests.

The velocity values that were obtained in CFD correspond to vectors of a velocity field solved in the Navier Stokes Equation (3) for turbulent flow (*k-ε* model). However, to obtain an average value, COMSOL can estimate averages in a specific area (Figure 7). The selected mesh (FINE grade) approached the area of a contour (the outlet) of the pipe with good approximation (0.641% with respect to the real geometry obtained by a circumference).

**Figure 7.** Example of a velocity contour plot.

#### *3.3. Residence Time Distribution Curves*

Figures 8–11 show the RTD curves, experimental and numerical, at different inlet velocities (V1 to V4). The current intensity value (*milliAmperes*), recorded in time intervals, varies by the exchange of copper ions in the electrode cells for the experimental RTD. On the other hand, the numerical RTD is the tracer concentration through time. The parameterization was made according to the Equation (14). It is important to emphasize that the CFD curves maintain a constant value most of the time, which is the stabilization of the received signal giving the supplied potential. Once the curves of the experimental and simulated tests are derived, a comparison is made between the curves to analyze the representation of the simulation. The dimensions of the axes of the graphs are: The dimensionless time (θ) on the horizontal axe and the normalized function for the concentration E (θ) on the vertical axe; in all the cases, there are diverse values indicating a not complete mixing. Close agreement between simulation and experimental tests were attained. The shape of the RTD curves evidences that the mixture in the pipe junctions does not obey a perfectly mixed model, owing to the fact that the maximum of the peak is far from θ = 0. The maximum of the curves was found at values of θ, higher than of the unity, what indicates that the fluid elements loss kinetic energy in the cross-junction due to the coalition of the hydrodynamic streams and the bifurcation of the fluid flow.

**Figure 8.** Residence Time Distribution (RTD) curves for V1.

**Figure 9.** RTD curves for V2.

**Figure 10.** RTD curves for V3.

**Figure 11.** RTD curves for V4.

The variations between the experimental and theoretical RTD curves may be caused by the variability of the manual injection of the tracer, and tracer detection may be compared with the numerical convection-advection model on CFD. In almost all the cases, the peak of the experimental model is higher than in the CFD, resulting on the maximum concentration obtained in the experimental scenarios. Specially the Figure 11, where the RTD curve on the E outlet is too short when compared to the RTD curve on the S outlet, this is due to the outlet velocities that are registered on the experimental test, 1.25 m/s on the E outlet and 0.43 m/s on the S outlet, the tracer goes faster for the first velocity. The curves that were obtained for the validation tests V1 to V4 show that the tracer is following a tracer flow pattern that is similar to the experimental injection. To quantify these variations, an average error was calculated by the Root Mean Square Error (RMSE), Table 4.

Despite these differences, even the trend maintains a good proximity in terms of the distribution of the tracer, which shows that the proposed mathematical model is sufficiently calibrated model for describing the experimental mixture phenomenon in cross-junctions.


**Table 4.** Root Mean Square Error (RMSE) values for the four scenarios to evaluate variances in the RTD curves.

#### *3.4. Variation of ScT Coefficient*

A variation of the turbulent diffusivity is presented, in the search to verify if the model, under these conditions of pressure and velocity, varies significantly if the *ScT* coefficient is modified. The turbulent diffusivity was simulated in different values, using assignments to the Schmidt Turbulent *ScT* number. The range of values for this parameter is 0.61, 0.71, 0.81, and the value that was obtained by the Kays-Crawford model [41], which amounts to 0.5666.

The results indicate that the model is not sensitive to the variation of these *ScT* coefficients, since the parameterization of the curves reflects practically homogeneous results. This conclusion is verified, if a RMSE for each case is obtained again (Table 5).


**Table 5.** Approximation of RTD curves varying *ScT* values.

Table 5 shows that the variation is verified of the value used of *ScT* was based on the model of Kays-Crawford [41], and it was the most appropriate in most cases. The values of *ScT* have a very low impact in scale for the simulated model. In Figure 12, it is observed that the diffusion increases in the zones of greater recirculation, and there is less speed in the model. This should have a more significant influence, according to the bibliography. For the moment, it can be determined that the speed and pressure range far exceed the studied cases, and this tells us that at high levels of speed and pressure, the turbulent diffusivity does not have a significant impact on the transport of solute.

**Figure 12.** *ScT* and eddy diffusivity plots based on the Kays-Crawford [41] model (Equation (12)). High values appear in the turbulent flow zones.

#### *3.5. Incomplete Mixing Simulations*

Once the CFD model has been calibrated, the mixing grade can be evaluated for diverse scenarios with the tracer injection at the two inlets, N and W. To represent the development of the mixing graphically, it was introduced by two parameters: Inlet relation (CN/CW); in the same way, Outlet relation (CE/CS). The diverse scenarios implemented were made by the combination of concentration at the inlets, in intervals of 0.25 mol/L, in a range Inlet relation = [0.00–2.00]. In the case of Inlet relation = 0.00, it means that CN = 0.00 mol/L and CW = 5 mol/L. At the same way, for an Inlet relation = 1.50, it means that CN = 7.5mol/L and CW = 5mol/L. The scenarios are represented in Figure 13, with the results for the outlet relation after the cross-junction mixing.

**Figure 13.** Solute mixing for variable concentration at the inlets.

The four curves V1–V4 that are shown in Figure 13 represent the concentration relation between the two outlets E and S. Each curve (apparently straight lines) has different slopes. It can be observed that curves V1 and V2 have a similar development pattern, and are different from the curves V3 and V4. This is because V1 and V2 have closer input and output flow velocities, around 1 m/s for each one. However, curves V3 and V4 correspond to the input and output velocities with a wide degree of difference, between 0.427 to 1.302 m/s. This indicates that the transport of solute is by convective flow. This difference between velocities causes the tracer to flow at a higher rate at one output than at the other one. Therefore, the greater the difference between the input and output velocities, the greater the change in the slope of the generated curve.

Due to the behavior of the velocities mentioned above, the mixing behavior is maintained in a different proportion to the outlets in practically all of the cases. The complete mixing is not carried out on any scenario. Even the scenario with the same concentration in both inlets (Inlet Relation = CN/CW = 1) might not be a complete mixing. That is, since the result contains the same concentration at the inlets, it does not mean that a total mixture of the two incoming flows has been made. The complete mixing model assumes that "Outlet relation" must be equal to one in all cases, by the reason that the estimation assigns the same concentration at the two outlets. No one scenario has the behavior like a complete mixing.

#### **4. Discussion and Conclusions**

The physical characteristics of the experimental network, the pumping equipment, and the measuring devices used, allowed for achieving higher conditions of turbulence in the experimental model. To analyze the mixing that occurs in the cross-junction under these turbulence conditions, four scenarios were generated for flow velocities between 0.43 and 1.53 m/s; and for pressures around the cross-junction between 1.56 and 1.96 bar, we considered four kinds of similar velocities for the

inlets and outlets boundaries around the cross-junction. It is common for these conditions to occur at these variations in real drinking water distribution networks.

The RTD curves present a novel and an adequate approximation, having RMSE values that are between 0.0072 and 0.0509. It means that the average distance of separation in some points of the curves carries this value. These separations occurred for the most part in the peaks of the curves. However, in the ascending and descending curvatures, they present a very precise similarity. In most of these regions, slopes are well represented, in all validation tests.

The variation of *ScT* solved many questions about the role of turbulent diffusion and the impact that it generates in these operating conditions. The RTD curves showed a minimal change between them. The change is almost inappreciable graphically. Therefore, more checks were obtained using the RMSE error for each curve, and the variation is seen up to the fourth decimal place of precision of the RMSE. This allows for concluding that convection is the main transport in diluted species, and that the diffusion does not affect much in the simulation of the tracer course through the cross-junction. Even with this, it was found that the model that was proposed by Kays-Crawford [41] was the adequate parameter when presenting the minimum RMSE values with respect to the turbulent Schmidt variation.

The mesh selection and the geometric construction were adequate, because the values that were derived from flow and average velocity in a contour surface correspond approximately to those that were calculated with analytical formulas. Also, the simulation times did not require extensive periods, and convergence was achieved without the presence of broken curves in the graphics. The graphics could be processed with a reasonable degree of affinity in the surface of symmetry. The hydrodynamic approach retains a low degree of error with respect to experimental measurements. The errors that were obtained at the inlet boundaries, for the pressure was between 0.502% to 1.478% on the four scenarios and at the outlets boundaries, for the velocity, was between 0.133% and 1.720%. It not only depended on the simulation in CFD and the mesh, but also on the continuity adjustments that were made in the experimental tests. These adjusted tests ensure that the entry costs were equivalent to the costs of exit in the trials, which was also kept in the CFD environment since one of the principles of its coding is based on the continuity equation for three dimensions.

Finally, the mixing grade that was evaluated for the scenarios with tracer injection at the two inlets, varying nine combinations, shows that the mixing behavior is maintained in a different proportion to the outlets in practically all the cases. This is due to the behavior of the velocities that are mentioned above; the complete mixing is not carried out on any scenario and it was verified with the RTD analysis, where they show us on the graphics of the RTD curves how the tracer was obtained on the outlet boundaries at different times and peaks. It also could be verified on the scenarios with the same concentration at the inlets, obtain the same concentration, but it does not mean that a total mixture of the two incoming flows has been made, and that was validated with the same RTD curves. No one of the scenario has any behavior like complete mixing.

**Author Contributions:** For research articles with several authors, a short paragraph specifying their individual contributions must be provided. The following statements should be used "Conceptualization, D.H.-C., J.L.N. and J.M.R.; Methodology, D.H.-C., J.L.N., P.A.L.-J. and J.M.R. ; Software, J.L.N. and J.M.R.; Validation, D.H.-C., M.R., J.L.N., P.A.L.-J. and J.M.R.; Formal Analysis, D.H.-C., X.D.-G., J.L.N., P.A.L.-J. and J.M.R.; Investigation, D.H.-C. and J.M.R.; Resources, X.D.-G., J.L.N., P.A.L.-J. and J.M.R.; Data Curation, D.H.-C. and M.R.; Writing-Original Draft Preparation, , D.H.-C. and J.M.R.; Writing-Review & Editing, D.H.-C., M.R., X.D.-G., J.L.N., P.A.L.-J. and J.M.R.; Visualization, D.H.-C., X.D.-G. and J.M.R.; Supervision, J.L.N., P.A.L.-J. and J.M.R.; Project Administration, J.M.R.; Funding Acquisition, X.D.-G., J.L.N., P.A.L.-J. and J.M.R.", please turn to the CRediT taxonomy for the term explanation. Authorship must be limited to those who have contributed substantially to the work reported.

**Acknowledgments:** To CONACYT for the Master and Ph.D. scholarships (417824 and 703220) to D.H.-C. and the Ph.D. scholarship (294038) to M.R.; To Universidad de Guanajuato for the financial support of the project No. 100/2018 of J.L.N.; To Engineering Division, Campus Guanajuato and Geomatics and Hydraulics Engineering Department for the financial support of this project; and finally, to SEP-PRODEP and UG for the financial support to publish this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Computational Determination of Air Valves Capacity Using CFD Techniques**

#### **Salvador García-Todolí 1, Pedro L. Iglesias-Rey 1, Daniel Mora-Meliá 2,\*, F. Javier Martínez-Solano <sup>1</sup> and Vicente S. Fuertes-Miquel <sup>1</sup>**


Received: 7 September 2018; Accepted: 8 October 2018; Published: 12 October 2018

**Abstract:** The analysis of transient flow is necessary to design adequate protection systems that support the oscillations of pressure produced in the operation of motor elements and regulation. Air valves are generally used in pressurized water pipes to manage the air inside them. Under certain circumstances, they can be used as an indirect control mechanism of the hydraulic transient. Unfortunately, one of the major limitations is the reliability of information provided by manufacturers and vendors, which is why experimental trials are usually used to characterize such devices. The realization of these tests is not simple since they require an enormous volume of previously stored air to be used in such experiments. Additionally, the costs are expensive. Consequently, it is necessary to develop models that represent the behaviour of these devices. Although computational fluid dynamics (CFD) techniques cannot completely replace measurements, the amount of experimentation and the overall cost can be reduced significantly. This work approaches the characterization of air valves using CFD techniques, including some experimental tests to calibrate and validate the results. A mesh convergence analysis was made. The results show how the CFD models are an efficient alternative to represent the behavior of air valves during the entry and exit of air to the system, implying a better knowledge of the system to improve it.

**Keywords:** air valve; CFD; hydraulic characterization; entrapped air

#### **1. Introduction**

One of the main problems related to the operation and start-up of water distribution systems is the presence of air inside the pipes. There are many causes giving rise to the presence of air pockets: filling and emptying operations, temporary interruptions of water supply, vortexes in pumping feed tanks, air inlet in points with negative pressure, inflow in air valves during the negative pressure wave of a hydraulic transient, release of the dissolved air in the water, etc.

The presence and movement of air in water distribution pipes causes problems in most cases. Air pockets inside pipes can generate the following issues: reduction of pipe cross section, even blocks; the generation of an additional head loss, which increases energy consumption of pumping groups; decrease of pumps' performance; loss of filters' efficiency; noise and vibration problems; corrosion inside pipelines by the oxygen transported by the air; important errors in flow meters; etc. These problems originate an irregular system operation. However, the most important effect derived from the presence of air inside the pipes is the surge pressure caused by entrapped air pockets. These high pressures have their origin in the low inertia of the air and its high compressibility. The movement of

water, with much more inertia, can compress the air pockets and generate unacceptable pressures in the system.

Air movement inside the pipes results from a combination of different effects like the drag of water in its movement, the natural movement of the air pockets within the water caused by the difference of density, and the resistance of the air both by frictional effects and the effects of surface tension. This movement is a complex problem. Martin [1] developed the first study about important pressure surges than the air pockets can generate inside the pipes. He pointed out that the slug flow regime is the only biphasic flow that can lead to significant pressure surges. Later, Liou [2] developed one of the first models of filling pipes with irregular profiles. Previously, most models had been developed for very simple geometries: horizontal, vertical or inclined with constant slope. A more general model than the one proposed by Liou, including the presence of air valves, was developed by Fuertes [3].

One of the current research lines on the presence of air pockets in hydraulic transients is the modelling of filling and emptying processes of pipelines. During the filling process the volume of air pockets reached extreme values. In these conditions, significant increases in pressure can be generated. Hence the need to study models that represent the filling of pipelines [4–6].

More recently, the problem of the presence of air in the pipeline has been extended to drainage and sewer systems. The increasing intensity of rainfall makes rain drainage systems exposed to high water floods. These large flow rates result in extremely rapid filling of large pipes, many of which are not designed to support important internal pressures [7–9].

The most characteristic device for controlling the air in pipes is the air valve [10]. This device can perform three basic functions. The first one is to release small amounts of air that accumulate at high points of the system during normal operation. The control of accumulated air in the highest points of the pipes prevents reducing the cross section and the potential reduction of transportation capacity of the system. A second function is related to system ventilation during the filling and emptying processes of the pipelines. These processes require that important air flows are admitted or expelled from the system and so the air valves are the most suitable devices to perform this function. Finally, the third function is related to the protection of facilities against transient phenomena.

The hydraulic transient is crucial in the manufacture, design and installation of pipelines, but to date, these effects have often been overlooked or poorly studied, mainly due to the difficulty of their evaluation (only through the use of complex mathematical models). The devices generally used in pressurized water pipes to manage the air inside the pipe are the air valves, because under certain circumstances, these can be used as an indirect control mechanism of a hydraulic transient. However, specific technical manuals do not abound and there is no firm policy in this regard to help engineers make decisions to select the most appropriate air valve. Consequently, it is necessary to develop tools to understand the behaviour of an air valve in an installation.

Most previous work on air valves in water facilities has focused on their ability to reduce the effects of hydraulic transients. In this regard, Stephenson [11] relates the correct selection of air valve size and standpipe used with water hammer minimization; Bianchi et al. [12] propose several simplified (experimentally validated) formulas for dimensioning the required air valve size during filling of a system; De Martino et al. [13] study the transient caused by the expulsion of air through an orifice and deduce a simple relationship to predict the maximum pressure during the transient which agrees well with experimental data; and Fuertes et al. [14] present a methodology for dimensioning air valves to control transient phenomena, but considering the strong transient effects caused by the closing of the air valve's float [13]. Finally, some researchers have studied the use of aeration devices in other applications of hydraulic engineering. In this regard, Bhosekar et al. [15] study the use of aerators in spillways in an application similar to [11], but in another field.

In short, the analysis of the literature does not allow see any previous study related with the use of the computational fluid dynamics (CFD) techniques for the characterization of the behavior of the air valves. Unfortunately, one of the major limitations is the reliability of information provided by manufacturers. This information can only be validated at the present through trials [16]. Since air valve tests are extremely expensive, the search for alternative techniques for characterizing these elements becomes a key point in this scientific field. Therefore, the goal of this work is to determine the capacity of CFD techniques to predict the behavior of commercial air valves. The work was divided into two parts: the numerical study made using a CFD code developed by ANSYS Fluent and the experimental study of the behavior of the air valves. In both cases, it is intended to obtain an experimentally determined relationship between the mass flow *G* and the pressure *pt* inside the pipe where the air valve is connected. The final comparison of both studies (numerical and experimental) finally allows to see the validity of the proposed methodology as an alternative to the testing of this type of elements.

#### **2. Characterization of Air Valves**

The characteristic curve of an air valve is the relationship between the inhaled or exhaled mass flow *G* and the pressure *pt* inside the pipe in the point where it is connected. There are different models to characterize the behaviour of an air valve, but they all assume that the air inside the air valve has compressible fluid behaviour. Regardless of the chosen model, it is necessary to determine certain parameters of the air valve, which vary according to the mathematical expression used for the characterization.

Traditionally, one of the most commonly used models resembles the behaviour of the fluid inside the air valve to the isentropic flow that occurs in a convergent nozzle (Figure 1). Wylie and Streeter [17] made the analogy between nozzles and air valves, where the air behaviour into the pipe is considered isothermal, while the airflow throughout the valve (both inlet and outlet flow) is assumed to be adiabatic.

**Figure 1.** Isentropic flow in a convergent–divergent nozzle.

In this model, the air outlet mass flow *G* of an air valve, when the flow is subsonic is:

$$G = \mathbb{C}\_{exp} A\_{exp} p\_t^\* \sqrt{\frac{7}{RT\_t} \left[ \left( \frac{p\_{atm}^\*}{p\_t^\*} \right)^{1.4286} - \left( \frac{p\_{atm}^\*}{p\_t^\*} \right)^{1.714} \right]} \tag{1}$$

In the Equation (1) *Aexp* is the air valve output cross section; *p*∗ *<sup>t</sup>* is the absolute pressure in the pipe upstream the air valve, representing the input pressure *p*<sup>0</sup> in Figure 1; *R* is the air characteristic constant in the classical thermodynamic formulation of an ideal gas; *Tt* is the air temperature inside the pipe; *p*∗ *atm* is the absolute value of the atmospheric pressure; and *Cexp* is the outlet characteristic coefficient of the air valve. This coefficient takes values less than one, and the lower the value the higher the airflow resistance.

When the air velocity is greater than the speed of sound (Mach number greater than 1), volumetric flow rate is blocked because speed cannot increase more. Thus, the mass flow rate can be higher since an increase in pressure causes an increase in the density of air. In these supersonic conditions, the mass flow *G* is:

$$\mathbf{G} = \mathbf{C}\_{\exp} A\_{\exp} \frac{0.686}{\sqrt{RT\_t}} p\_t^\* \tag{2}$$

A similar approach can be carried out in the case of air inlet. In this case, there are two differences. On one side, the pressure in the inlet section is constant and equal to the atmospheric pressure. On the other side, the pressure in the output section is variable. This output section corresponds with the point of the pipe connected with the air valve. Therefore, the air inlet mass flow *G*, when the flow is subsonic is:

$$\mathbf{G} = \mathbb{C}\_{\text{adm}} A\_{\text{adm}} \sqrt{\mathcal{T} p\_{atm}^\* p\_{atm} \left[ \left( \frac{p\_t^\*}{p\_{atm}^\*} \right)^{1.4286} - \left( \frac{p\_t^\*}{p\_{atm}^\*} \right)^{1.714} \right]} \tag{3}$$

In Equation (3), *Cadm* is the inlet characteristic coefficient of the air valve, having the same considerations as *Cexp* in the output process; *Aadm* is air valve input cross section; and *ρatm* is the air density at atmospheric conditions.

If the flow is supersonic, the volumetric flow is kept constant. In this case, as the inlet pressure is constant (atmospheric pressure) the mass flow rate also remains constant. Thus, the flow is blocked. This means that even if the pressure inside the pipe decreases more, the amount of air admitted will not increase. Thus, in the conditions of supersonic flow *G* is:

$$G = C\_{adm} A\_{adm} \frac{0.686}{\sqrt{RT\_{adm}}} p\_{atm}^\* = constant \tag{4}$$

Expressions (1)–(4) are theoretical formulations of the potential behaviour of an air valve. The real behaviour of the air valve must be obtained through tests, and this information should be provided by the manufacturer. As an example, Figure 2 shows a graph built from the technical information provided by a manufacturer.

**Figure 2.** Air valves technical information provided a manufacturer.

Often, simplified expressions are used to treat numerically the characteristic curves of an air valve instead of the previous theoretical equations. The need to simplify these equations is related to decreasing the number of parameters on which the behaviour of the device depends. In this regard, Boldy [18] proposes a representation based on the equations of incompressible flow. Later, Fuertes [3]

performed a comprehensive review of the different models representing the behaviour of an air valve. All simplified models seek to reduce the number of parameters needed to characterize this, setting a relationship between the mass flow admitted or expelled and the pressure inside the pipe.

Normally, the manufacturers of air valves present a graph of the characteristic curve. This is the ratio between the admitted/ejected airflow and the difference in pressure between the inside and outside of the air valve. This curve is obtained with experimental tests in all possible operating regions. It is possible to obtain some equations from it, which relate the airflow to the difference in pressures. Unfortunately, the conditions under which these tests are performed are often not referenced in the technical information provided by the manufacturers, so it is nearly impossible to reproduce these tests in a laboratory. In fact, previous studies [3] have already revealed the discrepancies between the commercial data provided by the manufacturer and the actual data obtained by testing. Figure 3 represents the differences between the experimental data obtained in the laboratory and the curve obtained from the technical information provided by the manufacturer.

**Figure 3.** Comparison between laboratory tests and information provided by a manufacturer.

Results such as those in Figure 3 shows the need to find methodologies that allow characterizing the behaviour of air valves with sufficient reliability. This work proposes the use of CFD programs as an alternative, because these techniques are now considered as standard tools to predict the fluid flow behaviour. Calculations obtained using CFD will be compared with experimental tests performed in the laboratory.

The main problem of testing the characteristics of an air valve is the volume of air required. A 80 mm air valve can require a flow rate of 6200 m3/h measured under standard conditions. If the size is 100 mm the required flow is 9700 m3/h, while for 300 mm it may be necessary up to 87,000 m3/h.

Currently there are two main techniques to test an air valve. The first is based on the storage of large quantities of air in high pressure tanks. During the test the air is released gradually through a system that reduces the pressure to the usual operating pressures of these elements. The difficulty of this system is the tank volumes required are very high (32 m3 for 100 mm and almost 300 m3 for a 300 mm). An alternative option is to use a blower capable of providing the flow and pressure necessary to perform the test. The problem with these installations is the cost associated with the construction of this type of equipment. Simply as a reference data, to test a 300 mm air valve capable of dealing 24 m/s with a pressure of 0.9 bar, requires an approximate power of 1.4 MW. This data is an indicator of the complexity of this type of facility and justifies the fact that there are few places to carry out this type of test.

The experimental study was conducted at the Air Valves Test Bench of Bermad (Israel). This installation has a blower of 315 kW with capacity of 16,000 m3/h of air measured in standard conditions a 0.6 bar of differential pressure. At each test point the flow stability conditions were verified and subsequently the mass flow rate and the pressure at the inlet were recorded. The measuring devices have been previously calibrated guaranteeing errors below 0.5%. The system complies with all European standards and allows the testing of air valves between 2 and 12 inches in a continuous process that allows the measurement of air flow over the entire range of pressures considered. More details about the test bench and the collected data can be found in [16].

In this case, the tests were carried out on more than 20 different models from 13 different manufacturers, with flow ranges between −2500 and 3500 m3/h and pressure ranges between −0.57 and 0.58 bar. All the elements tested were air/vacuum valves 80 mm diameter, capable of exhausting air during pipeline filling and admitting large amounts of air if pressure drops below atmospheric.

#### **3. Computational Fluid Dynamics (CFD) Modelling**

CFD means the use of a computer-based tool for simulating the behaviour of systems involving fluid flow, heat transfer, and other related physical processes. CFD works by solving the Navier–Stokes equations over a region of interest, describing how the different properties (velocity, pressure, temperature, density, etc.) of a moving fluid change. CFD is used by engineers and scientists in a wide range of fields, including motor industry (combustion modelling and aerodynamics), buildings (thermal comfort, fire effects, or and ventilation), electronics (heat transfer within and around circuit boards) or medical applications (cardiovascular medicine). In recent years, CFD models have been applied for solving several hydraulic engineering problems. In this way, CFD has been used for a variety of purposes, as simulations of flow transient in pipes [19,20], behaviour of different control valves under different conditions [21–24], mixing models in water distribution systems [25–28], or for the characterization of elements in open channels [29–34].

CFD allows the analysis of flow patterns that are difficult, expensive or impossible to study using experimental techniques and although they cannot completely replace the measurements, the amount of experimentation and the overall cost can be significantly reduced. In addition, the level of detail that can be achieved with these techniques is very high, generating a lot of additional information without added costs. This allows more complex and complete studies of those that can be performed experimentally.

CFD techniques also have a number of disadvantages. On the one hand, the results of a CFD simulation are not 100% reliable, since it is necessary to simplify real physical systems and properly model complex phenomena such as turbulence. On the other hand, the CFD techniques require computers with great calculation capacity. Since computer capabilities have increased in the last years, CFD techniques provide a tool for determining the pneumatic characteristics of an air valve.

To perform a CFD analysis, there are three main stages: pre-processing, simulation and postprocessing. The pre-processing step includes geometry definition, mesh generation and boundary conditions definition. Once the physics problem has been identified, the simulation stage consists in solving the governing equations related to flow physics problems. Finally, the post-processing is the analysis of the results.

The geometry definition consists of specifying the shape of the physical boundaries of the fluid. Among other issues, it is necessary to define whether the computational model will be two-dimensional or three-dimensional. In this work, a standard computer-aided design (CAD) program is used to model a 3D commercial air valve. It should be noted that it is not necessary to include all the details of the air valve. Thus, given that the nominal diameter of the element is 80 mm, details of the geometry of size less than 1 mm have not been considered. Additionally, the geometry of the air valve (Figure 4) has been considered, so only half of it is represented. Likewise, in order to guarantee that the inflow into the air valve has been established, a straight pipe of the length equal to four diameters has been defined at the entrance.

**Figure 4.** Geometry and meshing structure of a commercial air valve for ANSYS Fluent analysis.

The second step is mesh generation, which consists in dividing the domain of the fluid into a number of smaller cells. The solver used (ANSYS Fluent) is based on the finite volume method, where the domain is divided into a finite set of control volumes. The general conservation equations for mass and momentum are solved on this set of control volumes.

For the representation of the behaviour of the air valves, a structured tetrahedral mesh has been selected, since this type of mesh allows a better adjustment than the hexahedral meshes, especially in the curved areas of the interior of the valve. In order to ensure proper convergence of the model, a preliminary analysis of the minimum mesh size required was performed. This analysis is performed considering the following boundary conditions: the mass flow rate of air is fixed in the inlet while the pressure remains constant in the output. Specifically in this study the required pressure at the inlet of the valve to generate a flow of 2000 m3/h of air under standard conditions, i.e., a mass flow of 0.667 kg/s, has been calculated.

All the simulations have been carried out by means of an implicit formulation in steady state. The resolution algorithm uses a density-based method, which has a coupled formulation of the equations of continuity, momentum and energy. This method of resolution allows better representation of the flow than pressure-based algorithms since the latter erroneously characterize negative pressure gradients. The numerical results of the convergence analysis of the mesh are collected in Table 1, where the value of the pressure P required at the inlet of the suction cup is shown based on the reference size used in the mesh and the number of elements it contains. As can be seen in Figure 5, mesh sizes around 4 mm are enough to analyze this type of problem with an adequate level of accuracy. Consequently, considering at least 300,000 cells guarantees an analysis with sufficient precision, as can be seen in Figure 6.


**Table 1.** Results of the mesh convergence analysis.

**Figure 5.** Analysis of the influence of the size of the mesh (size of the cells).

**Figure 6.** Analysis of the influence of the size of the mesh (number of cells).

Since a tetrahedral mesh has been used instead a hexahedral mesh, it has been necessary to review the aspect ratio of the different mesh elements. For the cell reference size finally selected (4 mm), 95% of the elements have a ratio of less than 2.5. Although only 5% of the cells have higher values, considering that less than 0.5% of the elements have an aspect ratio greater than 4.

The equations that describe the properties and motion of the fluid are numerically solved in each of the defined mesh cells. In general, these equations are obtained by applying the equations of continuity and momentum. The application of these laws allows obtaining the Navier–Stokes equations, which for a compressible flow, such as the air inside an air valve, can be expressed in a simplified way such as:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho \stackrel{\rightarrow}{\boldsymbol{\nu}}\right) = 0 \tag{5}$$

$$
\rho \frac{\partial \stackrel{\rightarrow}{\upsilon}}{\partial t} + \rho \left( \stackrel{\rightarrow}{\upsilon} \nabla \right) \stackrel{\rightarrow}{\upsilon} = -\nabla p + \rho \stackrel{\rightarrow}{\mathfrak{g}} + \nabla \cdot \mathfrak{r}\_{\overleftarrow{ij}} \tag{6}
$$

where *ρ* is the density of the fluid, <sup>→</sup> *v* is the velocity field, *p* is the pressure field, *g* are the forces per unit volume to which the fluid is subjected and *τij* is the stress tensor.

The Navier–Stokes equations are analytical and in order to solve them with a computer it is necessary to convert them into an algebraic system. This process is known as numerical discretization and there are different methods, being the most used finite difference, finite elements and finite volume.

Another of the basic aspects to consider when performing CFD analysis is the definition of the turbulence model. The objective of the turbulence models for the Reynolds-averaged Navier–Stokes (RANS) equations is to compute the Reynolds stresses and this is one of the basic aspects to consider when analyzing a flow through CFD.

There are several different formulations for solving turbulent flow problems such as Spalart–Allmaras, *k*-epsilon, *k*-omega, Menter´s Shear Stress Transport models, etc. All of these models increase the Navier-Stokes equations with an additional viscosity term. The differences between them are the use of wall functions, the number of additional variables solved and what these variables represent.

This work uses a *k*-*ε* two-equation turbulence model introduced by Launder and Spalding [35], where *k* is the turbulent kinetic energy and *ε* is the rate of dissipation of kinetic energy. The *k*-*ε* model has great applicability to turbulent air and water flows and is very popular for practical engineering applications, due to its good convergence and relatively less computational requirements, obtaining results with comparable accuracy to higher order models.

Actually, the *k*-epsilon model would be more properly called a family of models, because some variants have been developed for so many specific flow configurations. In this work, a feasible *k*-*ε* realizable model has been used, where the term realizable means that the model satisfies certain mathematical constraints on the Reynolds stresses, consistent with the physics of turbulent flows. The modelled transport equations to find the terms *k* and *ε* are:

$$\frac{\partial}{\partial t}(\rho k) + \frac{\partial}{\partial x\_{\dot{j}}}(\rho k u\_{\dot{j}}) = \frac{\partial}{\partial x\_{\dot{j}}} \left[ \left( \mu + \frac{\mu\_l}{\sigma\_k} \right) \frac{\partial k}{\partial x\_{\dot{j}}} \right] + G\_k + G\_b - \rho \varepsilon - Y\_M + S\_k \tag{7}$$

$$\frac{\partial}{\partial t}(\rho \varepsilon) + \frac{\partial}{\partial x\_j}(\rho \varepsilon u\_j) = \frac{\partial}{\partial x\_j} \left[ \left( \mu + \frac{\mu\_t}{\sigma\_\varepsilon} \right) \frac{\partial \varepsilon}{\partial x\_j} \right] + \rho \mathbb{C}\_1 \mathbb{S}\_\varepsilon - \rho \mathbb{C}\_2 \frac{\varepsilon^2}{k + \sqrt{v\varepsilon}} + \mathbb{C}\_{1\varepsilon} \frac{\varepsilon}{k} \mathbb{C}\_{3k} \mathbb{G}\_b + \mathbb{S}\_\varepsilon \tag{8}$$

In Equations (7) and (8) *uj* represents velocity component in corresponding direction (*xj*); *μ<sup>t</sup>* is the absolute dynamic viscosity of the air; *Gk* is the turbulent kinetic energy generated by the variations of the mean flow velocity components; *Gb* represents the kinetic energy generated by the buoyancy effects; *YM* is the contribution of the pulsatile expansion associated with the compressible turbulence; *C*1*ε*, *C*<sup>2</sup> and *C*3*<sup>ε</sup>* are constants; *σ<sup>k</sup>* and *σε* are the Prandl numbers for *k* and *ε* respectively; and *Sk* and *S<sup>ε</sup>* represent a global variation in time of the parameters *k* and *ε*, being defined independently of the rest of variables. Likewise, the constants *C1*, *η* and *S* are defined by:

$$C\_1 = \max\left[0.43, \frac{\eta}{\eta + 5}\right] \tag{9}$$

$$
\eta = S \frac{k}{\varepsilon} \tag{10}
$$

$$S = \sqrt{2S\_{ij}S\_{ij}}\tag{11}$$

On the other hand, turbulent viscosity *μ<sup>t</sup>* is calculated by a combination of *k* and *ε* by the equation:

$$
\mu\_t = \rho \mathbb{C}\_{\mu} \frac{k^2}{\varepsilon} \tag{12}
$$

where *Cμ* is a constant.

In this model, the constants *C*1*ε*, *C*2, *C*3*ε*, *σ<sup>k</sup>* and *σε* are adopted from the values proposed by Launder and Spalding [35], since this proved that they are very effective for turbulent flows with a wide range of Reynolds number for both water and air. The realizable model has shown some

improvements over the standard *k*-*ε* model, because it predicts more accurately the spreading rate of both planar and round jets. Likewise, it has also been shown that it performs better in flows involving rotation, boundary layers under strong adverse pressure gradients, separation and recirculation.

For a complete treatment of this problem, properties, initial and boundary conditions of the flow in space and time, would need to be specified. Regarding the properties of the fluid, it is necessary to specify viscosity, density and thermal properties. Regarding the flow inside the air valve, the viscosity of the air is considered constant, while the density varies admitting an air behaviour as if it were a perfect gas.

As initial conditions initial values for the variables are considered, from which the iterative process begins. The closeness of these initial values to the final solution has an influence on the process convergence time. In this work the initial pressure value is equal to the atmospheric pressure, while the initial air velocity is zero at all points.

The boundary conditions control the value of the variables or their relationship in the domain limits analyzed. It is basically about setting fixed values of pressure, velocity and temperature. On the model of an air valve of this work, the pressure is kept constant at the entrance and the exit. To analyze the air expulsion, this fixed inlet pressure will change in the different simulations, while the pressure in the outlet section will be constant and equal to the atmospheric pressure.

Once the model input values are completed, the software can solve the equations for each cell until an acceptable convergence is achieved.

The convergence of the method is based on the analysis of different criteria. Thus, as a general rule, the error is required to be less than 10−<sup>3</sup> in the continuity equation, in all the Navier–Stokes momentum equations for each direction (*x*, *y*, *z*), and in the model-specific equations of turbulence (*k* and *ε* in this case). Additionally, given the compressible nature of the fluid, the use of the energy equation is required. In this case, the degree of convergence is more demanding (10<sup>−</sup>6) than in the rest of convergence equations. Likewise, the fulfillment of a certain value of the residuals of the equations is not enough. Therefore, in order to ensure the convergence of the method, stability is required in the values of both the pressure at the inlet of the air valve and the mass flow that circulates inside the valve.

Finally, a post-processing software included in ANSYS is used to analyse the results generated by the CFD analysis. The data obtained can be analysed both numerically and graphically. In general, the most interesting graphic results that can be obtained is the distribution of pressures inside the air valve (Figure 7).

**Figure 7.** Pressure distribution in the midplane of the air valve.

Figure 7 allows analyzing in detail the different parts of the interior of the air valve where greater energy is lost. From the point of view of the design of the device, these zones are enlisted in areas where the pressure gradients are greater. Undoubtedly, improving the capacity of admission or expulsion of an air valve involves improving the design of these points. However, the objective of this work is not the design of this type of device, but to determine the capacity of the CFD techniques to predict their capacity of admission and expulsion.

#### **4. Results**

The application of the described analytical methodology allows comparing results obtained computationally with those obtained experimentally. For the particular case of one of the air valves studied, the results considering air as expulsion flow are shown in Figures 8 and 9.

**Figure 8.** Comparison of the curve of an air valve: experimental data vs. computational fluid dynamics (CFD) analysis.

**Figure 9.** Comparison of the curve of an air valve: experimental data vs. CFD analysis.

The Figure 8 considers the representation of pressure and flow values. This representation demonstrates the good correspondence between the results obtained through the CFD simulation and the results obtained through laboratory tests. Additionally, Figure 9 represents the pressure values of the CFD results against the laboratory tests. This allows comparing both methodologies in a more detailed way. The results demonstrate the validity of the methodology based on CFD to be able to predict the behavior of an air valve.

Even considering the goodness of the results in Figure 9, to validate the CFD-based methodology of this work as an adequate technique for representing air inlet and outlet characteristics in an air valve, two additional questions need to be answered. On the one hand, it is necessary to validate the methodology for other geometries. On the other hand, it is also required to validate the method when it is applied in the case of air intake from the outside to the inside of the pipeline.

Related to this, the CFD analysis was performed on six different models of air valves of the same nominal diameter (DN 80 mm) but with different internal configurations and different air outlet mechanisms (down, side, mushroom). The results of these analyses are presented in Figure 10.

**Figure 10.** Comparison between laboratory tests and CFD simulations of six different models of air valve. Air expulsion flow.

The results show a good fit between CFD simulations and laboratory tests. However, all tests are performed up to a certain maximum pressure value. This maximum pressure is conditioned in each case for different reasons. In some cases, some of the air valves did not allow further extension of the air expulsion tests, since the flow that causes the dynamic closing of the valve was reached. That is to say, the air valve closes without the presence of water simply by the efforts caused by the air current that overcome the weight of the float. In other cases, the analysis is not extended because the maximum capacity of the test equipment is reached, which is 3000 m3/h of free air or 0.5 bar of differential pressure.

Once the capacity of the CFD model to represent the behavior of the air valve in ejection is verified, the process is repeated in the case of air intake. It should be noted that the compressibility effects of air are more noticeable when the pressures are lower. In this sense, it is important to highlight

the capacity of the CFD model to represent the sonic block that is reached in the air valve when the pressures approach values close to −0.5 bar.

Figure 11 shows the results of the comparative study of six models of air valves in the intake phase. As can be observed, there are hardly any differences between the results obtained by laboratory tests and by CFD techniques.

As in the study of the expulsion of air, there is an analysis limit of 0.5 bar differential pressure between the outside of the pipeline and the inside of the pipeline. There are two reasons for this limitation of the study. The first is the maximum test capacity of the test bench used. The second is that for differential pressures greater than 0.5 bar below atmospheric pressure, the so-called sonic block is produced. In other words, the volumetric flow admitted by the suction cup reaches its maximum value. Consequently, as shown in Figure 11, both the results of the CFD model and the laboratory tests tend asymptotically to a value of maximum intake flow and this value will not increase although the depression inside the conduction increases.

**Figure 11.** Comparison between laboratory tests and CFD simulations of six different models of air valve. Air intake flow.

#### **5. Conclusions**

CFD techniques provide a tool for determining in detail the behavioral characteristics of hydraulic elements. However, CFD requires a proper calibration of the models, so it will be essential to validate the results by tests in the laboratory. Specifically, this work shows the effectiveness of these techniques to represent the behavior of air valves that are installed in water supply networks. According to the results, it is possible to state the following:


closure has only been validated through laboratory tests. This is the reason why in some cases the curves of the air valves obtained through CFD cover a range of flow greater than those obtained experimentally.


To summarize, it can be concluded that the CFD models are an efficient alternative to represent the air valves during the entry and exit of air to the system. Given the difficulty and the costs of testing this type of element, the application of these techniques is an adequate alternative for the characterization of these elements. In the same way, the methodology could be very useful to verify the technical information provided by the manufacturers.

**Author Contributions:** All authors contributed extensively to the work presented in this paper. S.G.-T., P.L.I.-R. and D.M.-M., contributed to the subject of research, the modelling, the data analysis, and the writing of the paper. F.J.M.-S. and V.S.F.-M. contributed to the experimental measures and the manuscript review.

**Funding:** This research was funded by the Program Fondecyt Regular, grant number 1180660.

**Acknowledgments:** This work was supported by the Program Fondecyt Regular (Project 1180660) of the Comision Nacional de Investigación Científica y Tecnológica (Conicyt), Chile.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

CFD Computational Fluid Dynamics RANS Reynolds-averaged Navier–Stokes

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Case Report*
