Assessment of Different Pressure Drop-Flow Rate Equations in a Pressurized Porous Media Filter for Irrigation Systems

: The small open area available at the slots of underdrains in pressurized granular bed ﬁlters for drip irrigation implies: (1) the existence of a region with non-uniform ﬂow, and (2) local values of modiﬁed particle Reynolds number >500. These ﬂow conditions may disagree with those accepted as valid for common pressure drop-ﬂow rate correlations proposed for packed beds. Here, we carried out detailed computational ﬂuid dynamics (CFD) simulations of a laboratory ﬁlter to analyze the results obtained with ﬁve different equations of head losses in porous media: (1) Ergun, (2) Darcy-Forchheimer, (3) Darcy, (4) Kozeny-Carman and (5) power function. Simulations were compared with experimental data at different superﬁcial velocities obtained from previous studies. Results for two silica sand media indicated that all equations predicted total ﬁlter pressure drop values within the experimental uncertainty range when superﬁcial velocities <38.3 m h − 1 . At higher ﬂow rates, Ergun equation approximated the best to the observed results for silica sand media, being the expression recommended. A simple analytical model of the pressure drop along ﬂow streamlines that matched CFD simulation results was developed.


Introduction
World annual freshwater extraction for municipal, industrial, and agricultural needs is approximately 3928 km 3 y −1 [1]. Some 44% of this water (1716 km 3 y −1 ) is consumed, mainly for irrigating farmland (38% of freshwater extraction), and the remainder (56%, 2212 km 3 y −1 ) is primarily released to environment in the form of wastewater, industrial effluents or agricultural drainage water [1]. The expected change in future rainfall scenarios as well as the increase of water demand for socio-economical needs require the implementation of adaptation measures in some areas, as in the Mediterranean Basin [2].
Water reuse from wastewater treatment plants appears as a suitable strategy to mitigate the stress on natural resources, with a water quality level depending on the end-user types [2]. Reclaimed water for crop irrigation must meet a minimum of quality standards to avoid any health hazard, but the cost of additional treatments and wastewater discharge taxes may compromise its financial viability [3]. However, drip irrigation methods are very efficient in water consumption terms [4] and avoid the contact between plant leaves or fruits and wastewater, which is advisable to further reduce contamination risks in some crops [5]. Therefore, drip irrigation with reclaimed water might be an appropriate method to cope with forthcoming water scarcity issues.
One of the main problems of the drip irrigation technique is the emitter clogging that raises maintenance costs and reduces crop productivity [5]. The filtration process of treated water may effectively prevent drippers to clog, thereby becoming essential for sustainable irrigation with non-conventional water resources [6]. The removal of particles and microorganisms and the decrease in turbidity may be satisfactorily achieved with pressurized porous media filters working under deep bed filtration conditions [7]. These porous media filters use a given amount of granular material like sand or crushed glass so as to create a bed inside a pressurized tank [8]. Underdrain elements with slots smaller in width than the grain size are located at the base of the packed bed to avoid loss of granular material [9]. In filtration mode, pumping work is needed to achieve a nominal flow rate, being proportional to the pressure drop through the filter [10]. Total filter pressure drop increases as particle retention in the packed bed progresses, eventually reaching a threshold value beyond which the backwashing regime begins [11]. In the latter mode, pumping work more intense than that of filtration is required to develop fluidization of the granular bed assuring particle removal [9]. As a consequence, an entire filtration cycle involves energy consumption that becomes a setback to drip irrigation implementation. The reduction of filter head losses would imply longer filtration mode runs and, hence, improve energy efficiency.
Filter pressure drop is the addition of purely hydraulic head losses (caused by major flow friction losses and minor losses from auxiliary elements such as diffuser, underdrains, etc.) plus the energy losses through the porous media [12]. Redesigns of current filters to enhance energy efficiency require models that correctly capture the flow behavior in all regions and, especially, in the packed bed. Erdim et al. [13] carried out a comprehensive study of the prediction capacity of 38 pressure drop-flow rate correlations found in the literature that applied to porous media formed by spheres of single diameter. The range of applicability of those expressions varied depending on the modified particle Reynolds number Re ε , with many of them sharing Re ε intervals. Erdim et al. [13] experimentally found that the well-known Ergun equation failed to predict the energy losses at Re ε > 500, and proposed a new correlation that better fitted data.
A relevant point related with many of those previous expressions was that the pressure drop per unit length was expected to be valid for a region with uniform flow within the packed bed. However, this is not always accomplished in commercial filters since the effect of having underdrain slots with an open area lower than that of the crosssectional filter implies a flow concentration towards the drainage [14]. Arbat et al. [12] and Pujol et al. [15] aimed to solve this problem by including non-uniform flow corrections to analytical pressure drop equations within the granular packed bed. The modifications of the equations took into account the variation in the number of water channels available when approaching the underdrain element. However, those expressions were not validated locally but, once added to the rest of the momentum sinks, compared globally with measured data of the total filter pressure drop.
Besides analytical expressions, several authors have applied numerical methods to simulate the flow behavior of pressurized porous media filters. The computational fluid dynamics (CFD) technique based on the finite volume method uses a momentum sink to determine the pressure drop per unit length within the packed bed [16]. This methodology has allowed the analysis of both laboratory and commercial filter type designs. Bové et al. [17] employed the CFD technique to corroborate the superior performance of a modified underdrain that reduced the pressure drop mainly by increasing the zone with uniform flow within the granular bed. The same authors used simulations to discuss the reasons of having different data values of the total filter pressure drop of two underdrain designs [18]. Mesquita et al. [19] confirmed by means of CFD analyses of three different types of underdrains that the best hydraulic behavior corresponded to the design with the highest flow uniformity within the porous media. More recently, Mesquita et al. [20] modeled different designs of a diffuser plate, finally proposing one that increased the flow uniformity above the sand surface. The validity of this design was experimentally confirmed [20]. Also based on conclusions from CFD analysis, Pujol et al. [21] proposed a redistribution of underdrains in a commercial sand media filter that reduced the total filter pressure drop by 5.8%.
However, a drawback of previous CFD studies is that they relied on expressions to estimate the momentum sink within the packed bed extracted from experimentally adjusted data that, in some cases, were not totally independent. For example, [17,18] assumed a pressure drop per unit length in the porous region that followed a quadratic function in terms of superficial velocity and whose coefficients were determined from pressure measurements in a slab within the granular bed under those conditions simulated. Analyses found in [14,21] also calibrated the momentum sink coefficients of the pressure drop equation employed in the CFD for the porous media by minimizing the error with respect to experimental data. On the other hand, though numerical analyses in [14,[17][18][19][20][21][22] were correctly set up for the objectives proposed, the level of discretization applied in the numerical models was somehow limited to include the recommended growth of layers of prismatic elements at certain walls only, ignoring those small walls of the underdrain element.
From the above, several questions may arise: can common pressure drop-flow rate correlations for packed beds correctly predict the total filter pressure drop even though there exist regions near the underdrain with high Re ε values? To obtain a successful CFD simulation, is the common procedure of fitting coefficients of empirical pressure drop equations for the porous media a consequence of not using a massive mesh? Is there a physically based equation with coefficients a function of flow and granular media properties that can accurately describe the pressure drop within the entire packed bed? Can we derive a simple analytical (i.e., closed-form) expression to determine the pressure drop in the non-uniform flow region of the packed bed? The present work was intended to answer the previous questions by analyzing the results of a numerical model of a pressurized packed bed filter that used a commercial type underdrain and comparing them with experimental data.

Experimental Data
Experimental data used in the present study was obtained in [23], who analyzed how different types of granular media modified the pressure drop of a pressurized filter. Information of the experimental set up can be found in [23] and it is not reproduced here since new data was not collected. However, the design of the laboratory filter is introduced in detail as it was the basis of our numerical analysis. The dimensions of the laboratory filter can be found in Figure 1. It consisted of a steel cylinder of height 750 mm and inner diameter D f = 200 mm (filter cross-sectional area A f = πD 2 f /4 = 31,416 mm 2 ). Vertical inlet and outlet steel pipes of 32 mm inner diameter were welded to both upper and bottom cylinder covers. An inner plate, 10 mm thick, was welded at a distance 133 mm from the bottom end of the cylinder. A single commercial pod-type underdrain (Regaber, Parets del Vallès, Catalonia, Spain) was installed at the center of the inner plate. This pod was assembled from three individual components (see Figure 1). It had 45 rectangular slots of dimensions 0.45 × 30 mm distributed in a truncated conic surface and oriented in vertical planes. The total slot open area of the underdrain was A i,u = 607.5 mm 2 . The inner diameter of the underdrain conduit that discharged into the bottom water-only chamber was 16 mm, with a total outlet area of the underdrain A o,u = 201.1 mm 2 . The diffuser plate was a circular element of 103 mm diameter located 105 mm below the cylinder top. It was supported by three rods welded to the inner side of the top cover.
Pressure data was measured at five points with digital manometers Leo 2 (Keller, Winterthur, Switzerland). Two manometers were installed at both inlet and outlet pipes (at 51.5 mm from covers, see pI and pO locations in Figure 1) and three manometers with 100 mm vertical gap at the side of the filter main body, with the first being located 100 mm above the underdrain base plate (see p1, p2 and p3 locations in Figure 1). Laboratory experiments were carried out with different porous media filled up to a depth equal to Water 2021, 13, 2179 4 of 25 117 mm or 317 mm above the base plate [23]. The 117 mm height was chosen for analyzing the performance of shallow filtration depths, since some studies have been carried out with small media depths [6]. Thus, from 1 to 3 pressure data values in the porous media bed were recorded. Experimental volumetric flow rates values Q varied from 0.092 L s −1 to 0.799 L s −1 , corresponding to superficial velocities v s (= Q/A f ) ranging from 10.5 m h −1 to 91.6 m h −1 . Note that rapid granular bed filters work with superficial velocities around 20 m h −1 [10] though much higher velocities (i.e., 61 m h −1 ) have been suggested [6].

Porous Media
The granular media employed were two types of silica sand and one type of microspheres [23]. The later were produced from selected flat glasses (Sovitec, Felurus, Wallonia, Belgium), whereas silica sand (Sibelco Hispania, Bilbao, Spain) was extracted from a quarry, being washed and dried at reception. Grain size ranges for the three media are listed in Table 1. These were obtained after repeatedly sieving the raw materials and their products with appropriate stainless steel screens. The average purity of the final fractions was expected to be 96.6%, varying from 93.1% to 98.1%.
Two of the pressure drop equations analyzed (Ergun and Kozeny-Carman, see Section 2.3) required the physical characteristics of the porous media. The procedures to obtain the porosity, equivalent diameter and sphericity are here summarized, shown in Table 1 and explained in detail in [23]. Porosity ε values, defined as the ratio of voids volume to bulk volume, were derived from both the real media density ρ r and the bulk density ρ b [23]: Bulk density was obtained after weighting a known volume (400 mL) of porous media. Real density was derived by calculating the volume of solids of the previous media by adding 400 mL of water and measuring the total volume achieved in a graduated cylinder. These procedures were repeated a minimum of three times and the averaged values are reported in Table 1. The equivalent diameter D eq was defined as the diameter of a sphere with a volume equal to the average volume of the media grains [24]. The procedure consisted of measuring the mass of 1000 grains. This value divided by 1000 gave an average mass per unit grain. The entire process was repeated three times to obtain the final average value of the mass per unit grain m g , from which the equivalent diameter was calculated: The sphericity coefficient ψ was defined as the ratio of the surface area of the sphere with equivalent diameter D eq to the actual surface area of the grain. The calculation of ψ was based on the indirect method of Soyer and Akgiray [24]. In essence, the method assumed that the Ergun equation correctly predicted the pressure drop ∆p in a region of depth ∆L with uniform flow. Thus, the pressure drop per unit length within the porous media followed the expression: where µ and ρ were the water absolute viscosity and water density, respectively, and v s was the superficial velocity. The sphericity value ψ was calculated as the value that maximized the Nash-Sutcliffe efficiency coefficient [23] defined as: with O i the measured data of pressure drop per unit length, and P i the predicted data from Equation (3) of the pressure drop per unit length, both at superficial velocity v s,i with i = 1, . . . , N the number of experimental points (N = 34 for SS2 and N = 36 for SS1 and MS media). Data used ∆L = 200 mm and ∆p as the pressure difference measured between heights 300 mm and 100 mm above the nozzle base plate (manometers p3 and p1 in Figure 1). In Equation (4), O is the average value of the measured data. We point out that the dataset employed to estimate ψ coincided with that used in the study of cases with porous media height equal to 317 mm above the sand base plate. Therefore, the sphericity coefficient ψ for the Ergun equation was not really independent of the dataset for both silica sand cases in the study of a porous media height equal to 317 mm, as discussed in Section 4. For glass microspheres we assumed ψ = 1 and the previous procedure was not applied. We point out that other authors have reported values of the coefficient of sphericity for silica sand (particle fractions from 0.60 mm to 0.85 mm) ranging from 0.72 [25] to 0.93 [26], which includes the value here found (=0.89). The indirect method of averaging the solution of the quadratic equation for ψ by rearranging Equation (3), as in [24], predicted the same result (±0.01) as that reported in Table 1.

Momentum Source Terms for Porous Media
The effect of the granular bed on the flow behavior was added as a source term S i to the momentum Navier-Stokes equations. This source term was, indeed, a momentum sink since the tortuosity and frictional effects of the small water channels within the packed bed produce flow energy losses. As pointed out in the introduction, several expressions for this sink term have been proposed, whose range of validity depends on the flow regime and packed bed characteristics (see the comprehensive study carried out in [13]). The dimensionless particle Reynolds number Re p : as well as the modified particle Reynolds number Re ε : are often used to assess the range of applicability of these equations [13]. For example, Ergun equation fails to predict pressure drop in porous media with Re ε > 500 [13]. Of course, such high values of Re p and Re ε are not expected in the uniform flow zone of packed bed filters for drip irrigation since superficial velocities are very modest. However, local values of Re p and Re ε may be over 500 near the underdrain slots. From Table 1 and experimental values of superficial velocity, the standard (modified) particle Reynolds number of datasets for SS1 was 2.1 < Re p < 18.6 (3.6 < Re ε < 32.1), for SS2 was 3.1 < Re p < 23.3 (5.2 < Re ε < 38.8), and for MS was 1.3 < Re p < 9.5 (2.2 < Re ε < 15.3). Flow inertial effects are expected to be of minor importance in filtration beds up to Re p about 1 [27]. In this case, the main contribution to the momentum sink S i is the viscous term, being known as Darcy equation: where v s.i is the ith component (i = x, y, z) of the velocity based on a filter empty of granular media (i.e., superficial velocity concept) and the constant α is the permeability term whose value can be obtained by fitting experimental data. An expression with the same linear relationship between pressure drop and superficial velocity (or, equivalently, flow rate) as in (7) can be derived from the Hagen-Poiseuille equation including a tortuosity correction factor [28]: being known as Kozeny-Carman equation and expected to be valid for Re ε < 2.
In scenarios with larger flow rates through the packed bed (i.e., Re p > 1), the inertial term becomes relevant and it is added to (7), now being known as Darcy-Forchheimer Equation [29]: where |v S | is the magnitude of the superficial velocity and C 2 is the inertial resistance factor. Indeed, the Ergun Equation (3) used to determine the sphericity coefficient follows a Dary-Forchheimer formulation, since it contains both linear and quadratic terms on flow velocity. The Ergun Equation (3), expressed in terms of a momentum source term for non-uniform flow, is: and it can be derived from adding the Burke-Plummer empirical correlation of pressure drop in packed beds for Re p > 1000 [27] to the Kozeny-Carman Equation (8) with the 180 constant term modified. Note that the second term on the right hand side of Equation (10) follows the expected S i vs. v 2 relationship of the pressure drop for turbulent flows. Finally, we also analyzed the power function model (see, e.g., [30]): where C o and C 1 are two constants. The value of 1/α for the Darcy Equation (7) was extracted from a linear regression fit through the origin (null pressure drop at zero flow) of experimental data for the case of bed height h s = 317 mm as a function of v s . The source term was calculated as S i = ∆p s,200 /∆L with ∆p s,200 the difference of the pressure drop between manometers located at 100 mm and 300 mm above the sand base plate and ∆L = 200 mm. Similarly, values of 1/α and C 2 for the Darcy-Forchheimer Equation (9) were derived from a parabolic regression fit through the origin. Finally, coefficients C o and C 1 of the power function Equation (11) were obtained after a linear regression fit of log (−S i ) (see Table 2 and Figure 2). Most of the power functions proposed in the literature (see [13,30]) modify the v s exponent of the Darcy-Forchheimer inertial term while maintaining the permeability one linear with v s (i.e., S i = Mv s + Nv λ s with M, N and λ constants). In these cases, λ values vary from 1.80 to 1.95 [13,30] (λ = 2 in Equation (9)). In flows where the viscous linear term is important, the exponent C 1 of a single term power function as in Equation (11) must be lower than λ, and a value of C 1 = 1.5 for Re ε > 675 has been suggested (see the discussion in [13]). Here, the C 1 value was smaller (see Table 2) since Re ε < 38.8 and inertial effects were moderate.
Note that the prediction of the Ergun equation was not independent of the experimental data shown in Figure 2 either. As explained in the previous section, these data were used to determine the sphericity coefficient of the porous media by applying the Ergun equation and the Nash-Sutcliffe efficiency coefficient.

Model Setup
The numerical model was developed with ANSYS-Fluent (Canonsburg, Pennsylvania, USA) [16]. This commercial software has successfully been applied in the analysis of pressurized porous media filters for drip irrigation (see [12,14,17,18,21,22]). This CFD code uses the finite-volume technique to solve the flow continuity and momentum equations, which included the sink term of the previous Section applied to the packed bed region. To save computational resources, the simulation domain consisted of half the filter (as in Figure 1a) since the model had symmetry along a vertical centered plane. The simulation domain was divided into three zones: water only upper region (from filter inlet to the upper surface of the packed bed), porous media region (packed bed; momentum sink term on), and water only lower region (from underdrain slots to filter exit). Table 2. Values of the inverse of the permeability α, inertial resistance factor C 2 , and power function coefficients C 0 and C 1 fitted with experimental data for Darcy, Darcy-Forchheimer and power function equations.  The finite-volume technique required the discretization of the domain into small control volumes named elements that formed the mesh. Since dimensions of interest varied several orders of magnitude (e.g., from filter diameter of 200 mm to slot width of 0.45 mm), an unstructured mesh carried out with ANSYS-Meshing software (Canonsburg, PA, USA) was chosen to discretize the three regions. A very refined mesh was used to study the flow behavior at the underdrain zone. In all regions, triangular elements were adopted to mesh surfaces and tetrahedrons to mesh the volumes except near all walls, in which five layers of triangular prisms were defined.

Porous Media Darcy (Equation (7)) Darcy-Forchheimer (Equation (9)) Power Function (Equation (11))
The maximum size of the elements at the underdrain slots was 0.1 mm only, with elements of maximum size of 0.2 mm in the underdrain side walls and 0.3 mm in the underdrain top and bottom faces (Figure 3). The maximum size of volume elements within the underdrain was 0.6 mm, whereas it was 5 mm in other regions except at both inlet, diffuser plate and outlet zones in which the maximum size of elements was 2 mm. In the sand region, near the pod, elements had a characteristic size of 0.5 mm. A total amount of 16.07 × 10 6 elements was required to mesh half of the filter, though other meshes were developed to analyze the sensitivity of the results to discretization changes, as explained in the next Section. The mesh metrics parameters were 0.899 for the maximum skewness value, 43.6 for the maximum aspect ratio value and 0.10 for the minimum orthogonal quality value. In comparison with previous studies that simulated the behavior of pressurized filters for drip irrigation, our study substantially improved the mesh applied. Thus, it was the first time that all existing walls in the simulation domain had inflation layers (i.e., layers of extruded prisms from the attached surface), this being very important to correctly simulate the flow through the underdrain element. Previous studies only considered meshes with inflation layers attached to the filter outer walls ( [14] in multi pod filters) or simply ignored them ( [18], who numerically analyzed the same laboratory filter). For example, the mesh here defined had 12 or 13 elements across the width of the underdrain slots (see Figure 3), whereas it was limited to 4 only in [18]. And more important, the layer of prisms near all the inner walls of the system, with element heights as low as 0.018 mm, let to y+ values below 3.5 even in the high speed conduit at the pod exit.
Turbulence effects were considered in water only regions by adopting the shearstress tensor (SST) k − ω model, in which k was the turbulent kinetic energy and ω the turbulent specific dissipation rate. This model has been proven to correctly simulate the turbulence in a wide variety of flows [31,32], including granular media filters (e.g., [14,21]). Laminar conditions were enforced in the packed bed region. Fluid with constant density (1000 kg m −3 ) and absolute viscosity (0.001 Pa s) was used. Boundary conditions were fixed absolute pressure (=141,325 kPa) at filter outlet and fixed velocity at filter inlet, where turbulent intensity was equal to 5% and turbulent viscosity ratio equal to 10. Four different flow rates were investigated (0.2, 0.4, 0.6 and 0.8 L s −1 ) that implied superficial velocities from 22.8 m h −1 to 91.7 m h −1 . Two packed bed heights h s were simulated: h s = 117 mm and h s = 317 mm. ANSYS-Fluent introduces the sink term of the momentum equations for porous media through either Darcy-Forchheimer or porous-law Equation (11) (see [16]). Thus, input data for porous media were, besides porosity ε, either the inverse of permeability 1/α and the resistance factor C 2 or the values of the power function coefficients C o and C 1 .This meant that equivalent values of the inverse of permeability 1/α and resistance factor C 2 were derived for both Kozeny-Carman (8) and Ergun (10) Equations (see Table 3). Table 3. Equivalent values of the inverse of the permeability (=1/α) and inertial resistance factor (C 2 ) for both Kozeny-Carman and Ergun equations required to set up the porous media in ANSYS-Fluent.

Porous Media
Kozeny-Carman (Equation (8)) Ergun (Equation (10) Finally, the solution of the governing equations involved an iterative process until results converged up to a defined tolerance. At each iteration, both continuity and momentum equations were solved simultaneously by applying the coupled numerical algorithm [16]. Second order schemes were chosen for all the spatial discretization terms. The absolute value of the residuals was fixed to 10 −5 for all variables. The numerical model was initially run in single precision for a minimum of 4000 iterations under steady conditions were reached. After that, the simulation was run in double precision a minimum of 400 additional iterations. A computational time of approximately 60 h per simulation was needed in an Intel Xeon W-2155 CPU with 128 GB of RAM. Pressure and mass flow rates at both inlet and outlet plus pressure at locations corresponding to both inlet and outlet tube manometers were monitored. Reported average data was calculated from the last 150 iterations of simulation datasets.

Mesh Sensitivity Study
The discretization error was estimated by investigating meshes with different size with and without inflation layers (see Table 4). All investigated meshes were unstructured with non-uniform element size inside the computational domain. The element size was smaller near the regions expected to have non-uniform flow. Characteristic sizes employed in those regions described in Section 2.4 for Meshes #1 to #3 followed the ratio of the maximum element size at the underdrain slot with respect to Mesh #4 found in Table 4. For example, the maximum element size at the underdrain side walls for Mesh #3 was 1.5·0.2 = 0.3 mm since 1.5 (=0.15/0.10) was the Mesh#3/Mesh#4 ratio of the first row in Table 4 and 0.2 mm was the value adopted for these surfaces in Mesh #4 (see Section 2.4 above). Note that commercial filters based on pod-type underdrains include a considerable number of drainage units (e.g., 12 pods in a sand filter of 500 mm inner diameter [22]). Therefore, the level of detail accomplished with Mesh #4, which was applied in the results section, is hardly achieved when studying commercial filters due to computational limitations. In this case, grids similar to either Mesh#1 or Mesh#2 without inflation layers at the underdrain units are more likely to implement (e.g., [12,17,18,21]). Table 4. General characteristics of the meshes analyzed (with inflation layers). Similar meshes without inflation layers were also adopted. Filter pressure drop ∆p f values for these different grids with a packed bed height of h s = 317 mm, using SS2 media with the Ergun Equation (10) and a volumetric flow rate value equal to 0.4 L s −1 are shown in Figure 4. The grid convergence index GCI 21 fine = 0.3% (order of convergence p = 4.7) was computed from the results of meshes #2 to #4 (including inflation layers) with ∆p f chosen as the objective value. We found discrepancies of the ∆p f value obtained in meshes with and without inflation layers. Differences larger than 2% were observed between the prediction of Mesh #4 with inflation layer and that of Mesh #2 without inflation, the latter being a characteristic grid of CFD studies of commercial filters. This variation increased above 3.5% when analyzing granular bed heights of h s = 117 mm for the same flow conditions. Therefore, the adoption of meshes without a high level of refinement may have caused that some authors discarded the use of the Ergun equation to represent the porous media behavior in commercial filters due to the apparent inaccuracy of the simulated results and adopted the Darcy-Forchheimer expression adjusted to reproduce experimental data [14,21]. We point out that the results from Mesh #4 were considered as very robust since differences of the ∆p f value less than 0.2% were found including the surface roughness of the walls, ten layers of prisms instead of five, or changing the turbulence model to the k − ε, with ε the turbulent dissipation rate. Other turbulence models like large eddy simulation (see, e.g., [33]) were note explored.

Filter Pressure Drop
Total filter pressure drop ∆p f (= pI − pO, see Figure 1a) took into account not only head losses in the packed bed but also those in water-only regions. Simulations of ∆p f for different values of superficial velocities confirmed that both Darcy and Kozeny-Carman expressions underestimated the pressure drop in the granular region at high velocities (see Figures 5-7). Particle Reynolds number Re p at v s = 22.9 m h −1 (=0.2 L s −1 ) ranged from 4.1 (microspheres) to 5.9 (silica sand 2), which were above the expected limit (Re p ≤ 1) assumed to ignore inertial effects in porous media equations. However, even at v s = 45.8 m h −1 (Re p varying from 8.3 to 11.7 depending on the granular media), ∆p f values from Darcy and Kozeny-Carman models were almost the same as those from Ergun, Darcy-Forchheimer and power function versions. Nevertheless, the assumption of momentum sink terms linear with superficial velocity clearly failed to reproduce experimental data at higher flow rates (see Darcy and Kozeny-Carman at v s = 68.7 m h −1 and 91.7 m h −1 in Figures 5-7).   Despite having a sink term linear with velocity, ∆p f for Darcy and Kozeny-Carman models did not behave linearly as a function of v s . Instead, a parabolic expression ∆p f vs v 2 s better fitted simulated data, especially when analyzing cases with h s = 117 m. This was due to the relevance of hydraulic losses mainly through the underdrain, in which the flow regime was clearly turbulent and, hence, the head losses were proportional to v 2 s . Parabolic fits to simulated total pressure data were enforced to go through the origin, with coefficients of determination R 2 > 0.999. Differences of these parabolic fits to simulation results with experimental data were reported in terms of root mean square error (RMSE) and Nash-Sutcliffe efficiency coefficient (NSE; Equation (4)) to assess the predictive behavior of models ( Table 5).
The numerical model that used the Ergun equation for calculating the pressure drop in the packed bed showed the best results in terms of RMSE and NSE for both SS1 and SS2 media (NSE > 0.978 and RMSE < 1.3 kPa for any packed bed height, with NSE > 0.99 in three of the four sand cases analyzed). Results with the Darcy-Forchheimer expression were consistently higher, particularly in the SS2 media in which RMSE values doubled those found in almost all cases when using Ergun. Among the linear approximations, the Darcy expression behaved better than Kozeny-Carman for all scenarios analyzed, as expected since the former used coefficients obtained after fitting observed head losses in a 200 mm slab within the packed beds. The power function expression gave better results than the linear expressions (Kozeny-Carman and Darcy), especially in the SS2 media, being superior than the Darcy-Forchheimer in terms of small values of RMSE and high values of NSE (Table 5). Results for microspheres, however, almost reversed the above conclusions, since second-order polynomial expressions (Ergun and Darcy-Forchheimer) had the worst indicators, with the power function results being the most similar to data. We tackle this issue in the discussion section.

Pressure Drop Per Filter Zones
As pointed out in Section 2.4, the porous media filter from the hydraulic point of view may be understood as formed by three main zones: (1) upper water-only domain from the filter inlet to the free surface of the packed bed; (2) porous media domain; and (3) lower water-only domain from the underdrain inlet to the filter outlet. The contribution of these zones to the total filter pressure drop ∆p f greatly varied as can be seen in Figure 8 for the case of SS1 media with h s = 317 mm at flow rate equal to 0.8 L s −1 . We point out that although the contribution ∆p/∆p f changed between models, head loss values ∆p for water-only regions were almost identical for a fixed value of Q. For example, at flow rate equal to Q = 0.8 L s −1 (scenario analyzed in Figure 8), we found ∆p = 17.2 kPa from the underdrain inlet to the filter exit (cyan in Figure 8) and ∆p = 0.04 kPa from the filter inlet to the top of the packed bed (blue in Figure 8, not visible due to its small contribution) for all porous media types (SS1, SS2 or MS) and porous media Equations (7)- (11). This meant that the hydraulic behavior of water-only regions were almost independent on the treatment of the porous media, being only a function of the volumetric flow rate Q. Pressure drop values in the water only region from pod inlet to filter exit were ∆p = 1.1 kPa, 4.3 kPa and 9.7 kPa for Q = 0.2 L s −1 , 0.4 L s −1 and 0.6 L s −1 , respectively. Since the flow was turbulent in this region (flow Reynolds number Re > 2500), hydraulic head losses were proportional to v 2 s or, equivalently, to Q 2 , as it can readily be inferred from the previous values of ∆p at Q.
In Figure 8, the analysis of the contribution of the porous media to the filter pressure drop was subdivided into three zones: (1) the upper zone from the free surface of the packed bed (h s = 317 mm) to z = 117 mm (pale orange in Figure 8), (2) from z = 117 mm to z = 57 mm (orange in Figure 8), and (3) from z = 57 mm to pod inlet (red in Figure 8). Region 1 was null in scenarios with packed bed height h s = 117 mm. At the first 200 mm within the granular media (from 117 < z < 317 mm), the flow behaved uniformly, since simulated pressure drop values differed less than 0.3% from those calculated by using analytical Equations (7)-(10) with superficial velocity v s = Q/A f and ∆L = 200 mm. However, discrepancies with uniform flow predictions were observed in region 2, from z = 117 mm to z = 57 mm, since all simulations predicted head losses 10% higher than those of Equations (7)-(10) with constant superficial velocity v s = Q/A f for a slab of 60 mm height. This was a consequence of the flow convergence towards the open area of the underdrain, as later discussed from streamlines data. This effect was much more pronounced in the third region analyzed (from z = 57 mm to pod inlet). This phenomenon was more intense for those models that included the inertial effect (Ergun and Darcy-Forchheimer). In these cases, the contribution of this zone to the head losses in the porous media was higher than that observed in the whole 260 mm thickness slab above. In Figure 8, for example, the small region from z = 57 mm to pod inlet accounted for 38% of the total filter pressure drop when adopting the Darcy-Forchheimer equation, reducing to 17% in case of using the Darcy one (viscous term only). Since the flow was uniform in the first 200 mm within the packed bed, ∆p values in regions from z = 117 mm to z = 57 mm and from z = 57 mm to pod inlet obtained from simulations with h s = 117 mm were the same as those with h s = 317 mm. The contribution of the bottom water-only region to the total pressure drop was also very remarkable for all cases (cyan in Figure 8). This was caused by the hydraulic loss through the complex inner parts of the underdrain. Therefore, this pressure drop reduced as a function of the square of the volumetric flow rate. For example, in the Darcy model for h s = 117 mm, the bottom water-only region accounted for 64% of the total filter pressure drop at Q = 0.8 L s −1 , whereas it was 47% at Q = 0.4 L s −1 .

Averaged Pressure Vertical Profile
Instead of analyzing pressure vertical profiles along vertical lines (as in [17,18]), we reported values of differences with respect to the inlet value of the averaged pressure on the horizontal plane XY at height intervals z of 10 mm (Figure 9). Averaged values within the granular bed ignored the contribution within the underdrain so as to better take into account the media effect. The vertical profile through the packed bed had a constant slope in the uniform flow region, being maximum for the Darcy-Forchheimer formulation and minimum for the Kozeny-Carman one. Averaged pressure values smoothed the model response to changes in flow velocity within the sand, leading to almost constant slopes up to z = 40 mm below which the trend was clearly non-linear. This apparently indicated small variations of pressure values between pressure drop models in porous media near the underdrain region (see, e.g., the variation obtained at z = 10 mm in Figure 9). However, pressure values at the underdrain slots substantially differed between models. Indeed, the pressure differences at the underdrain slots were equal to those observed at the bottom chamber (region starting at z = −10 mm in Figure 9), being invariable through all this water-only region. Vertical profile data for case h s = 117 mm almost exactly matched those for h s = 317 mm but neglecting the 117 < z < 317 mm contribution of the packed bed (see Figure 9). For SS1 media, the model with the Ergun equation better reproduced the experimental data ( Figure 9) for all values of flow rate simulated. The same was concluded for SS2 but not for MS, in which the power function data better approximated the observations.

Flow Uniformity
A descriptive analysis of flow uniformity within the porous media was carried out by means of discretizing the filter cross-sectional area into ten sectors of equal surface area. The innermost of these sectors was circular with radius r 1 such that A 1 = πr 2 1 = A f /10. The rest of the sectors were annular with radius r i that satisfied the condition A i = π r 2 i − r 2 i−1 = A f /10 for i = 2, . . ., 10. A perfect uniform flow condition would imply Q i /Q = 0.10 = 10%, for i = 1, . . ., 10, with Q i the volumetric flow rate through area A i , and Q that through A f . These ten sectors were also defined at different heights within the packed bed.
For the conditions analyzed in Figures 8 and 9b and referring to simulations with the Ergun model only, the box plot of Q i /Q values (in %) at different heights are shown in Figure 10. Small hollow squares indicate the 10% position (mean value of all datasets). At level z = 117 mm, differences between sectors were very small, pointing out a high flow uniformity at that height, as it was already confirmed in the previous analysis of Figure 8. As the surface approached the position of the underdrain, variations between annular sectors increased, though moderately till z = 70 mm. In cases with z ≥ 50 mm, the sector that contributed the less to the flow rate was the outermost annular ring A 10 . On the other hand, the sector that contributed the most at z = 30, 50 and 70 mm was the innermost annular ring A 2 since the circular central sector (A 1 ) did not have a full open area beneath it but the pod cover. Conclusions from other model equations of porous media were similar to those extracted from the Ergun case. Thus, Figure 10 explained how the flow uniformity broke down approximately below a horizontal plane at z = 70 mm above the underdrain base plane.

Flow and Pressure Behavior in the Underdrain Zone
Flow velocity vectors through the center of the slot located at an angle of 44 • with respect to the symmetry plane of Figure 1a are shown in Figure 11 (case h s = 317 mm, SS1 granular media, Ergun equation and flow rate 0.8 L s −1 ). There was a high contrast between the flow speed in the packed bed region and that found inside the underdrain. In this case, the superficial velocity in the uniform flow zone within the sand was 91.7 m h −1 (=0.025 m s −1 ) reaching a maximum value of 4740 m h −1 (=1.35 m s −1 ) at the slots, though it was only 0.10 m s −1 at a distance approximately equal to 2.8 mm from them. Inside the underdrain, we clearly observed the formation of two vortices in both upper and lower regions (see Figure 11) originated by the intensity of the water jet at the pod inlet. At the same time, there appeared secondary vortices at the beginning of the exit underdrain pipe, mainly formed by the sharp edges of the geometry. This effect implied a flow concentration near the center of the exit pipe, where the flow reached velocities up to 6.7 m s −1 .
All the previous effects involved primary and minor energy losses leading to large pressure drop values through the underdrain. Figure 11b clearly indicated the abrupt decrease in pressure (>16 kPa) along the entrance zone of the underdrain exit pipe. This value turned out to be larger than that attributed to that of the packed bed at the entry to the pod (on the order of 8 kPa). The phenomena inside the underdrain were purely hydraulics, and several improvements were previously suggested such as enlarging the diameter of the underdrain exit pipe [17] or avoiding minor losses by utilizing the entire filter cross-section as the water drainage area [34]. The pressure drop within the non-uniform zone in the granular bed, however, was not straightforward to interpret. Therefore, we defined ten streamlines in this region starting at the top of the granular bed located at radial intervals of 10 mm (see Figure 11a). The pressure drop through streamlines #4 and #8 are carefully analyzed in the next section.

Discussion
Figures 8 and 10 confirmed the departure from uniform conditions as the flow approached the underdrain element. Arbat et al. [12] and Pujol et al. [15] proposed a simplified pressure drop model to take this effect into account. In their analytical model, the granular bed was divided into two regions: region 1 far from the underdrain in which the flow was essentially uniform, and region 2 close to the underdrain with a non-uniform flow behavior. The pressure drop formulation applied to region 1 followed Ergun Equation (10) being understood as the sum of the Poiseuille friction loss term (linear with Q) as fluid flowed through the water channels in the packed bed plus a minor loss term (proportional to Q 2 ) due to sudden expansions and contractions as channels were not straight (see [27]). In region 2, both major and minor loss terms were modified according to a reduction in the number of channels available. In region 2, [12,15] assumed a linear decrease in the available diameter for the water flow from D f (=200 mm; filter diameter) to D s = (4A i,u /π) 1/2 = 27.2 mm; equivalent diameter of the slots open area). The extent of region 2 was calculated by matching the virtual surface area of a scaled underdrain to the filter cross-sectional area A f . The perpendicular distance from this virtual surface to the real underdrain was 18 mm for the design here employed.
However, the information from velocity vectors and streamlines nearby the pod observed in Figure 11 suggested that the flow was initially directed towards the region occupied by the entire pod surface area and not specifically towards the region occupied by the slots only. Flow trajectories deviated to reach the slots open area at the very end of their paths. Thus, we assumed that in the non-uniform zone, the fluid did not flow through a region with cross-sectional areas linearly decreasing from D f to D s but linearly decreasing from D f to D u along a distance d nu,2 and linearly decreasing from D u to D s along a distance d nu,1 , with D u (=94.2 mm) the equivalent diameter of the lateral surface of the pod-type underdrain (see Figure 12). Therefore, the effect of the non-uniform flow zone in the packed bed was calculated by expressing the pressure drop per unit length dp/dL in terms of the volumetric flow rate Q, since it was a magnitude conserved along any cross-sectional area of Figure 12. Here, the length L corresponded to the distance traveled along the streamline starting at the slots and ending at the top of the sand region. Thus, from Equation (9): where D was the diameter of the available cross-section, being calculated as: with D 1 = D s , D 2 = D u and D 3 = D f , and L 1 = 0 m, L 2 = d nu,1 and L 3 = d nu,1 + d nu,2 (see Figure 12). Therefore, Equation (13) for i = 1 and i = 2 applied to regions A and B in Figure 12, respectively. By substituting Equation (13) into Equation (12) we found: Water 2021, 13, 2179 20 of 25 whose solution was: with p i the pressure at L i , v s (= 4Q/ πD 2 ) the magnitude of the superficial velocity at L (or, equivalently, at D) and r D the diameter ratio defined as: Note that in the limit r D = 1, Equation (15) reduced to Equation (9) valid for uniform flow, as it should.
Results of the analytical expression (15) applied to evaluate the pressure drop through the path of a streamline compared with those extracted from the simulation and with the prediction of the analytical model proposed in [12,15] are shown in Figure 13. Results referred to streamlines #4 and #8 in Figure 11 with the parameters corresponding to the Ergun model with Q = 0.2 L s −1 or Q = 0.8 L s −1 for SS1 and MS media and two different heights h s of the granular bed. The d nu,1 value was chosen as half the minimum lateral distance between two consecutive slots, being d nu,1 = 1.5 mm. The d nu,2 value was calculated in a similar way as the non-uniform length L nu in Pujol et al. [15] but scaling the pod in the radial direction only since the commercial underdrain here used did not have slots in its top surface. This procedure assumed that flow non-uniformity was essentially caused by the effect of the pod openings, reaching a distance of influence equal to d nu,2 = 30 mm in our case. However, for streamline #4 in Figure 11, we added an additional 10 mm to the previous value of d nu,2 to take into account that its path was not perpendicular to the slots (and, therefore, it did not travel a distance d nu,2 = 30 mm within the non-uniform region) but clearly inclined (d nu,2 > 30 mm).
The analytical results from Equation (15) remarkably reproduced the pressure data along the entire path of the streamlines for all packed bed heights, flow rates and porous media types. In contrast, the analytical model in [15] only explained the behavior found in the uniform flow zone. Thus, simulations were almost matched by the analytical model (Equation (15)) that constrained the fluid to flow through a porous domain consisted of a cylindrical region (uniform zone) plus two consecutive contracted cones till reaching the slots (non-uniform zones; see Figure 12). The analytical formulation saved all computational costs related to CFD modeling (time and economical). However, its main shortcoming was the requirement of properly estimating the critical distances d nu,1 and d nu,2 . Pressure drop values from other underdrain designs would be needed to evaluate these distances in order to confirm the validity of the criterion applied here, which was based on geometrical considerations only.
On the other hand, from Figures 5-7, we observed that filter pressure drop ∆p f using the Ergun equation for porous media correctly reproduced experimental data for the whole range of superficial velocities tested when using silica sand but not microspheres. In the latter case, simulations substantially overpredicted ∆p f at high flow rates. The comprehensive study carried out in [13] about the pressure drop-flow rate correlations for packed beds of spheres concluded that Ergun equation overestimated the pressure drop when Re ε > 500. For the microspheres used in the present study, this condition was attained when v s > 0.47 m s −1 . At Q = 0.8 L s −1 , that value of superficial velocity was not reached in the uniform zone (C in Figure 12), but it was in the non-uniform zone (A in Figure 12). Therefore, we simulated the microsphere case with a momentum sink term that followed the pressure drop correlation proposed by Erdim et al. [13]: This expression was found to better represent the head loss of flows through packed spheres at high particle Reynolds number than the Ergun equation [13].
in the uniform flow zone. Thus, simulations were almost matched by the analytical m (Equation 15) that constrained the fluid to flow through a porous domain consisted cylindrical region (uniform zone) plus two consecutive contracted cones till reaching slots (non-uniform zones; see Figure 12). The analytical formulation saved all comp tional costs related to CFD modeling (time and economical). However, its main shortc ing was the requirement of properly estimating the critical distances , and Pressure drop values from other underdrain designs would be needed to evaluate t distances in order to confirm the validity of the criterion applied here, which was b on geometrical considerations only. On the other hand, from Figures 5-7, we observed that filter pressure drop ∆ ing the Ergun equation for porous media correctly reproduced experimental data for whole range of superficial velocities tested when using silica sand but not microsphe In the latter case, simulations substantially overpredicted ∆ at high flow rates. comprehensive study carried out in [13] about the pressure drop-flow rate correlation packed beds of spheres concluded that Ergun equation overestimated the pressure d when > 500. For the microspheres used in the present study, this condition wa tained when > 0.47 m s −1 . At = 0.8 L s −1 , that value of superficial velocity was reached in the uniform zone (C in Figure 12), but it was in the non-uniform zone (A Figure 12). Therefore, we simulated the microsphere case with a momentum sink t that followed the pressure drop correlation proposed by Erdim et al. [13]: This expression was found to better represent the head loss of flows through pac spheres at high particle Reynolds number than the Ergun equation [13].
Simulations carried out at = 0.8 L s −1 and ℎ = 117 mm using Equation (17 stead of Equation (10) obtained a reduction in ∆ of 0.50 kPa, though closer to obse tions than Ergun model, still higher than the measured value (see Figure 7). Figure 8, and the streamlines analysis above, confirmed the strong contributio Simulations carried out at Q = 0.8 L s −1 and h s = 117 mm using Equation (17) instead of Equation (10) obtained a reduction in ∆p f of 0.50 kPa, though closer to observations than Ergun model, still higher than the measured value (see Figure 7). Figure 8, and the streamlines analysis above, confirmed the strong contribution of the non-uniform flow in the packed bed to the filter pressure drop. In this region, flow velocities were higher than superficial velocities employed in Figure 2. Thus, values of pressure drop per unit length in the non-uniform region would correspond to extrapolated values (at higher velocities) of the regression fits shown in Figure 2. A qualitative observation of the predicted curves in Figure 2 for the quadratic expressions (Ergun and Darcy-Forchheimer) indicated that they almost exactly followed measured data for silica sand media (SS1 and SS2) but had a trend with a higher slope than the experimental data (though still within the error bars) for microspheres (MS). Therefore, the deviation when extrapolating the data trend would be expected to increase, and so the head loss predicted in the non-uniform zone.
A careful analysis of microspheres data obtained in the laboratory revealed some discrepancies between pressure drop per slab of height equal to 100 mm within the sand region (case h s = 317 mm). As it was shown in Figure 1a, the laboratory filter had pressure sensors vertically separated 100 mm between them (i.e., sensor p1 at height above the base plate h 1 = 100 mm, sensor p2 at h 2 = 200 mm and sensor p3 at h 3 = 300 mm). Figure 2 took into account pressure drop data between sensors p3 and p1. Figure 14 show data of pressure drop between sensors p3 and p2 and between sensors p2 and p1. We expected similar pressure drop data for both slabs since, from Figure 10, the level at h 1 = 100 mm was lying in the uniform flow zone. However, Figure 14 indicated lower pressure drop values in the slab from 200 mm to 100 mm than in the slab from 300 mm to 200 mm at high flow rates. This could be a consequence of an incipient particle retention in the uppermost packed bed layer. Though the test was carried out with tap water and porous media were washed prior to use, it might have occurred that some small particles were remaining in the system. This phenomenon was not observed in SS1 and SS2 cases, in which pressure drop data series for both slabs (h 2 − h 1 ) and (h 3 − h 2 ) were almost indistinguishable. The above point might be the reason that Darcy-Forchheimer equations overestimated the filter pressure drop for microspheres. However, Ergun equation for microspheres did not rely on any fitting with experimental data since we assumed a coefficient of sphericity equal to one. From Figure 2, the Ergun prediction for the momentum sink tended to overestimate the observed slope at higher velocities, which magnified the contribution of the non-uniform flow zone to the total pressure drop. An assumption of a sphericity value less than one for microspheres would not necessarily improve the prediction from the Ergun Equation. However, the use of the Ergun equation for the porous media correctly predicted the total filter pressure drop for the range of superficial velocities ( 50m h −1 ) in which the measured pressure difference at the two 100 mm vertical slabs ( Figure 14) were equal (overlapped error bars).
pressure sensors vertically separated 100 mm between them (i.e., sensor above the base plate ℎ = 100 mm, sensor p2 at ℎ = 200 mm and sensor p3 at Figure 2 took into account pressure drop data between sensors p3 and p1. Fi data of pressure drop between sensors p3 and p2 and between sensors p2 expected similar pressure drop data for both slabs since, from Figure 10, th 100 mm was lying in the uniform flow zone. However, Figure 14 indicated lo drop values in the slab from 200 mm to 100 mm than in the slab from 300 m at high flow rates. This could be a consequence of an incipient particle re uppermost packed bed layer. Though the test was carried out with tap wate media were washed prior to use, it might have occurred that some small p remaining in the system. This phenomenon was not observed in SS1 and which pressure drop data series for both slabs (ℎ − ℎ ) and (ℎ − ℎ ) we distinguishable. The above point might be the reason that Darcy-Forchheim overestimated the filter pressure drop for microspheres. However, Ergun eq crospheres did not rely on any fitting with experimental data since we assu cient of sphericity equal to one. From Figure 2, the Ergun prediction for th sink tended to overestimate the observed slope at higher velocities, which m contribution of the non-uniform flow zone to the total pressure drop. An ass sphericity value less than one for microspheres would not necessarily imp diction from the Ergun Equation. However, the use of the Ergun equation f media correctly predicted the total filter pressure drop for the range of supe ties (≲ 50 m h −1 ) in which the measured pressure difference at the two 100 slabs ( Figure 14) were equal (overlapped error bars). Finally, Table 6 summarizes the main advantages and disadvantages applying the momentum sink equations to evaluate the pressure drop in p All comments are exclusively related to the findings of the present study. T ues of the superficial velocity beyond which either Kozeny-Carman or Da underpredicted the pressure drop is not included in Table 6 since it was porous media type and packed bed height (see Section 3.1). Finally, Table 6 summarizes the main advantages and disadvantages found when applying the momentum sink equations to evaluate the pressure drop in porous media. All comments are exclusively related to the findings of the present study. Threshold values of the superficial velocity beyond which either Kozeny-Carman or Darcy equations underpredicted the pressure drop is not included in Table 6 since it was a function of porous media type and packed bed height (see Section 3.1).

Conclusions
We carried out CFD simulations of a pressurized granular bed filter with a commercial underdrain unit with the finest discretization found in the literature (16 × 10 6 elements required to mesh half of the laboratory filter). A minimum number of 12 elements across the slots (0.45 mm width) of the pod-type underdrain was used. Five types of equations to evaluate the pressure drop in the porous media were investigated: linear in v s (Darcy and Kozeny-Carman expressions), quadratic in v s (Darcy-Forchheimer and Ergun expressions) and a power function model in v s . Results were compared with experimental data obtained at different superficial velocities for three porous media (microspheres and two types of silica sand) and two different heights of the granular bed.
Total pressure drop results when applying all the previous porous media equations were very similar and lied within the experimental uncertainty range at least up to 38.3 m h −1 for both silica sand media with h s = 117 mm and h s = 317 mm. This superficial velocity corresponded to Re p values ranging from 6.9. to 982 depending on the media used, being Re p figures well above the recommended values for applying the linear Equations (Darcy and Kozeny-Carman).
Filtration modes with higher superficial velocities increased the discrepancies between models. The Ergun equation, based on flow and granular media properties, better followed the observed trend for both silica sand media and packed bed heights (NSE coefficient > 0.978 and RMSE < 1.3 kPa with respect to measured data) than using momentum sink expressions with terms fitted from experimental data. However, microsphere datasets were not accurately described by adopting pressure drop models with quadratic terms at high superficial velocities, likely due to some data inconsistencies for that range.
A detailed analysis of simulated data revealed important head losses inside the underdrain, mainly due to the sharp entrance into its exit pipe and the high velocity reached. For packed beds with low height h s = 117 mm, these minor and major hydraulics losses can represent more than 40% of the total filter pressure drop, so the inner design of the underdrain element is of paramount importance to reduce filter energy requirements.
Simulations also reported the influence of the non-uniform flow zone within the porous media. In quadratic models (Darcy-Forchheimer and Ergun), this region accounted for more than 40% of the total filter pressure drop at low granular bed heights. This confirmed the relevance to correctly design the water drainage so as the maximize the flow uniformity within the porous media.
An analytical model that took into account this non-uniform region within the packed bed was developed. The model divided the granular media intro three regions. The upper part with uniform flow, a second non-uniform zone where the flow was directed to the underdrain and a third non-uniform zone at very short distance from the slots where the flow finally moved towards the underdrain openings. Pressure drop values along streamlines obtained from the analytical model almost exactly matched CFD results, substantially improving previously proposed equations. The extent of non-uniform regions was derived from geometrical considerations of the underdrain element. However, further analyses with different underdrain designs would be needed to generalize the validity of the tree-zone analytical model.