Incomplete Mixing Model at Cross-Junctions in Epanet by Polynomial Equations

Abstract

In Water Distribution Networks (WDN), the water quality could become vulnerable due to several operational and temporal factors. Epanet is a hydraulic and water quality simulation software, widely used, to preserve the control of chemical disinfectants in WDN among other capabilities. Several researchers have shown that the flow mixing at Cross-Junctions (CJs) is not complete as Epanet assumes for the cases of two contiguous inlets and outlets. This paper presents a methodology to obtain the outlet concentrations in CJs based on experimental scenarios and a validated Computational Fluid Dynamics (CFD) model. In this work, the results show that the Incomplete Mixing Model (IMM) based on polynomial equations, represents in a better way the experimental scenarios. Therefore, the distribution of the concentration could be in different proportions in some sectors of the network. Some comparisons were made with the complete mixing model and the Epanet-Bulk Advective Mixing (BAM), obtaining relative errors of 90% in some CJs.

1. Introduction

Chlorine disinfection is one of the most important potabilization techniques in drinking water treatment processes. This chemical agent is mainly used to destroy pathogenic organisms, but not only at the treatment, it also remains active throughout the network. However, it is also true that chemical-based disinfection has caused unwanted risks due to byproducts, such as trihalomethanes (THM) (among others; cf. [1,2,3,4]) generated by reactions with organic matter present in water. In order to control these and other situations, engineers employ diverse tools, including hydraulic and water quality simulation programs. Epanet is one of the most widely used software to simulate the behavior of disinfectant distribution through WDN [5]. However, it is very important to have an accurate model for predicting water quality so that it is as reliable as possible.

1.1. The Complete Mixing Model

Epanet, as in the most hydraulic simulation software, assumes that the mixing at pipe junctions is complete and instantaneous [2,5]. Epanet calculates the concentration of a substance that leaves the mixing junction by the weighted average of the concentrations of the inlet pipes (1), as described in Figure 1.

For a $k$ junction:

$$C_i=\frac{{\sum }_{j\in I_k}{Q}_j{C}_j+Q_{k,ext}{C}_{k,ext}}{{\sum }_{j\in I_k}{Q}_j+Q_{k,ext}}$$

(1)

• C_i represents the concentration at the start of link i (mg/l).
• Q_j is the flow rate at link j (l/s).
• Q_k,ext is the external source flow entering the network at node $k$ (l/s).

• C_k,ext is the concentration of the external flow entering at node $k$ (mg/l).

• I is the link with the flow leaving node $k$.

• I_k is the set of links with the flow into $k$ (in Figure 1, I_k are the pipes with labels Q_j C_j).

In cases where CJs have two contiguous inlets and outlets, several researchers have shown in physical models and CFD numerical simulations that the mix is not complete and instantaneous. Most of the investigations were based on the use of tracers for the distribution of concentrations.

Two of the first researches described physical and CFD scenarios (2D and 3D modeling for crosses of 50.8 mm diameter and 0.78 m/s of flow velocities), obtaining that the tracer can derivate only 12–14% to one inlet, contrary to Epanet predictions of 50% in each outlet [6,7]. This research concluded that a higher simulation precision is obtained using 3D models solved by the Reynolds Averaged Navier Stokes (RANS) equations compared to the 2D modeling. Apart from this, the solute transport in cross and double-tee junctions under Reynolds (Re) of 20,000 was solved by a Large Eddy Simulation with pure advection [8]. Under this approach, it is not necessary to consider the turbulent Schmidt number and it is better to place a special emphasis on the meshing of the model. A similar conclusion was mentioned in [9]: Convective transport in high turbulence conditions does not have a significant influence by the variation of the Turbulent Schmidt Number.

Contrasting flows (Re_inlet1: 1500/Re_inlet2: 40,000, equal flow at the outlets) were studied in [10]. In these cases, the solute mix can be almost null, from values of 4% to values closer to 0% of the concentration at one pipe outlet. Their best results were obtained by refining the mesh on the model walls and adjusting internal cross geometries. A CFD simulation analysis was made in [11], considering the influence of the turbulent diffusivity parameters for different velocity settings (3, 4, and 5 m/s) in a 0.040 × 0.040 m diameter cross. These results showed that the mix was imperfect in all the scenarios analyzed. Mixing effects can be greater when the turbulent Schmidt is lower. In [12,13], some mixing scenarios were carried out using chlorine in different cross configurations. Incomplete mixing was done where the inlet and outlet boundaries were opposite only. A similar conclusion was written in [14], reporting that contiguous inlet and outlet cases distributed chlorine about 20% in one outlet. However, all the scenarios were carried out with flow rates lower than 0.5 l/s.

1.2. Incomplete Mixing Models (IMMs)

There are several researches that propose mathematical functions to estimate the solute mixing phenomena. Complementary codes have also been developed for Epanet in order to carry out water quality simulations considering incomplete mixing.

In [15], the authors investigated the solute transport through cross and tee-junctions in a scaled 3 × 3 laboratory pipe network (0.0127 m in diameter and 0.91 m in length) using a conductivity tracer. The formulated mixing models were implemented in the Epanet-BAM that can solve the water quality simulation based on mixing coefficients that derivate the incomplete solute distribution through cross outlets. The simulations reached flow rates up to 3.85 l/s, with Re around 24,000. The pollutant intrusion from a source was applied to evaluate the impact for the possible consumers affected by the contaminated water. Due to the hydrodynamic and the pollutant concentrations, there were more affected consumers by the Epanet complete mixing model than the Epanet-BAM model. However, a contrasting work was described in [16], using two real high interconnection WDNs and a totally meshed hypothetical network (over 780 CJs). These researchers applied Epanet-BAM and evaluated them together with a “combined” mixing model (combination of complete and incomplete mixing mediated by a scale factor “s”), described in [17]. The results of [16] for both real networks did not exceed the 5% of the average of absolute difference in respect to the complete mixing. For the hypothetical network they could reach differences up to 50% when the water sources were relatively close. However, if they were positioned in opposite places of the network, the differences were uniformed to 10%, so they concluded that in most cases there is no significant difference between the complete and IMMs.

In [18], some CFD experiments were carried on CJs with different pipe diameters and variable flows. They developed incomplete mixing equations for three diameter ratios combinations, fitted with the Particle Swarm Optimization. They concluded that the dominant parameters are Reynolds ratios to the inlets and outlets and diameter ratios for the inlets. They mentioned that the tracer inlet concentration is mostly conserved towards its contiguous outlet. This conclusion is also cited in [14]. These authors performed experimentation in the same experimental assembly, but in this case, they worked directly analyzing different flows on the same cross (0.0254 m diameter) rather than the varying Reynolds number ratios with different pipe diameters. In the results exclusively for CJs, triphasic equations were formulated to obtain the concentration at the outlets. Once again, they concluded that the tracer from one inlet was mostly conserved to the adjacent outlet. In [19], they focused their studies on laminar and transitional flows. They determined three zones as a function of the Reynolds number, for Re greater than 3000, the mix is constant and invariant under flow variations, so that at higher Reynolds values, the degree of mixing decreased.

Researchers from the University of Arizona, in collaboration with the University of Korea, assembled a laboratory pipe model to demonstrate the differences between simulations with complete and incomplete mixing [20,21]. They built a 5 × 5 square network with nine CJs and 32 pipes of 0.016 m diameter. Three scenarios were configured achieving an incoming total flow in the network of 1.2 l/s, with average Re values of 9000, having laminar, transitional, and turbulent flows in the three scenarios. An error prediction study was carried out by monitoring the concentration at different measurement points. The model with incomplete mixing resulted with a 15% prediction error, while Epanet with the complete mixture exceeded 66%. The complete mixing model generated a more uniform distribution of the concentrate throughout the network, speculating that the solute remains equivalent in most points. Their results with the IMM generated a “diagonal” narrow plume, where in some areas the substance was more concentrated. From these experiments AZRED II emerged, a software derived from Epanet that allows simulating while considering the IMM [22].

Table 1. Experimental (EXP) and numerical (CFD) parameters used in different mixing studies.

REF *	Diameter (mm)		Velocity (m/s)		Flow (l/s)		Reynolds		EXP/CFD	Authors
	Min	Max	Min	Max	Min	Max	Min	Max
[6]	50.8		0.87		1.76		44,111		EXP	Van Bloemen et al., 2005
	50.8		0.78		1.58		44,388		CFD
[7]	50.8		0.79		1.595		40,000		CFD	Webb and Van Bloemen 2006
[23]	19.1		0.38	2.53	0.11	0.72	7151	48,122	EXP	Romero et al., 2006
	50.8		0.87		1.76		44,000		CFD
[10]	25.4		0.06	1.58	0.03	0.80	1500	40,000	CFD	Webb 2007
[17]	12.7		0.24	0.40	0.03	0.05	3100	5013	EXP/CFD	Ho et al., 2007
[24]	12.7	50.8	0.87	1.73	0.11	3.51	11,000	88,000	EXP	Romero et al., 2008
	305		0.03	0.07	2.05	5.26	8561	21,698	CFD
[15]	203.2		0.09	0.119	3.09	3.85	19,304	24,157	CFD	Ho and Khalsa 2008
[21]	16		0.14	1.63	0.03	0.33	2295	26,147	EXP	Song et al., 2009
[25]	25.4	50.8	0.16	0.45	0.08	0.23	4000	13,000	EXP	Ho and O’ Rear 2009
[26]	76.2	152.4	0.01	0.03	0.12	0.26	2100	-	CFD	Andrade et al., 2010
[11]	40		3.00	6.00	3.77	7.54	120,000	240,000	CFD	Liu et al., 2011
[18]	250		0.24	0.72	11.78	35.34	60,000	180,000	CFD	Yu et al., 2014
[14]	25.4		0.61	3.49	0.310	1.770	15,545	88,722	EXP/CFD	Shao et al., 2014
[13]	12.7	38	0.79	1.34	0.17	0.9	17,043	30,172	EXP/CFD	Mompremier et al., 2015
[12]	17.5	43.8	0.27	0.55	0.11	0.21	6050	12,100	EXP/CFD	Mompremier 2015 (PhD Thesis)
[19]	16	50	0.03	0.07	0.01	0.14	500	3500	EXP/CFD	Yu Shao et al., 2019
[8]	50.8		0.4		0.81		20,000		CFD	Luka et al., 2020

* REF values are the numbered references in this article.

shows a summary on the variables used in experimental and numerical scenarios of mixing at CJs from 2005. In this research, experimental and CFD scenarios in two stages were proposed. First, an experimental research with four scenarios with a CJs of 0.1016 m diameter was made, using an electrochemical tracer. The formulated CFD model was validated in [9]. In this paper, eight more scenarios were presented, modeling the hydraulics in a range of velocities of 0.43 to 2.48 m/s, reaching flows from 6 to 12 l/s, and Reynolds from 80,000 to 250,000. The tracer transport was simulated with the CFD validated model of the first four scenarios, to obtain the IMM by polynomial equations. Considering that these dimensions were not experimentally presented before in 0.10 m diameter pipes, the IMM presented in this work could be represented in real WDN under similar operating conditions.

Improving the precision of water-quality simulation could increase the capabilities of other applications. Researchers in [27] developed a chlorine optimization model based on genetic algorithms to determine the best location of booster chlorination stations, in order to optimize the use of disinfectants in WDN. In [28], the authors presented a study in which they determined the strategic areas for the location of sensors of contaminants, in order to make an early detection of a possible accidental or intentional distribution of harmful chemical agents to prevent the risk of drink contaminated water. Both examples predicted the best points for the location of water quality control devices. However, they depended on the way in which Epanet simulates the chemical distribution throughout the network. The predictions could not be accurate enough, as there are diverse CJs in the studied networks.

1.3. Goals and Improvements of the Current Work

This paper presents an analysis of mixing at CJs using a validated CFD model, presented in [9]. This analysis includes the variation of concentrations at the cross inlets in twelve scenarios at different operating conditions (flow velocity and pressure). With these CFD simulations, mixing parameters are obtained and they allow estimating the concentration at the outlets that occurs in the CJ. The lack of precision in water quality simulation can potentially lead into an improper design of monitoring systems, and the location of re-chlorination systems in unsuitable areas. Improving the accuracy in these models is also necessary to simulate the spatio-temporal dispersion of chemical and microbial agents during accidental or intentional contamination events. The advantage of this research is the implementation of IMM by polynomial equations on Epanet, using a code that allows obtaining the concentrations at the outlets of the CJs in each quality time step of the Epanet simulation. Most of the previous researches were experimental and numerical models that simulated the mixing on specific CJs. Only two previous works managed to simulate their results on Epanet (Epanet-BAM and AZRED II). In those cases, the concentrations at the outlets are based on fixed initial values during all the simulation. The IMM based on polynomial equations will be compared with the complete mixing model of Epanet and the Epanet-BAM. For the second one, Rossman uses it as an option to include non-ideal mixing at pipe junctions [29].

2. Materials and Methods

2.1. Experimental Cross and the CFD Model

The experimental model was carried out on a 0.1016 m galvanized iron cross from an experimental WDN of the University of Guanajuato. Electromagnetic flow meters and a digital storage oscilloscope with four measuring channels were used to monitor the flow and pressure values. The WDN can conduct flows up to 30 l/s at pressures of 1.78 bar. The CFD model of the cross was formulated considering the main elements of the experimental model. The experimental model, the geometry, and mesh specifications of the CFD model are shown in Figure 2. The CFD module of the COMSOL Multiphysics software, version 4.3b, was used to solve the RANS equations according to the standard k-ε model, which works with the high Reynolds number turbulence. This model was validated by the stimulus-response technique, generating Residence Time Distribution curves according to the experimental tracer measurements of four scenarios. The formulation and validation details, for the hydraulic and the outlet concentrations, are described in [9].

In this research, eight scenarios were proposed with different hydraulic conditions, all performed in the experimental WDN and in the validated CFD model in order to obtain twelve different scenarios for the IMM by polynomial equations. In these new scenarios, only the hydraulic validation was performed due to the results obtained in [9], showing that the convection is the main transport of diluted species for CJs, and the diffusion does not affect the tracer trajectory. This research analyzes the cases of mixing from two contiguous inlets and outlets. The four boundaries are named following the cardinal directions: North (N), West (W), East (E), and South (S), where N and W are the inlets, and E and S, are the outlets. The velocity and pressure conditions of the scenarios are described in

Table 2. Velocity (v, m/s) and pressure (P, bar) values on the inlet and outlet boundaries for the twelve scenarios proposed (S1–S12).

Boundary		S1		S2		S3		S4		S5		S6
		v	P	v	P	V	P	v	P	v	P	v	P
Inlet	N	1.12	1.60	0.93	0.90	1.13	1.92	1.11	1.88	1.18	1.60	1.41	1.72
	W	1.27	1.60	1.05	0.90	1.04	1.91	1.06	1.88	1.22	1.60	0.82	1.72
Outlet	E	1.06	1.56	1.03	0.90	1.20	1.91	1.00	1.88	1.07	1.60	1.11	1.72
	S	1.32	1.56	0.96	0.90	0.96	1.91	1.17	1.88	1.33	1.60	1.12	1.72
Boundary		S7		S8		S9		S10		S11		S12
		v	P	v	P	V	P	v	P	v	P	v	P
Inlet	N	0.81	1.57	1.34	0.60	1.07	1.60	2.23	0.66	0.67	* 1.96	1.71	0.32
	W	1.53	1.56	2.06	0.60	1.26	1.60	1.33	0.66	0.99	* 1.96	2.06	0.32
Outlet	E	1.04	1.56	1.14	0.57	0.99	1.60	1.08	0.63	1.25	* 1.96	1.97	*0.28
	S	1.30	1.56	2.26	0.53	1.33	1.56	* 2.48	0.62	* 0.43	* 1.96	1.80	0.29

* V varies from 0.43 to 2.48 m/s and 0.28 to 1.96 bar for P.

.

For each simulation of Si, the inlet concentrations will be varied numerically on the CFD model at the cross-boundaries N and W according to

Table 3. Concentration variation at the inlets C_N and C_W (mol/m³) for the proposed scenarios.

CN	0	62.5	125	187.5	250	312.5	375	437.5	500
CW	250	250	250	250	250	250	250	250	250
IN = CN/CW	0.00	0.25	0.50	0.75	1.00	1.25	1.50	1.75	2.00

, in order to obtain 108 diverse scenarios.

The IN coefficient is the ratio of the concentrations at the inlets, Equation (2):

$$IN=\frac{C_N}{C_W}$$

(2)

where C_N and C_W are concentrations at the N and W inlets, respectively. Similarly, the OUT coefficient is the outlet concentrations ratio. It will be obtained by the following equation:

$$OUT=\frac{C_E}{C_S}$$

(3)

where C_E and C_S are concentrations at the East and South outlets, respectively. These coefficients are the objective of the IMM by polynomial equations to find the concentrations at the outlets.

2.2. Concentration Outlets Using the OUT Coefficient

The OUT coefficient can be estimated analytically. This coefficient will be used to find the concentrations at the E and S outlets by the following conservation equation:

$$C_N{Q}_N+C_W{Q}_W=C_E{Q}_E+C_S{Q}_S$$

(4)

where C and Q denote the concentration and flow, respectively and according to their boundaries N, W, E, and S.

From Equation (3):

$$C_E=\left(OUT\right)C_S, C_S=\left(OU{T}^{-1}\right)C_E$$

(5)

From Equations (4) and (5):

$$C_E=\frac{C_N{Q}_N+C_W{Q}_W}{Q_E+(OU{T}^{-1})Q_S}$$

(6a)

$$C_S=\frac{C_N{Q}_N+C_W{Q}_W}{\left(OUT\right)Q_E+Q_S}$$

(6b)

Using Equations (6a) and (6b) and the OUT coefficient, the concentrations at the outlets can be computed using the resulted flows in the four boundaries as well as the concentration of the inlets from every CJ with contiguous two inlets and two outlets of the Epanet simulations.

2.3. The IMM by Polynomial Equations

To apply the IMM by polynomial equations, the operating flow conditions at each cross are compared with one of the twelve scenarios previously proposed (Table 2) and selecting the one that best approximates the data. The rescaled proportion of the inlet and outlet flows (Q_r) of the cross are obtained as follows:

$$Q_{rIN}=\frac{Q_N}{Q_W}$$

(7a)

$$Q_{rOUT}=\frac{Q_E}{Q_S}$$

(7b)

These values are also calculated for the twelve scenarios (S1 to S12) of this paper, named $Q_{rIN,Si}$ and $Q_{rOUT,Si}$ (the added suffix Si corresponds to the twelve scenarios of this study, i = 1, 2… 12). For example, converting the four boundary velocities of scenario S1 (Table 2) to l/s, the flow rates are: Q_N = 9.048, Q_W = 10.292, Q_E = 8.610, and Q_S = 10.730 (l/s). From Equations (7a) and (7b), $Q_{rIN,S1}$= 0.879 and $Q_{rOUT,S1}$= 0.802. The ratios $Q_{rIN,Si}$ and $Q_{rOUT,Si}$ for the twelve scenarios are shown in

Table 4. Ratios Q_rIN,Si and Q_rOUT,Si for the twelve proposed scenarios.

	S1	S2	S3	S4	S5	S6
QrIN,Si	0.879	0.882	1.049	1.085	0.962	1.730
QrOUT,Si	0.802	1.069	0.861	1.247	0.806	0.986
	S7	S8	S9	S10	S11	S12
QrIN,Si	0.527	0.652	0.851	1.681	0.680	0.832
QrOUT,Si	0.799	0.502	0.746	0.436	2.921	1.098

. In general, the variables $Q_{rIN}$ and $Q_{rOUT}$ describe the proportion of the flows at the inlets and outlets of the CJ.

The scenario to be studied will be associated with one of the twelve scenarios proposed using Equation (8):

$$R_i=\left|Q_{rIN}-Q_{rIN,\mathrm{S}i}\right|+\left|Q_{rOUT}-Q_{rOUT,\mathrm{S}i}\right|$$

(8)

where Q_rIN and Q_rOUT are the ones of the Epanet simulation and Q_rIN,Si and Q_rOUT,Si will be the results according to each scenario S1 to S12. The minimum value of the R_i gives the most accurate scenario Si (similar in its proportion of the inlet and outlet flows). Then, the IN coefficient is calculated from Equation (2) and is used in the corresponding polynomial equation resulting from the lowest value of Equation (8), in order to obtain the OUT coefficient. Therefore, C_E and C_S can be given by Equations (6a) and (6b).

Pseudo algorithm steps to approximate C_E and C_S concentrations:

• To solve the hydrodynamics of the network in the current hydraulic time step, in order to identify the CJ with two contiguous inlets and two contiguous outlets. For each quality-time step, the initial Q and C values will be: Q_N, Q_W, Q_E, Q_S and C_N, C_W, respectively. If the algorithm does not detect the two contiguous inlets and outlets at the CJ, then the concentrations are calculated using the Epanet complete mixing model.
• To calculate the Q_rIN and Q_rOUT ratios with Equations (7a) and (7b) for each CJ in step 1. They are compared to Q_rIN,Si and Q_rOUT,Si from the twelve scenarios and the chosen scenario Si will be the one with the lowest value of R_i (Equation (8)).
• The IN coefficient (Equation (2)) is calculated and evaluated in the selected polynomial equation from the previous step to obtain the OUT coefficient.
• Equations (6a) and (6b) are applied to determine the concentrations at the cross outlets C_E and C_S.
• To solve all the quality-time steps involved in the hydraulic time step and repeat from step 1 until the simulation is over.

Computational Algorithm to Apply the IMM by Polynomial Equations

A programming code was developed using the Epanet Matlab Toolkit [30], which locates and analyzes the distribution of flow rates at CJs. Auxiliary junctions and bypass pipes of short length are incorporated into the located CJs, maintaining the same parameters of their components (elevation, diameter, roughness, etc.). These auxiliary junctions and bypass pipes do not interfere with the hydraulic operation of the network. In each Epanet quality-time step, the C_E and C_S are calculated using the methodology described in the previous section. The resulting outlet concentrations are assigned, by setpoint boosters added by this programming code, at the respective outlet boundaries.

2.4. Application Examples

2.4.1. Tests in a Proposed Network

The main model network is built with 47 pipes connected with 30 junctions, nine of these are CJs with diverse pipe diameters (Figure 3). All pipes are 600 m length and 1.5 × 10⁻⁵ m of Darcy-Weisbach roughness. The example network operates as an hydraulic-stationary regime. This means that there are no changes in node demands, consequently, there is no variation in the flow rates or pressures, but it simulates several hours for modeling the water quality.

Three water quality cases, EX1, EX2, and EX3 are shown in Figure 4. In all cases, two Reservoirs (RVs) with a source concentration are located in different positions. The EX1 scenario was performed to analyze the distribution of a non-reacting chemical during the flow simulation, in order to represent a possible contaminant event in WDNs. In the EX1 scenario, the reservoir RV1 supplies 1.5 mg/l of a substance, while RV2 provides clean water (Figure 4). EX1 has two demand nodes located at the lower corners of the network, with demands of 10 l/s (left) and 20 l/s (right), to generate various operational mixing flows on the nine CJs with two contiguous inlets and outlets without demand.

The EX2 and EX3 scenarios were carried out to analyze the differences between the mixing models representing the disinfection on WDN considering water consumption in all the nodes. The concentrations of the chemical in the RVs represent the values of chlorine supply by the water treatment plants (according to Mexican Normative) and the reaction coefficients are applied to simulate the decay of residual chlorine. The World Health Organization recommends that the concentrations for chlorine on all nodes should be between 0.2 to 0.5 mg/l, although in Mexico government regulations specify a maximum value of 1.5 mg/l. Subfigures EX2 and EX3 have different positions for the RVs. All the junctions have a demand of 1 l/s. The RVs concentrations are 0.65 and 1.20 mg/l, and the reaction coefficients kb and kw used are −1.512 d⁻¹ and −0.275 m/d, respectively. The kb coefficient was obtained in a real network from the city of Guanajuato, Mexico [31] and kw is an average of the findings of some authors [32,33,34].

2.4.2. Epanet Example Network NET3.net

The IMM by polynomial equations is applied in the example network Net3.net proposed by Rossman [5]. This network is a more complex model trace and its principal characteristics are: (a) two RVs, three tanks, two pumps, 117 pipes, and 92 junctions (nine of them are CJs); (b) one general demand pattern, four patterns at certain junctions; (c) Hazen-Williams headloss formula with coefficients from 110 to 199; (d) diameters from 0.2032 to 2.5146 m; (e) total pipe length: 65,478 m, and (f) total base demand: 0.192 m³/s.

This network is adapted for simulating chlorine distribution in a scenario proposed by [27]. Reaction coefficients kb and kw are −1.35 d⁻¹ and −1.05 m/d, respectively. The total time simulation was extended to 240 h to balance tank reactions.

The results of the scenarios EX1, EX2, EX3 are compared with the complete mixing model of Epanet and with the Epanet-BAM using a mixing parameter of 0.5, which is the best prediction for the experiments according to [17]. The Net3.net scenario is compared only with the complete mixing model of Epanet since its configuration and control rules have a complexity that generates a different hydraulic scenario on Epanet-BAM, which does not allow comparing the concentrations results.

3. Results

3.1. OUT Coefficient and Its Concentrations Results

All the scenarios were executed in a particular CJ of an experimental network of the University of Guanajuato, of Mexico described in [19], in order to compare with the CFD simulation results.

Table 5. Results of experimental data (EXP), CFD simulations, and the relative error between them. For N and W inlets, the pressure (P) is shown and for E and S outlets, the velocity (v) is shown.

		S1			S2			S3			S4
		EXP	CFD	error	EXP	CFD	error	EXP	CFD	error	EXP	CFD	error
P (bar)	N	1.600	1.577	1.437%	0.900	0.913	1.444%	1.920	1.923	0.158%	1.880	1.893	0.670%
	W	1.600	1.576	1.478%	0.900	0.913	1.403%	1.910	1.923	0.703%	1.880	1.893	0.677%
v (m/s)	E	1.062	1.060	0.211%	1.025	1.019	0.627%	1.201	1.193	0.665%	1.005	0.987	1.720%
	S	1.323	1.322	0.133%	0.959	0.960	0.035%	0.963	0.967	0.472%	1.167	1.180	1.130%
		S5			S6			S7			S8
		EXP	CFD	error	EXP	CFD	error	EXP	CFD	error	EXP	CFD	error
P (bar)	N	1.600	1.616	1.005%	1.720	1.732	0.698%	1.570	1.578	0.502%	0.590	0.600	1.695%
	W	1.600	1.616	0.989%	1.720	1.734	0.811%	1.560	1.575	0.982%	0.590	0.600	1.695%
v (m/s)	E	1.072	1.079	0.592%	1.107	1.116	0.845%	1.040	1.044	0.381%	1.136	1.130	0.524%
	S	1.330	1.319	0.763%	1.122	1.108	1.249%	1.302	1.294	0.593%	2.263	2.264	0.026%
		S9			S10			S11			S12
		EXP	CFD	error	EXP	CFD	error	EXP	CFD	error	EXP	CFD	error
P (bar)	N	1.600	1.618	1.099%	0.660	0.649	0.879%	1.960	1.971	0.575%	0.320	0.319	1.483%
	W	1.600	1.617	1.048%	0.660	0.653	0.889%	1.960	1.971	0.557%	0.320	0.317	2.129%
v (m/s)	E	0.994	0.994	0.010%	1.080	1.081	0.141%	1.246	1.237	0.744%	1.974	1.981	0.360%
	S	1.331	1.327	0.279%	2.477	2.470	0.292%	0.427	0.421	1.192%	1.798	1.785	0.734%

shows the relative error percentages between the experimental tests and CFD simulations. Velocity errors have a range from 0.01% to 1.72%, while the error range relative to pressure is from 0.158% to 2.129%. These errors are sufficient to consider that the CFD simulations approach the real experimental tests with accuracy. Once the model is reliable from a hydrodynamic point of view, it can be confident that the velocity vectors adequately simulate the transport of diluted substances passing through the CJ.

An example surface graph with the velocity vector field of scenario S5 is shown in Figure 5. It can be verified that the vectors move in a collision-course at the inlet flows N and W. This causes an instantaneous contact in a reduced fraction of the flows, but the occurrence of a perfect mixing in this flow contact is not possible. Similar observations were discovered in [6,7]. There are areas with low velocity and even recirculation of the flow just at the beginning of the outlet pipes E and S. This is also due to the thrusts generated by the collision of the flows at the cross.

Figure 6 shows the concentrations in color surface plots through CJ for the operational conditions of scenario S10 (the highest flow shove from the N inlet to the S outlet) and all nine concentration assignations from Table 2. Due to the flow rate at N being greater than the one at W, the inlet flow impact causes a predominating outlet flow rate at the S outlet (Table 2 and Table 4). In the first plot (IN = 0.00), concentration C_S decreases rapidly due to the null concentration of C_N. The C_E = 0.113 is obtained due to the same reason. It must be noted that the greatest increase in concentration is C_S in each of the nine cases of Table 3 for scenario S10. Nevertheless, the concentration C_E is almost conserved in all cases, but there is always a concentration contribution from C_W. This fact is mentioned in the conclusions of [13,18]. Therefore, complete mixing is not attained in any case, even if the IN = 1.00 (all boundaries with the same concentrations still do not mean that a perfect mixing has occurred).

A graph of IN and OUT coefficients is shown in Figure 7. All of them show an upward trend, but none resemble the one of complete mixing. If this were to occur, the OUT coefficient would always have the value of 1. It is possible to get this result, but it is not due to the mixing. In fact, it is possible when both concentration inlets, N and W, are practically the same.

The Polynomial Equations

Each curve from the twelve scenarios in Figure 7 was approximated by a Polynomial Equation (PEi) for every scenario, as a function of the IN coefficient (set of Equations (9)), which adjusts the correlation coefficient R2 to be practically 1.

$$\begin{gathered} \mathit{PE}\mathbf{1}\left(IN\right)=-0.052841I{N}^2+1.0049IN+0.035921 \\ \mathit{PE}\mathbf{2}\left(IN\right)=-0.01706I{N}^2+0.90467IN+0.11139 \\ \mathit{PE}\mathbf{3}\left(IN\right)=-0.10307I{N}^2+1.0834IN+0.018478 \\ \mathit{PE}\mathbf{4}\left(IN\right)=-0.025456I{N}^2+0.94364IN+0.082754 \\ \mathit{PE}\mathbf{5}\left(IN\right)=0.007894I{N}^3-0.10606I{N}^2+1.0727IN+0.025038 \\ \mathit{PE}\mathbf{6}\left(IN\right)=-0.02166I{N}^3+0.048757I{N}^2+0.71965IN+0.27756 \\ \mathit{PE}\mathbf{7}\left(IN\right)=0.00048463I{N}^3-0.0087194I{N}^2+0.7664IN+0.24077 \\ \mathit{PE}\mathbf{8}\left(IN\right)=-0.08936I{N}^2+1.0723IN+0.016051 \\ \mathit{PE}\mathbf{9}\left(IN\right)=-0.068973I{N}^2+1.0368IN+0.030728 \\ \mathit{PE}\mathbf{10}\left(IN\right)=-0.077967I{N}^4+0.46528I{N}^3-1.1642I{N}^2+1.7759IN+0.00083556 \\ \mathit{PE}\mathbf{11}\left(IN\right)=-0.0012052I{N}^2+0.54835IN+0.4556 \\ \mathit{PE}\mathbf{12}\left(IN\right)=-0.0094417I{N}^2+0.85981IN+0.14881 \end{gathered}$$

(9)

With these PEi, the OUT coefficient can be analytically estimated in order to obtain the concentrations on the outlets boundaries by the methodology described in Section 2.3.

3.2. Results of the Application Examples

3.2.1. Results from the Tests in the Proposed Network

EX1 is the scenario that represents a contaminant event, with a chemical entering the WDN by the RV1. The range of flow rates on the cross-pipes varies from 0.37 to 16.60 l/s. Applying the steps to calculate the incomplete mixing (Section 2.3), relative differences greater than 90% in respect to the complete mixing were obtained. The Epanet-BAM results were closer, about 44% of the mean difference. In the complete model, the concentration remained more homogeneous in the network. In the Epanet-BAM and IMM by polynomial equations, there was a greater concentration in the southwest area, while the east area only retains a minimal concentration (Figure 8).

Some extreme results on this EX1 network were registered. For example, the first CJ, below RV2 (marked in Figure 8), has 0.00 and 1.50 mg/l for C_N and C_W inlets, respectively. The complete mixing of Epanet assumes that each outlet C_S and C_E should leave 0.40 mg/l. Applying the IMM by polynomial equations, the concentrations outlet are 0.8371 and 0.0007 mg/l for the C_S and C_E, respectively. This can happen since there is a very high difference between the inlet flows at Q_N and Q_W, about 16.60 and 6.05 l/s, respectively. The differentiated flow-impact generates a minimum concentration contribution from the W inlet (1.50 mg/l) and in the same way, the North concentration (0.00 mg/l) is conserved in a higher proportion since it has the predominant flow rate (Q_N is near 300% in respect to the Q_W). These results are consistent with [14,16,18]. The Epanet-BAM and the IMM by polynomial equations result in similar contour plots, which in most of the nodes the concentrations are above 0.8 mg/l. However, the east sector has an even lower concentration in the IMM by polynomial equations. EX1 has a sector where several nodes do not reach the value of 0.20 mg/l (this is the minimum concentration in Epanet-BAM). This shows that in contamination events, the IMM by polynomial equations could obtain more severe scenarios than the complete mixing model.

The Root Mean Square Error (RMSE) of the concentrations of all junctions in the complete and IMM by polynomial equations was 0.327 for the EX1. The average relative errors in EX1 is 56.3%, and the maximum change in concentrations in EX1 is up to 94% in comparison to the complete mixing model.

Scenarios EX2 and EX3 represent the disinfection on WDN considering the limits of chlorine concentrations in real WDN. The concentration results on the nodes have a notorious similarity for the three mixing models (complete, IMM by polynomial equations, and Epanet-BAM) (Figure 9). The average relative errors between the three mixing models in both scenarios did not exceed 17.4%. This occurs since the reservoirs are in opposite areas, which causes the flow path from both reservoirs to coincide in several junctions, generating a more balanced mixing and, therefore, the concentrations become more uniform. Unlike the EX1, the sources are relatively close, which causes the water to flow through similar paths towards the lower zones of the network and does not produce a considerable mixing as in EX2 and EX3 scenarios. However, in disinfection simulations, the importance of the IMM by polynomial equations with opposite points of the disinfectant injection will only be in specific nodes, as in this case, the maximum difference in the concentrations was around 31% in two nodes for the EX2 compared with the complete mixing model.

3.2.2. Results by the Epanet Example Network NET3.net

To describe the results of the IMM by polynomial equations, the contour plot of concentration for the simulation hour was selected with the lowest pressure in the network (Figure 10). The zones near the CJs are more influenced by the incomplete mixing, but the farthest junctions receive similar concentration rates from their incoming pipes. However, various sectors have more concentrations from the flows passing through CJs. Figure 10 shows the differences between the three mixing models. The complete mixing (Epanet) and the IMM by polynomial equations have a good agreement with their contour plot (Figure 10). Even so, it can be observed that the zones close to the CJs have different concentration ranges between both mixing models. Therefore, the proposed model could define the areas where the disinfectant is really sufficient and, with this, know its availability and distribution to the rest of the network.

The RMSE of the concentrations of all junctions in the complete and IMM by polynomial equations was 0.1577 for the NET3.net (Figure 11.) However, some nodes achieved absolute errors close to 82% (Figure 11, nodes 38 and 55). These nodes are in central areas of the network, close to the pipelines with high flow rates (429 to 445 l/s). In sets of nodes, for example, several nodes, from ID 10 to 20 and from ID 50 to 60 have absolute errors of 28 and 31%, respectively and also in areas close to the CJs.

3.3. Model Capabilities and Limitations

It is recommended that the flow combination of the CJ to be studied does not have extreme differences between the inlet and outlet boundaries. The initial calculations of this model are based on the coefficients Q_rIN and Q_rOUT, the proportion should be as close as possible to the range: 0.4 to 2.9 (similar Q_r values from Table 4). Values far from these limits do not show a sufficient similarity to scenarios S1 to S12 and the results could not be appropriate.

Another consideration is that the polynomial equations were generated with an IN coefficient between 0.00 to 2.00. The concentration values that exceed twice from one inlet to the other will be evaluated in a suitably quasilinear polynomial equation but could be outside the scope of this research.

This model is suitable for application in the CJ of different combinations of diameters. In all the examples of networks used for the application of this IMM by polynomial equations, CJs of configurations of different diameters were found. Therefore, the response of the model to these changes in geometry was evaluated for crosses from 0.0762 × 0.0508 m² to 0.762 × 0.3048 m², building them in CFD, maintaining the hydraulic structure and numerical parameters for their accurate representation.

Table 6. Outlets Concentrations at the East (C_E) and South (C_S) on Cross-Junctions (CJs) from the Incomplete Mixing Model (IMM) by polynomial equations and simulation in CFD. Percentages are absolute error values.

		EX1			EX2			EX3				Net3.net
C Outlet	Cross Diameter (m)	INC	CFD	Error	INC	CFD	Error	INC	CFD	Error	Cross Diameter (m)	INC	CFD	Error
CE	0.1524 × 0.1016	0.810	0.836	2.6%	0.902	0.868	3.3%	0.739	0.749	0.9%	0.762 × 0.3048	1.527	1.441	8.6%
CS		1.150	1.138	1.2%	0.980	0.978	0.2%	0.838	0.833	0.6%		1.434	1.416	1.8%
CE	0.1016 × 0.1016	0.060	0.040	2.0%	0.364	0.368	0.5%	0.546	0.539	0.7%	0.3048 × 0.3048	0.868	0.749	11.8%
CS		0.720	0.753	3.3%	0.643	0.600	4.3%	0.513	0.520	0.6%		1.119	1.109	1.0%
CE	0.1016 × 0.0762	0.070	0.060	1.0%	0.528	0.528	0.0%	0.397	0.394	0.3%	0.3048 × 0.2032	1.068	1.042	2.6%
CS		0.490	0.499	0.9%	0.509	0.508	0.1%	0.386	0.387	0.1%		1.195	1.204	0.8%
CE	0.0762 × 0.0762	0.890	0.876	1.4%	0.356	0.348	0.8%	0.519	0.527	0.8%	0.3048 × 0.2032	1.069	0.993	7.7%
CS		1.140	1.143	0.3%	0.458	0.480	2.2%	0.579	0.576	0.3%		1.087	1.082	0.6%
CE	0.0762 × 0.0508	1.140	1.138	0.2%	0.307	0.301	0.5%	0.509	0.511	0.2%	0.254 × 0.2032	0.554	0.523	3.1%
CS		1.150	1.142	0.8%	0.302	0.305	0.3%	0.533	0.532	0.1%		0.597	0.635	3.8%

shows a comparison of the outlet concentrations C_E and C_S on CJs applying the IMM by polynomial equations and CFD simulations for the four scenarios.

Table 6 shows 40 outlet boundaries, the errors vary from 0.1 to 11.8%. However, 72.5% of the outlets have an error below 2%, 17.5% of the outlets have an error between 2.2 to 4%, and 10% of the outlets have an error major of 4%. Therefore, the IMM by polynomial equations has an efficient representation for different configurations of cross dimeters, which are commonly used in WDN.

4. Discussion

The complete mixing model shows that the cross-leaving concentrations are of about equal magnitude, thus the OUT coefficient should have a value of 1.00 in any case. Results of this research show that the homogeneous mixing was not attained in any of the twelve scenarios analyzed.

The S10 and S11 scenarios could specially give a better reference in the convective solute transport through CJs. Their Q_rIN,Si and Q_rOUT,Si ratios have the most varied values in Table 4. Scenarios S10 and S11 have outlet flows in reverse proportion too. In other words, scenario S10 has the highest flow leaving through the S outlet, while in scenario S11, the highest flow rate is through the E outlet.

In the results depicted in Figure 6 for scenario S10, when IN = 0.00, (C_N = 0 and C_W = 250 mol/m³), the OUT value is approximately 0.001. This means that the Q_N flow shows the concentration from Q_W almost completely towards outlet S (135.17 mol/m³). Due to this, the concentration from C_W was greatly reduced, since the Q_N inlet contributes a lot of flow with no concentration. This also makes the E outlet have a very low concentration (0.110 mol/m³).

On the other hand, the S11 scenario has the highest outflow at the E boundary. A convective shove makes an important concentration contribution from both C_N and C_W inlets to the E outlet. The curves in Figure 7 for scenarios S10 and S11 should have different behaviors due to this. The trace of scenario S11 has the lowest slope, thus the trace of scenario S10 should be the highest slope. However, the behavior of the S10 trace took a decaying curved shape for values of IN>1. This happened due to the way the OUT coefficient is calculated. Emphasizing the first subfigures in Figure 6, for IN values less than 1, since E is not receiving a sufficient amount of W, the C_E concentration has values close to the adjacent inlet C_N (this also happens in [14,19] and is an important property of Epanet-BAM [15]), but when IN>1, both E and S outlets receive a lot of C_N concentrate due to C_N being greater than C_W.

The network in the EX1 example, representing a case of contamination, shows a large difference between the three mixing models. The larger diameter crossings provided greater flow towards the northwest and east areas of the network, while flows whose concentration is relatively zero (through RV2) continued to flow through the east area of the network. With complete mixing, the concentration in all the nodes is above 0.4 mg/l, and in Epanet-BAM, it is 0.2 mg/l, while in the incomplete mixing, the minimum concentration was 0.0007 mg/l. However, the contour concentration on the results of the IMM by polynomial equations has a greater area with a high concentration, which implies a greater quantity of nodes. The model shows a greater effect on the network in case of contamination.

EX2 and EX3 scenarios were closer to the complete mixing, possibly since the RV were placed in opposite areas of the network (this is a similar case to [16]) and this balanced the flow rates and allowed a more homogeneous concentration distribution. These cases are similar to the simulation of disinfection in WDN, that the IMM by polynomial equations has an effect on specific nodes.

Diversity in internal geometry and diameter changes in crosses was analyzed to be considered in the water quality simulation. The IMM by polynomial equations is the result of an experimental model with a cross of 0.1016 m in diameter, commonly used in WDN, also under commonly operated hydraulic conditions. The results are considered to be highly applicable on Epanet, according to the comparative mode in Section 3.3, where various concentration outlets for different crosses were validated and results with errors less than 4% were obtained in 90% of the cases.

The Epanet-BAM helped validate the IMM polynomial equations since it shows a range between the complete and incomplete mixings. However, the NET3.net scenario could not be compared since the configuration and controls have a complexity that generates a different hydraulic scenario on Epanet-BAM than the simulation in Epanet.

An important consequence of incorporating the incomplete mixing in water-quality simulation programs is the improvement of the calculation of temporary solute distribution through the network. With complete mixing, the concentrations are distributed in a more homogeneous way, thus the network can reach an average state at the end of the simulation. On the other hand, with the IMM by polynomial equations and depending on the location of the substance sources and the number of CJs, some areas will appear with concentrations in different proportions to the ones with complete mixing and could be a determinant mainly for simulations of contaminants than those of disinfection.

5. Conclusions

Twelve scenarios were simulated experimentally to validate a hydraulic CFD model of a cross with two contiguous inlets and outlets. Flow velocity conditions close to 2.5 m/s and pressures of 1.96 bar were reached. It is common for these conditions to occur in real WDN.

Concentration combinations were made to the inlets using the IN coefficient, covering a total of 108 cases of mixing to obtain a model for IMM by polynomial equations. In all of them, a different concentration can be seen at each outlet of the CJ. One of the factors that influence the disparity of outlet concentrations is the flow impact from both inlets of the cross. This causes an incomplete mixing due to the minimal flow contact. The outlet concentrations are due to the transport of solutes governed mainly by convection. The calculation of these concentrations for unexplored mixing cases maintains a mass-conservation derived from Equation (4).

To apply this IMM by polynomial equations directly to Epanet, a programming code was developed that processes hydrodynamics information of the network and identifies the mixing cases in CJs. The code inactivates complete mixing by reaction coefficients implemented in auxiliary nodes and bypass pipes in the CJs. The incomplete mixing is controlled by setpoint boosters located at the cross boundaries. With this, it is possible to continue analyzing the water quality with Epanet itself, considering the effects of incomplete mixing.

It is of paramount importance to consider case EX1, since there are different applications of the simulation of water quality in which it can find works that focus on avoiding situations of contamination induced through the water of the distribution networks. Pollutants can pass through some crossings and their concentration can cover different areas and affect large numbers of consumers.

The NET3.net network has a more complex hydraulic development, where the operation and direction of the flows at the cross boundaries varies at different time steps during the simulation. In this example, slight differences appeared in most nodes, as well as significant changes in the areas close to the CJs. The average RMSE value of 0.1577 represents the concentration variability that the network would have if this mixing model is considered. The absolute errors for a pair of node groups (Figure 11) ranged from 32 to 82%. These variations could improve the prediction of chlorine distribution in the network and make it possible to identify areas where the availability of chlorine concentration is insufficient to not carry out an adequate disinfection.

The IMM by polynomial equations was proposed in some conditions of the relations of inlet and outlet flows and of a relation of inlet concentrations. However, the outlet concentration of diverse CJs from the four scenarios simulated on Epanet, were evaluated with the CFD model, obtaining that 90% of the outlets had an absolute error below 4%, considering that the IMM by polynomial equations has an efficient representation for different configurations of cross commonly used in WDN.

Researchers have dedicated their studies to strategically position early-detection sensors to ensure the least harm to consumers. With a more accurate quality model, the best sensor positions can be found with better probability, and the water-quality simulation software makes better predictions on the distribution of a contaminant. Having a mixing model with greater precision will increase the reliability of the water-quality simulation software.