CFD Modeling of Hydrocyclones—A Study of Efficiency of Hydrodynamic Reservoirs

The dynamics of hydrocyclones is complex, because it is a multiphase flow problem that involves interaction between a discrete phase and multiple continuum phases. The performance of hydrocyclones is evaluated by using Computational Fluid Dynamics (CFD), and it is characterized by the pressure drop, split water ratio, and particle collection efficiency. In this paper, a computational model to improve and evaluate hydrocyclone performance is proposed. Four known computational turbulence models (renormalization group (RNG) k-ε, Reynolds stress model (RSM), and large-eddy simulation (LES)) are implemented, and the accuracy of each for predicting the hydrocyclone behavior is assessed. Four hydrocyclone configurations were analyzed using the RSM model. By analyzing the streamlines resulting from those simulations, it was found that the formation of some vortices and saddle points affect the separation efficiency. Furthermore, the effects of inlet width, cone length, and vortex finder diameter were found to be significant. The cut-size diameter was decreased by 33% compared to the Hsieh experimental hydrocyclone. An increase in the pressure drop leads to high values of cut-size and classification sharpness. If the pressure drop increases to twice its original value, the cut-size and the sharpness of classification are reduced to less than 63% and 55% of their initial values, respectively.


Introduction
Hydrocyclones are useful devices for removing particles from liquids via a centrifugal field. In general, larger particles are mostly trapped in the underflow, while the lighter ones go to the overflow with most of the working fluid exiting, as can be observed in Figure 1. Some advantages of these devices are the low operation and maintenance cost and no moving parts. They also have widespread applications in the mining and oil industry. In the past decade, the main challenge for these devices has been to minimize the pressure drop and maximize the separation efficiency to improve hydrocyclone performance.
The d 50 (cut size) parameter is a key factor in hydrocyclone performance. The cut size is defined as the size at which a particle has a 50% probability of leaving the hydrocyclone through either the underflow or overflow. This factor is correlated with the collection efficiency, which indicates the number of particles that move toward the underflow. Finally, another critical factor is the pressure drop, which is an indicator of how much energy the hydrocyclone wastes, as a measure of the main system's loss. The main assessment tools used to analyze hydrocyclone performance are laboratory measurements and numerical methods such as Computational Fluid Dynamics (CFD). Compared to experimental procedures, numerical methods cost less and save time. Moreover, CFD has proved to be a reliable tool when used with a suitable choice of turbulence models. The most accurate and frequently implemented models are renormalization group (RNG) k- [1], the Reynolds stress model (RSM) [2], and large-eddy simulation (LES) [1,3].
Hsieh [4] built a hydrocyclone prototype and proposed a computational model based on the finite difference method. This study is considered a benchmark for analyzing the role of pattern flows in particle classification, which makes it valuable for evaluating new work in this area. Later work such as that presented in [3,[5][6]] used Hsieh's data for the validation of their numerical models. These studies concluded that only two turbulence models-RSM and LES-are able to model the anisotropy of the turbulent flow in hydrocyclones and thus accurately capture its flow dynamics. These two models, combined with the volume of fluid (VOF) multiphase model, can replicate the velocity profiles, air core diameter, and split water flow typically seen in a hydrocyclone. VOF tracks free surface or fluid interfaces. This technique is based on an Eulerian approach, and it is suitable not only for capturing the interface between the air core and the water in the hydrocyclone, but it is also less diffusive than other multiphase models such as the mixture model.
Despite multiple studies on hydrocyclone modeling, there is no unified methodology to evaluate the performance of these devices. In an attempt to create a standardized computational model and experimental method to characterize performance, numerous researchers have continued validating and testing new approaches to model hydrocyclone performance [1,[7][8][9], However, the majority of the work has consisted of experimental procedures. Commercially available hydrocyclones are based on traditional models, such as the Rietema and the Bradley design [10][11]. These commercial designs propose relations between the inlet diameter, vortex finder diameter, and cone angle. Such relations have remained invariant over the last two decades. This has encouraged academics to explore whether these parameters are adequate to operate a hydrocyclone at its best efficiency point.
Tang et al. [12] showed that a decrease in the size of the inlet nozzle leads to an increase in separation efficiency and a decrease in d50, but with this change, there is an evident effect on the increase in pressure drop. A better separation efficiency results in higher pressure drop values. Sun et al. [13] performed a multi-objective optimization; the variables they optimized were the pressure drop and the collection efficiency. In that study, it was found that the inlet and vortex finder diameter The main assessment tools used to analyze hydrocyclone performance are laboratory measurements and numerical methods such as Computational Fluid Dynamics (CFD). Compared to experimental procedures, numerical methods cost less and save time. Moreover, CFD has proved to be a reliable tool when used with a suitable choice of turbulence models. The most accurate and frequently implemented models are renormalization group (RNG) k-ε [1], the Reynolds stress model (RSM) [2], and large-eddy simulation (LES) [1,3].
Hsieh [4] built a hydrocyclone prototype and proposed a computational model based on the finite difference method. This study is considered a benchmark for analyzing the role of pattern flows in particle classification, which makes it valuable for evaluating new work in this area. Later work such as that presented in [3,5,6] used Hsieh's data for the validation of their numerical models. These studies concluded that only two turbulence models-RSM and LES-are able to model the anisotropy of the turbulent flow in hydrocyclones and thus accurately capture its flow dynamics. These two models, combined with the volume of fluid (VOF) multiphase model, can replicate the velocity profiles, air core diameter, and split water flow typically seen in a hydrocyclone. VOF tracks free surface or fluid interfaces. This technique is based on an Eulerian approach, and it is suitable not only for capturing the interface between the air core and the water in the hydrocyclone, but it is also less diffusive than other multiphase models such as the mixture model.
Despite multiple studies on hydrocyclone modeling, there is no unified methodology to evaluate the performance of these devices. In an attempt to create a standardized computational model and experimental method to characterize performance, numerous researchers have continued validating and testing new approaches to model hydrocyclone performance [1,[7][8][9], However, the majority of the work has consisted of experimental procedures. Commercially available hydrocyclones are based on traditional models, such as the Rietema and the Bradley design [10,11]. These commercial designs propose relations between the inlet diameter, vortex finder diameter, and cone angle. Such relations have remained invariant over the last two decades. This has encouraged academics to explore whether these parameters are adequate to operate a hydrocyclone at its best efficiency point.
Tang et al. [12] showed that a decrease in the size of the inlet nozzle leads to an increase in separation efficiency and a decrease in d 50 , but with this change, there is an evident effect on the increase in pressure drop. A better separation efficiency results in higher pressure drop values. Sun et al. [13] Fluids 2020, 5, 118 3 of 19 performed a multi-objective optimization; the variables they optimized were the pressure drop and the collection efficiency. In that study, it was found that the inlet and vortex finder diameter were significant factors in hydrocyclone performance. This study showed that increasing the size of the vortex finder leads to a decrease in separation efficiency as well as an increase in pressure drop.
The work of Vakamalla et al. [1] proposed working with parabolic profiles instead of a conical section. In order to reduce the computational time in CFD simulations, the air core was replaced by a physical rod. This simplified the simulation by converging faster and therefore allowing more simulations to be performed. However, that approach led to inaccurate estimations of the velocity profiles, pressure drop, and separation efficiency.
From the literature reviewed, it was found that the amount of hydrocyclone parameters result in many configurations that need to be evaluated. It is necessary to understand the role of individual geometric dimensions and their interactions on hydrocyclone performance and hence, finding the optimal design for certain specific operating conditions. Using the preceding principles, the main goals of this paper are to: (1) develop a computational model that evaluates and potentially helps to improve the performance of hydrocyclones and (2) analyze how structures (saddle point and vortices) could affect the performance of hydrocyclones as they are evaluated in different geometrical configurations.

Computational Fluid Dynamics
Since the flow in the hydrocyclone is turbulent, it is necessary to determine what turbulence model is suitable to replicate the hydrodynamic behavior within the device. RNG k-ε, RSM, and LES are the turbulent models evaluated in the present work. To check their accuracy, the obtained results from those models are compared against the experimental data by Hsieh [4].
Based on Reynolds decomposition, it is possible to split the instantaneous velocity into an averaged velocity u i and the fluctuating component of the velocity u i . Hence, the momentum equation for an incompressible fluid can be written as follows: The term −ρu i u j from Equation (3) is known as the Reynolds stress tensor, which leads to the turbulent closure problem. In consequence, in order to complete the equations, the tensor must be modeled.
2.1.1. RNG k-ε RNG, also known as the renormalization group, is a variant of the original k-ε model that considers the effects of rotation in the closing equations of turbulence. Using this model, the authors of [5] and [1] have found that the predictive advantages of this model for modeling hydrocyclones are not sufficient to calculate the pressure gradient, velocity profiles, and air core formation accurately. However, RNG k-ε was found by the authors to be useful to initialize the flow in order to reduce the computational time the solution takes to converge in LES simulations.
It is well known that turbulence is affected by rotation. The CFD modeling software Ansys© (Fluent, 2006) [14] proposes that when using RNG k-ε to model swirling flows, the effective turbulent viscosity is calculated based on the rotational tensor, a swirl factor (α), and the turbulent viscosity (µ t0 ). This modification takes the following form: The software documentation recommends values not greater than 0.07 for weakly to moderately swirling flows [14]. However, the hydrocyclone develops a strong swirling flow. For this case, the software suggests using a larger value, but it does not specify a value. As a consequence, the model in this study was calibrated using a trial and error process. After this process, it was found that the value of the constant α must be around 0.7; otherwise, there would not be a formation of the air core.
RNG k-ε is considered a two-equation model, where turbulent kinetic energy and the dissipation rate of energy are modeled as follows: The terms on the left of Equation (5) describe the accumulation and convection of dissipation energy. On the right side, the first three terms represent the production, dissipation, and diffusion of dissipation energy. The latter term is known as a source term and is described by the following expression: In Equation (7), the terms c µ and β are constant. This source term is unique for RNG k-ε. As Stendal [15] pointed out, the source term results in less disruption of ε (energy dissipation). Thus, less energy will be dissipated into heat. This implies that the RNG k-ε model is less dissipative than the other two-equation models. As a result of this phenomenon, this turbulence model is suitable to describe swirling pattern flows.

RSM-Transport Equation for Reynolds Stress Model
The accuracy of prediction of the RSM for hydrocyclones was demonstrated by [16]. The RSM model considers the anisotropic characteristic of the flow due to the fluid movement inside the hydrocyclone being highly turbulent. This turbulence model does not include the Boussinesq assumption. On the contrary, in this model, the term (ρu i u j ) is expressed by the transport equation below: The term on the right-hand side is: Turbulence Diffusion Stress Production Fluids 2020, 5, 118

of 19
Pressure Strain Dissipation term Rotation Production Fi j = 2ρΩ k u j u m ε ikm + u i u m ε jkm (13) To determine the pressure strain term, the Quadratic Pressure-Strain model was enabled. In addition, the boundary condition at the walls was treated by the standard wall function [17].

Large-Eddy Simulation (LES)
In this model, the Navier-Stokes equations are computed directly without any approximation on the large turbulence scales; the smaller scales are modeled as described in Equations (14)- (17). LES is based on the principle that larger turbulence scales are those that carry most of the momentum. LES is expected to give more accurate results compared to any models of the RANS family since it involves less modeling.
The disadvantage of LES is its high computational cost compared to RANS models, as LES requires a much finer mesh. In theory, the size of the mesh elements ought to be equal or lower than the Kolmogorov [18] length scale. This is required, since the model filters the Navier-Stokes equations from a specific length scale, directly simulating the larger scales and modeling the smaller vortices.
The Navier-Stokes filtered equations are represented below: Momentum: LES requires a SGS (sub-grid-scale) model to calculate the momentum transfer of scales smaller than the mesh, as can be seen in Equation (14). In this work, the Smagorinsky Lilly SGS model was implemented, where the τ sgs ij is modeled as the product between the turbulent viscosity and the strain rate tensor, which is represented as: The turbulent viscosity is determined algebraically by the following expression: L s is usually calculated as the cubic root of the finite volume at each node of the mesh in the high turbulence regions. However, the Ansys© Fluent software recommends calculating it based on the minimum value between the distance from the nearest wall and by the local scale of the mesh, ∆.
C s is the Smagorinsky calibration constant, with a value of 0.1 used in the present paper based on the work presented in [1,6]. However, in the study of DeSouza and Neto [19], the authors suggested that the value of this constant for hydrocyclones should be 0.15, but the tangential velocity profile is slightly overestimated. On the other hand, authors such as Saidi and Farhanieh [20,21] have managed to replicate flow characteristics in hydrocyclones with a dynamic version of the constant. In conclusion, there is no consensus in the literature for the value of the Smagorinksy constant. This implies that the Smagorinsky constant must be calibrated in order to resolve the smallest turbulence length scale accurately.
It is worth clarifying that it is beyond the scope of this study to deepen the understanding of the physical assumptions used for modeling the SGS tensor. The objective is to validate and compare the results for the different turbulence models as commonly reported in published works.

Multiphase Model
The development of the air core is a phenomenon that is part of the hydrocyclone dewatering process. For this work, the VOF (volume of fluid) multiphase model was implemented. In previous works [3,22] this model has been shown to correctly predict the shape and development of the air core. VOF is a technique based on a full Eulerian approach, where the purpose is to trace the interface between two continuous phases. The interface of air and water of the hydrocyclone studied is modeled using transport Equation (18).
where α p represents the volume fraction of the nth phase. The average density is described as follows: The RANS-VOF simulation is achieved by Favre Decomposition, as shown in the work by Fan et al. [23], where a formulation for Equation (18) is obtained in terms of an average volume fraction. This formulation is not limited to a RANS turbulence model but can also be applied to LES.

Particle Tracking
When analyzing a hydrocyclone, several governing equations are necessary. The first is the governing equation of particles, and the other is related to the particle-fluid coupling. According to [24], the most prevalent forces in hydrocyclones are drag force, pressure gradient, and the Saffman lift. Hence, the motion equation of particles can be expressed by Equation (20) F D u − u p is the drag force per unit mass. This force balances with the centrifugal force. The smooth particles are simulated as small spheres, so the F D term can be written as: where the Reynolds number is: Then, the drag coefficient is calculated as follows: The coefficients a 1 , a 2 , and a 3 are described by the Morsi-Alexander drag law. Finally, the Saffman Lift [25] Equation (20) is represented by the formulation: where C LS is the Saffman Lift Coefficient [25].
On the other hand, due to particles moving in a fluctuating velocity field, the turbulence effect is considered by enabling the stochastic tracking model. In addition, the DRW (Discrete Random Walk) was implemented in order to calculate the instantaneous velocity field.
The instantaneous velocity is a function of a constant that obeys a Gaussian distribution and the root-mean-square velocity component, as seen in Equation (25).
The root-mean-square velocity component can be calculated as follows (assuming isotropic conditions): In the case of the RSM model, the above expression is modified, so each component depends on the RMS value of the principal diagonal term of the Reynolds Stress tensor.
Finally, the interaction of small vortices and particles is considered by defining the lifetime scale of the eddies and the time integration trajectory of particles. The equations of those parameters are described below: For LES, the lifetime is equivalent to the LES time scale.

Numerical Simulation
As mentioned before, hydrocyclones are complex devices that involve many different physical principles, which increase the level of complexity of a stable CFD simulation of them. Previous work [26] suggests a step-by-step methodology that guarantees not only numerical convergence but also enables a further computational model validation of simple to complex cases.
In this work, the steps followed to carry out the simulation are shown in Figure 2. The first step only considers the water phase, and the simulation is carried out in steady state. Once a negative gauge pressure core is developed, the air phase is enabled by setting backflow at both outlets of the hydrocyclone. In the second step, the simulation is allowed to run until it reaches stability. Once the simulation is stable, pressure drop, water separation, and velocity profiles were validated by comparing the results with data from Hsieh's work. Finally, particles of different sizes are injected over the inlet surface and traced to the underflow. Due to the high computational cost of LES, calculations for this model are initialized using the stable solution obtained using RSM. only considers the water phase, and the simulation is carried out in steady state. Once a negative gauge pressure core is developed, the air phase is enabled by setting backflow at both outlets of the hydrocyclone. In the second step, the simulation is allowed to run until it reaches stability. Once the simulation is stable, pressure drop, water separation, and velocity profiles were validated by comparing the results with data from Hsieh's work. Finally, particles of different sizes are injected over the inlet surface and traced to the underflow. Due to the high computational cost of LES, calculations for this model are initialized using the stable solution obtained using RSM.

Mesh and Geometry
A 75-mm Hsieh hydrocyclone was simulated; the geometrical details are provided in Table 1.

Mesh and Geometry
A 75-mm Hsieh hydrocyclone was simulated; the geometrical details are provided in Table 1.  Figure 3 shows the structured mesh used in the simulations. A mesh convergence procedure was followed. Five meshes were evaluated, each mesh increasing 1.2 times the number of elements from the previous one. The tangential velocity profile was chosen as the convergence parameter, as seen in Figure 4. A mesh of 505,100 elements was ultimately selected due to the balance between good prediction accuracy and relatively low computational cost. Table 2 summarizes the quality criteria for the selected mesh.   Figure 3 shows the structured mesh used in the simulations. A mesh convergence procedure was followed. Five meshes were evaluated, each mesh increasing 1.2 times the number of elements from the previous one. The tangential velocity profile was chosen as the convergence parameter, as seen in Figure 4. A mesh of 505,100 elements was ultimately selected due to the balance between good prediction accuracy and relatively low computational cost. Table 2 summarizes the quality criteria for the selected mesh.

Boundary Conditions and Setup
The inlet is set as "inlet velocity" at 2.28 m/s corresponding to a mass flow of 1.117 kg/s. The walls were treated as non-slip conditions, which means that the velocity value is zero on the walls. Both outlets were defined with the boundary condition of a "pressure outlet" at atmospheric pressure. The air phase is set to enter the domain through both outlets when a negative pressure is developed inside the hydrocyclone. Finally, the turbulence conditions imposed at the hydrocyclone inlet for turbulent intensity and hydraulic diameter are determined by Equations (32) and (33), respectively.
As for the outlets, the default values of 5% and 10% were used for the turbulent intensity and the viscosity ratio, respectively. For the LES turbulence model, the inlet perturbation was generated by using a spectral synthesizer in Ansys Fluent. The turbulence intensity was 5% as in the RANS models. For the particle tracking, the size distribution used was the same at Hsieh's [4]. Limestone spherical particles of various diameters (35 µm, 25 µm, 17 µm, 12 µm, 9 µm, 6 µm, 4 µm) were found using the Rosin-Ralmmer distribution [27]; the distribution parameters are the mean diameter size and shape factor, with values of 17.96 µm and 1.44, respectively. Those parameters were found by a regression model. Additionally, the inlet velocity of each particle was set to 2.28 m/s. The simulations were run using Ansys© Fluent v19; the SIMPLE algorithm was used for the coupling of pressure and velocity. For the spatial discretization schemes used, the reader is referred to Table 3. This table also lists the temporary discretization and VOF multiphase model schemes. A 0.9 restitution coefficient was implemented at the hydrocyclone walls. All simulations are initialized using 4000 iterations of only the water phase. Then, the simulation is changed to a transitory state, and a time step of 1 × 10 −4 s is implemented for LES and RNG k-ε. For RSM, the time step was set at 2.5 × 10 −4 s. For the three models evaluated, 20 iterations per time step were used throughout the simulations.
It was observed that once steady state was reached, the residuals are kept below 1 × 10 −5 for the RMS turbulence model, which indicates a quantitative convergence. For each time step, the residual of continuity starts with a maximum value of 1 × 10 −3 and drops slightly below 1 × 10 −5 by the end of the time step for LES and RNG k-ε. The other residuals are below 1 × 10 −5 . The residual of continuity is expected to be higher than the other residuals when there are significant changes in the flow direction, as is the case for this study. However, it is possible to rely on these results because of the high quality of the mesh, as reported in Table 2.
Finally, for each turbulence model, the computationally measured variables are averaged after sampling the flow field variables during 0.5 s of simulation. For the RNG k-ε, RSM, and LES turbulence models, the simulation took approximately 36 h, 50 h, and 96 h to converge, respectively. All the simulations were carried out on an 8 core Intel Xeon E5-2630 v3 workstation with 96 GB of RAM.

Modeling of the Continous Phase
In order to assess the performance of the computational model developed in this study, the experimental data published by Hsieh was used to validate the results. There are two different axial positions at which experimental data is available, 60 mm and 120 mm from the top of the cylinder of the hydrocyclone. The tangential and axial velocity profiles are analyzed at both locations. Figure 5a,c show the tangential velocity profiles at the 60 mm and 120 mm locations, respectively. While all models accurately predict the tangential velocity profile in the central forced vortex, both RNG k-ε and LES have difficulties in the free vortex region closer to the walls where the velocity is overestimated by these models. At both locations, the LES model is better at predicting the peak velocity, although its radial position is slightly offset. RSM underpredicts the peak velocity, while RNG k-ε overpredicts it at both locations.
The axial velocity profiles, shown in Figure 5b,d, show relatively good agreement between all models and the experimental data, although both RNG k-ε and LES present some disagreement in the free vortex region closer to the wall.
The RSM model is part of the RANS family of turbulence models. In such models, the turbulent viscosity tends to diffuse the most pronounced velocity peaks, which might help explain the observed results. Although the RNG k-ε model was not expected to predict flow characteristics so accurately, this model combined with the higher than normal swirl factor constant (α s ) yields good results.
LES accurately predicts values of the tangential velocity (within 2.86% compared to Hsieh) in the forced vortex region, including the peak velocity. Nevertheless, it was observed that LES overestimates the values of tangential velocity in the free vortex region closer to the wall. As for the axial velocity profiles, LES also has some problems closer to the wall. These results suggest that a finer mesh might be required in order to accurately capture the flow hydrodynamics near the wall when using a fixed Smagorinsky constant model.
Nevertheless, the mesh used in this work is finer than that used in previous works [3,6], where better results were obtained for the velocity profile at the free vortex region. A possible explanation is that for values of y + ranging from 40 to 100 (as in previous works), the inertial sublayer is well described by the logarithmic law of the wall [14], as usually used in Ansys Fluent for this range of values of y + . The y + value of 20 falls into the buffer layer, for which such law of the wall might not be adequate. A finer grid might have been needed to improve the model.
On the other hand, it is worth pointing out that studies such as Rudolf [22] are consistent with the results of this work in which LES overestimates the velocity in the free vortex region. This indicates that trying to describe the hydrodynamics behavior of hydrocyclones using the LES model is not recommended unless using a finer mesh, which is expensive in terms of computational cost. The pressure contours predicted by the different models are shown in Figure 6. Table 4 presents the net pressure drop calculated for each model, showing differences of 2.3%, 16.0%, and 1.9%, respectively, for RNG k-ε, RMS, and LES. By comparing with results from [4], it can be seen that high values of pressure drop are correlated with the gradient of the tangential velocity in the radial direction. This relation can be explained by the pressure gradients equalizing the centrifugal forces in the flow. The higher the pressure gradients, the higher the velocity gradient.
Air core diameter and the water split ratio are the last two factors to be validated for each turbulence model. Figure 7 shows that the three models can predict the formation of the air core. The RNG k-ε turbulence model accurately predicts this phenomenon, although no study mentions it previously. Table 4 quantitatively summarizes these results. The air core diameter is better predicted by RSM and LES, while the water split ratio is better predicted by RSM. For the purpose of iteratively The pressure contours predicted by the different models are shown in Figure 6. Table 4 presents the net pressure drop calculated for each model, showing differences of 2.3%, 16.0%, and 1.9%, respectively, for RNG k-ε, RMS, and LES. By comparing with results from [4], it can be seen that high values of pressure drop are correlated with the gradient of the tangential velocity in the radial direction. This relation can be explained by the pressure gradients equalizing the centrifugal forces in the flow. The higher the pressure gradients, the higher the velocity gradient.
Air core diameter and the water split ratio are the last two factors to be validated for each turbulence model. Figure 7 shows that the three models can predict the formation of the air core. The RNG k-ε turbulence model accurately predicts this phenomenon, although no study mentions it previously. Table 4 quantitatively summarizes these results. The air core diameter is better predicted by RSM and LES, while the water split ratio is better predicted by RSM. For the purpose of iteratively improving the design of the hydrocyclones, it was found that RSM is the best model, as it accurately predicts the main variables, while maintaining a relatively low computational cost. Therefore, this model will be implemented for particle tracking and the exploration of new geometric hydrocyclone designs.      Table 4 summarizes the results of the validation of turbulence models for the continuous phase.

Modeling of a Discrete Phase
The discrete phase validation was executed in a steady state from the Ansys© Fluent ® post-processor. One-way coupling with Lagrangian particle tracking was implemented given the low concentration of particles (4.88 wt%). Figure 8 shows the comparison of experimental results with respect to the simulated CFD analysis.   Table 4 summarizes the results of the validation of turbulence models for the continuous phase.

Modeling of a Discrete Phase
The discrete phase validation was executed in a steady state from the Ansys © Fluent ® postprocessor. One-way coupling with Lagrangian particle tracking was implemented given the low concentration of particles (4.88 wt %). Figure 8 shows the comparison of experimental results with respect to the simulated CFD analysis. The "Fish Hook" effect is evident in the portion of the curve from 1 to 5 μm. The used turbulence model tries to capture this phenomenon but fails to estimate the probability of separation in that region. However, the parameters obtained-namely, the sharpness of separation and cut diameter (d50)-are experimentally consistent with the CFD simulation, and they have an accuracy of 4% and 7%, respectively. The previous results confirm that the one-way model, with steady tracking, is suitable for predicting the efficiency curve.
In summary, the computational model offers reliable results enabling the understanding of flow characteristics, which are not easily identified in experimental procedures. In terms of cost, this model converges with a few iterations, so it is possible to test more than one geometrical configuration to meet the needs and performance required for industrial applications.

Performance Variables
In this work, three performance criteria are considered: pressure drop, the sharpness of classification, and the efficiency separation. The parameter performance is described as follows: The "Fish Hook" effect is evident in the portion of the curve from 1 to 5 µm. The used turbulence model tries to capture this phenomenon but fails to estimate the probability of separation in that region. However, the parameters obtained-namely, the sharpness of separation and cut diameter (d 50 )-are experimentally consistent with the CFD simulation, and they have an accuracy of 4% and 7%, respectively. The previous results confirm that the one-way model, with steady tracking, is suitable for predicting the efficiency curve.
In summary, the computational model offers reliable results enabling the understanding of flow characteristics, which are not easily identified in experimental procedures. In terms of cost, this model converges with a few iterations, so it is possible to test more than one geometrical configuration to meet the needs and performance required for industrial applications.

Performance Variables
In this work, three performance criteria are considered: pressure drop, the sharpness of classification, and the efficiency separation. The parameter performance is described as follows: 4.1.1. Pressure Drop: where P in is the pressure at the inlet and P out is the pressure at the outlet evaluated at the top of the vortex finder.

Classification of Sharpness
Also known as imperfection, this parameter depends on three magnitudes, d 25 , d 50 , and d 75 , which are the probabilities of separation at 25%, 50%, and 75%, respectively.

Cut-Size Diameter
The d 50 is the parameter estimating the same probability of particles leaving the reservoirs through the overflow or the underflow.

Geometric Considerations
A computational experiment was conducted to understand how changes in the geometric parameters influence the performance of a hydrocyclone. Table 5 shows the three factors (D i , D o , L v ) that are more likely to affect the separation efficiency and pressure drop. Those factors are non-dimensional with respect to the hydrocyclone diameter. The level of each factor was selected based on previous studies [13] and [28]. However, in these studies, there were several restrictions. In their respective works, they approached the optimal location. In Figure 9, the design space is represented by a cube; each point is the possible combination of the levels of each factor. All values are investigated in five levels: −α, −1, 0, 1, and +α. The parameter α is calculated with 2k/4. In the case of three factors, alpha is equal to 1.68, which indicates that each point of the equidistant design space is 1.68 units from the central point.
where is the pressure at the inlet and is the pressure at the outlet evaluated at the top of the vortex finder.

Classification of Sharpness
Also known as imperfection, this parameter depends on three magnitudes, d25, d50, and d75, which are the probabilities of separation at 25%, 50%, and 75%, respectively.

Cut-Size Diameter
The d50 is the parameter estimating the same probability of particles leaving the reservoirs through the overflow or the underflow.

Geometric Considerations
A computational experiment was conducted to understand how changes in the geometric parameters influence the performance of a hydrocyclone. Table 5 shows the three factors (Di, Do, Lv) that are more likely to affect the separation efficiency and pressure drop. Those factors are nondimensional with respect to the hydrocyclone diameter. The level of each factor was selected based on previous studies [13] and [28]. However, in these studies, there were several restrictions. In their respective works, they approached the optimal location.
In Figure 9, the design space is represented by a cube; each point is the possible combination of the levels of each factor. All values are investigated in five levels: −α, −1, 0, 1, and +α. The parameter α is calculated with 2k/4. In the case of three factors, alpha is equal to 1.68, which indicates that each point of the equidistant design space is 1.68 units from the central point.

Evaluating Hydrocyclone Performace
A total of four hydrocyclones including Hsieh's hydrocyclone were evaluated; the results of this process are summarized in Table 6. In some configurations, two factors were maintained constant in order to determine the influence of one single factor on the hydrocyclone performance. For instance, the inlet size varies from 0.15 to 0.25 in configuration HC3 and HC4, but the cone length and vortex finder diameter remain constant. The same methodology is applied in the other configurations. From Figure 10, the observed configurations, HC1, HC2, HC3, and HC4, are more efficient in terms of separation and cut diameter. Based on preliminary observation, it is possible to identify that the higher the inlet velocity, the greater the particle classification. As a result, the pressure drop is more significant. The previously mentioned results can be confirmed by comparing the HC3 and HC4 configurations, where the only geometric parameter that varies in these two configurations is the inlet diameter, see Table 6.

Evaluating Hydrocyclone Performace
A total of four hydrocyclones including Hsieh's hydrocyclone were evaluated; the results of this process are summarized in Table 6. In some configurations, two factors were maintained constant in order to determine the influence of one single factor on the hydrocyclone performance. For instance, the inlet size varies from 0.15 to 0.25 in configuration HC3 and HC4, but the cone length and vortex finder diameter remain constant. The same methodology is applied in the other configurations. From Figure 10, the observed configurations, HC1, HC2, HC3, and HC4, are more efficient in terms of separation and cut diameter. Based on preliminary observation, it is possible to identify that the higher the inlet velocity, the greater the particle classification. As a result, the pressure drop is more significant. The previously mentioned results can be confirmed by comparing the HC3 and HC4 configurations, where the only geometric parameter that varies in these two configurations is the inlet diameter, see Table 6.  Figure 11 shows the streamlines of the velocity field in an axial cross-sectional cut. It is possible to see several stagnation points, since the fluid is losing momentum in the axial direction. Therefore, saddle points and vortices start appearing in the domain. In hydrocyclones, the saddle point represents the place where the separation process takes place. The primary purpose of improving the classification efficiency is to displace saddle points downward, near the spigot, where the smallest change of momentum occurs in the axial direction, and consequently, the particles are more likely to follow a path to the underflow. The vortices present in the hydrocyclones are undesirable, since they increase the residence time of the particles.  Figure 11 shows the streamlines of the velocity field in an axial cross-sectional cut. It is possible to see several stagnation points, since the fluid is losing momentum in the axial direction. Therefore, saddle points and vortices start appearing in the domain. In hydrocyclones, the saddle point represents the place where the separation process takes place. The primary purpose of improving the classification efficiency is to displace saddle points downward, near the spigot, where the smallest change of momentum occurs in the axial direction, and consequently, the particles are more likely to follow a path to the underflow. The vortices present in the hydrocyclones are undesirable, since they increase the residence time of the particles.
In Figure 11, the configurations are analyzed in pairs to elucidate how the flow patterns change from one geometry to another. Notice from Figure 11 and Table 6 that both configurations are similar. Both have the same length dimension of the hydrocyclone. It was found that the pressure drop and the cut diameter decrease with an increase in slenderness.
Furthermore, the classification is constant when the hydrocyclone length changes. Moreover, by analyzing Figure 11, it was observed that the inlet size is a significant factor, which is inversely proportional to pressure drop and directly proportional to the sharpness of classification and cut diameter. The saddle point and vortices are less chaotic than in the other configurations. In Figure 11, the configurations are analyzed in pairs to elucidate how the flow patterns change from one geometry to another. Notice from Figure 11 and Table 6 that both configurations are similar. Both have the same length dimension of the hydrocyclone. It was found that the pressure drop and the cut diameter decrease with an increase in slenderness.
Furthermore, the classification is constant when the hydrocyclone length changes. Moreover, by analyzing Figure 11, it was observed that the inlet size is a significant factor, which is inversely proportional to pressure drop and directly proportional to the sharpness of classification and cut diameter. The saddle point and vortices are less chaotic than in the other configurations.
Ultimately, a possible interaction between the inlet and vortex finder may exist; an increment of those parameters reports a lower value in the sharpness of classification and the pressure drop. In general, every configuration has one common attribute: the air core formation does not block the exit of the particles through the underflow. The larger the diameter of the vortex finder, the larger the width of the air core and the fewer particles will be separated from the working fluid.
From the HC2 and HC3 configurations, it is evident that they perform better than Hsieh's hydrocyclone, which is not a commercial hydrocyclone. The HC2 and HC3 configurations are in the range of the geometrical parameters proposed by the Rietema and Bradley brands, see Table 7. It would be possible to find the optimum operation point by exploring new designs similar to the configurations HC2 and HC3. Ultimately, a possible interaction between the inlet and vortex finder may exist; an increment of those parameters reports a lower value in the sharpness of classification and the pressure drop. In general, every configuration has one common attribute: the air core formation does not block the exit of the particles through the underflow. The larger the diameter of the vortex finder, the larger the width of the air core and the fewer particles will be separated from the working fluid.
From the HC2 and HC3 configurations, it is evident that they perform better than Hsieh's hydrocyclone, which is not a commercial hydrocyclone. The HC2 and HC3 configurations are in the range of the geometrical parameters proposed by the Rietema and Bradley brands, see Table 7. It would be possible to find the optimum operation point by exploring new designs similar to the configurations HC2 and HC3.  Table 8 summarizes the performance criteria of a hydrocyclone with respect to its geometrical parameters. The "−"sign indicates a decrease in the efficiency with a reduction in the factor, while the "+" sign indicates an increase. Table 8. Relationship between geometric parameters and hydrocyclone performance.

Parameters
∆P I dc This relationship enables more efficient hydrocyclones to be designed where their sharpness classification and cut size diameter is greater at the expense of minimal pressure drop.

Conclusions
In the development of the computational model, three turbulence models were evaluated. It was found that RSM can reproduce most of the flow patterns of the continuous phase at a relatively low computational cost. RSM was more time consuming than the RNG k-ε, due to the seven fields resolved per iteration in the former, but it was able to more accurately predict the main variables of interest. On the other hand, LES will require a finer grid to reach convergence and it is overall a more computationally expensive model. The RNG method took 36 h to run compared to 96 h for LES. It was also found that RSM can predict the classification curve in hydrocyclones point by point, making it suitable for the characterization of hydrocyclone performance.
The last significant finding was that some vortices and saddle points move as the geometry changes. These stagnation points are responsible for the low efficiency volume inside hydrocyclones. In order to ensure better performance, it is necessary to work with slender hydrocyclones, lower inlet diameters, and a lower vortex finder diameter. With this, it is also guaranteed that the momentum interchange remains in the axial direction. Consequently, the particles can leave the domain through the spigot. Finally, with these geometric variations, it is possible to build more compact hydrocyclones that meet the space limitation in any hydraulic system. For future work, we recommend applying the response surface methodology to mathematically determine the dependence between the evaluated geometrical factors and the performance of these filters.