Performance of Reynolds Averaged Navier–Stokes and Large Eddy Simulation Models in Simulating Flows in a Crossflow Ultraviolet Reactor: An Experimental Evaluation

Abstract

Computational Fluid Dynamics (CFD) has been increasingly adopted as a design tool for the simulation of UV disinfection efficiency and the optimization of the configuration of a UV reactor. However, the performance of CFD with different turbulence closures may vary significantly. In the present study, an experimental evaluation was performed to assess the performance of CFD with five Reynolds Averaged Navier–Stokes (RANS) turbulence closures and three Large Eddy Simulation (LES) sub-grid scale (SGS) models. A simplified crossflow reactor with a single lamp sleeve was fabricated for the experimental measurements and numerical simulations. Overall, the superior performance of LES compared to RANS models in flow predictions within a complex configuration is demonstrated.

1. Introduction

Ultraviolet (UV) disinfection, as an alternative among the assortment of treatment methods, has shown high efficiency in water, wastewater, and ballast water treatments [1]. The microbial inactivation inside a UV reactor is governed by the total amount of UV light that the microorganism receives, namely the UV exposure dose. Thus, both the UV light intensity distribution and water flow prediction are needed for the calculation of the UV dose. In the past few decades, integral approaches to UV disinfection simulation have been commonly employed, which consist of a fluence rate model for UV light intensity distribution, a numerical model for flow prediction, and a microbial inactivation kinetic model for the fluence response of target microorganisms. These models have been extensively studied individually [2,3,4,5] and systematically [6,7,8,9]. Among these studies, Computational Fluid Dynamics (CFD) has been commonly adopted as a powerful numerical model for flow prediction to predict the performance of a UV reactor and to conduct the configurational optimization [6,7,8,9,10,11,12,13,14].

Despite the increasing adoption of CFD, it is worth noting that flow prediction is challenging due to the complex flow pattern inside compact UV reactors with many UV cylindrical lamps [7,15,16]. As far as we are aware, Liu et al. [7] first showed that the efficiency of a disinfection reactor is sensitive to the RANS turbulence model selection. Wols et al. [15] further investigated the performance of CFD simulations by comparing the standard k-epsilon turbulence model and Large Eddy Simulation (LES) with experimental measurements. LES was found in their study to predict disinfections better than RANS models. Zhang et al. [16] pointed out in their review work that unsteady flow structure effects and the complexity of reaction systems are key challenges for disinfection process simulations. Shah et al. [17] summarized the CFD performances in previous studies and again pointed out the reduced accuracy of RANS in flow predictions within UV reactors. Despite the reduced accuracy of RANS with certain turbulence models, studies using RANS models [11,12,18,19,20,21,22,23], however, remained the dominant approach for UV reactor simulations in recent decades due to their efficiency. We note particularly that correct simulations of flow separation around cylindrical UV lamps are key to the improvement of numerical accuracy. Recently, an inspiring study by Bose and Yeo [24] on atmospheric boundary layer flow problems again demonstrated the weakness of RANS models in predicting the flow separation behind a smooth hill, which is comparable to the flow behind a UV lamp. All in all, the usage of RANS models for UV disinfection simulations could result in the reduced accuracy of flow predictions as well as consequent inactivation efficiency. Therefore, it is necessary to evaluate the performance of the RANS and LES models for a satisfactory prediction inside complex UV reactors in order to develop an improved understanding of the tradeoff.

In the present study, we performed an experimental validation for five RANS turbulence closures (models) and three LES sub-grid scale (SGS) models, all of which are well-developed and commonly used models. A typical configuration of the UV disinfection for ballast water treatments is a closed crossflow reactor with UV lamp sleeves usually placed perpendicular to the flow with a staggered pattern, which involves the wall interaction. Therefore, a customized closed conduit in which the flow passed a horizontally orientated cylindrical sleeve was adopted for measurements of the velocity field with the approach of Particle Image Velocimetry (PIV). The experimental and numerical methodologies are described in Section 2 and Section 3, respectively. The results and discussion are presented in Section 4. Finally, conclusions are drawn in Section 5.

2. Experimental Methodology

2.1. Experimental Setup

A lab-scale model was fabricated (Zoncepz Solutions, Singapore, Singapore) using transparent acrylic sheets with a thickness of 10 mm to quantify the hydrodynamics behind the single lamp sleeve of crossflow UV reactors. As illustrated in Figure 1, the cylinder (sleeve) had an outer diameter (Ø_o) of 30 mm and was placed at the mid-depth of the model. The reactor was designed with a rectangular cross-section of 71 mm (height) × 50 mm (width), as it was more convenient for PIV measurements to be conducted with a flat surface rather than a curved surface (as with a circular cross-section) [25]. The reactor began with a circular inlet tube that had an inner diameter (Ø_i) of 24 mm and was followed by an abrupt change in the cross-sectional area of the conduit, which led to the flow separation. To achieve a uniform flow, a porous plastic block was placed near the inlet tube and an acrylic sheet with 56 holes was attached to the bounded wall to hold the porous block. The upstream region occupied approximately 1.3 m, within which the flow became fully turbulent.

The PIV technique was adopted for measurement, with the seeding particle size chosen to be 50 µm for all the experiments. To capture the motion of seeding particles in the model, a dual-cavity frequency-doubled pulsed Nd:YAG laser and a charge-coupled device (CCD) camera were combined as an integral system. The pulsed laser system emitted a light sheet with a wavelength of 532 nm for illumination. The CCD camera captured images with a resolution of 1600 × 1186 at a frame rate of 15 Hz. The specifications of the system were summarized in [26].

2.2. Experimental Procedures

During the experiments, the flow circulation was maintained by connecting the model with a reservoir tank. Water was pumped continuously from the reservoir, and the flow rate was controlled by a digital flow meter. The laser light sheet was adjusted to coincide with the center plane parallel to the vertical direction, and the camera was focused on the laser sheet to capture the scattered light. The seeding particles were added to the tank to obtain a fully mixed solution with a concentration of 0.1 g/L. The time interval between two frames was adjusted according to various flow rates and was determined through calibration runs. The measurements were conducted at three flow conditions, as shown in

Table 1. Flow conditions of experiments.

Cases	1	2	3
Flow Rate Q0 (L/s)	0.63	1.18	1.55
Cross-Sectional Averaged Velocity U0 (m/s)	0.176	0.332	0.437
Re	6600	12,400	16,400

.

To investigate the hydrodynamics around the cylinder, especially the variations in velocity profiles behind the cylinder where a recirculation zone (wake region) is usually present in turbulent flows [27], a total of 279 image pairs were recorded at 2 Hz during each measurement which lasted for 138 s. These image pairs were post-processed using the Dantec Dynamic Studio software for the time-average velocity field. The velocity field was first derived through “Adaptive Correlation”, based on recorded images, and then filtered by “Moving Average Validation”. A time-averaged velocity vector map was finally generated by averaging the vector fields at individual time instances. To examine the velocity field, velocity profiles were extracted along vertical lines at seven locations, which included 75 mm upstream, and then 50, 75, 100, and 125 mm downstream. The streamwise velocity profiles along the streamwise centerline (y = z = 0) were also extracted to investigate velocity variation along the centerline.

3. Numerical Methodology

3.1. Governing Equations

The RANS equations can be written in Cartesian tensor form [28,29] as a Continuity equation:

$$\frac{\partial \rho }{\partial t}+{\frac{\partial }{\partial x}}_i\left(\rho u_i\right)=0$$

(1)

and momentum equation:

$$\frac{\partial }{\partial t}\left(\rho u_i\right)+\frac{\partial }{\partial x_j}\left(\rho u_i{u}_j\right)=-\frac{\partial P}{\partial x_i}+\rho g_i+\frac{\partial }{\partial x_j}\left[\mu \left(\frac{\partial u_i}{\partial x_j}\right)\right]+\frac{\partial {\tau }_{ij}}{\partial x_j}$$

(2)

where $g_i$ is the gravitational acceleration in the i direction, $\rho$ is the fluid density, $P$ is the pressure, $t$ is the time, and $u$ is the velocity. Compared to the Navier–Stokes equations, the RANS equations involve additional terms, ${\tau }_{ij}=-\rho \overline{u_i^{\prime }{u}_j^{\prime }}$, which are the Reynolds stresses. These terms must be solved with additional equations to ensure the closure of equations.

For the k-ε-type and k-ω-type closures, the Boussinesq hypothesis, based on the assumption that the stresses can be incorporated through turbulent viscosity or eddy viscosity, is adopted in the following manner [30]:

$${\tau }_{ij}=-\rho \overline{u_i^{\prime }{u}_j^{\prime }}={\mu }_t\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\frac{2}{3}\left(\rho k+{\mu }_t\frac{\partial u_k}{\partial x_k}\right){\delta }_{ij}$$

(3)

where ${\mu }_t$ is the turbulent viscosity, $u_i$ is the mean velocity component, ${\delta }_{ij}$ is the kronecker delta, and k is the turbulent kinetic energy which is calculated by

$$k=\frac{\overline{u_i^{\prime }{u}_i^{\prime }}}{2}=\frac{1}{2}\left(\overline{u^{\prime 2}}+\overline{v^{\prime 2}}+\overline{w^{\prime 2}}\right)$$

(4)

Three k-ε closures, including the standard, Re-normalization Group (RNG), and realizable k-ε closures, together with two k-ω closures, including the standard and Shear-stress transport (SST) k-ω closures, are introduced to solve the turbulent viscosity. Details are described in Appendix A.1, Appendix A.2, Appendix A.3, Appendix B.1 and Appendix B.2.

In LES models, eddies are filtered into large and small parts based on the local grid sizes, and large eddies are computed directly by the instantaneous Navier–Stokes equations, while small ones are modeled based on the Boussinesq hypothesis. The filtered continuity, momentum, and concentration transport equations in Cartesian coordinates for LES models are as follows [31]:

Continuity equation:

$$\frac{\partial \rho }{\partial t}+{\frac{\partial }{\partial x}}_j\left(\rho {\tilde{u}}_j\right)=0$$

(5)

Momentum equation:

$$\frac{\partial }{\partial t}\left(\rho {\tilde{u}}_i\right)+{\frac{\partial }{\partial x}}_j\left(\rho {\tilde{u}}_i{\tilde{u}}_j\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\mu \left(\frac{\partial {\tilde{u}}_i}{\partial x_j}\right)\right]+\frac{\partial {\tau }_{ij}}{\partial x_j}$$

(6)

where the tilde indicates spatially filtered variables; ${\tau }_{ij}=-\rho \widetilde{u_i{u}_j}-\rho {\tilde{u}}_i{\tilde{u}}_j$ are the sub-grid scale (SGS) stresses.

The sub-grid-scale turbulence closures are typically based on the Boussinesq hypothesis, which calculates the sub-grid-scale turbulent stresses as follows [31,32]:

$${\tau }_{ij}-\frac{1}{3}{\tau }_{kk}{\delta }_{ij}=-2{\mu }_t{\tilde{S}}_{ij}$$

(7)

where ${\mu }_t$ is the SGS viscosity, ${\tau }_{kk}$ is the isotropic part of the SGS stresses which can be added to the filtered static pressure term or neglected, and ${\tilde{S}}_{ij}$ is the strain rate tensor given by

$${\tilde{S}}_{ij}=\frac{1}{2}\left(\frac{\partial {\tilde{u}}_i}{\partial x_j}+\frac{\partial {\tilde{u}}_j}{\partial x_i}\right)$$

(8)

In the present study, three SGS models, including Smagorinsky–Lilly, Dynamic Kinetic Energy (DKE), and Wall-Adapting Local Eddy-Viscosity (WALE), are introduced to solve the SGS viscosity, and details are described in Appendix C.1, Appendix C.2 and Appendix C.3.

3.2. Numerical Setup

Figure 2 shows the three-dimensional computational domain created in Gambit which was a mesh generator for ANSYS Fluent, where the origin of the coordinate was set to the center of the cylinder (tube sleeve). The computational domain was based on the physical setup in the experiments, as shown in Figure 1, with the same dimensions but a simplification of the upstream configuration. Within the physical setup shown in Figure 1, there was a substantial distance (over 1.5 m) upstream for the flow to develop into a fully turbulent state before reaching the cylinder. However, the region of interest was mainly the area behind the sleeve where flow separation happened. To save computational effort, only a 0.3 m length in the upstream domain was simulated, and it was numerically examined to be sufficient to cover the upstream region (about 0.05 m based on the experiments as presented in the following sections) affected by the presence of the cylinder.

3.3. Meshing and Boundary Conditions

The domain was meshed into structural meshes with 2.0 million cells. The grid size ranged from 0.5 mm to 2 mm, and the near-wall regions were further refined by adaptation. To check the grid independence, the results from three grid schemes were compared with each other. Figure 3 shows the velocity profiles of Re = 12,400 at x = 100 mm downstream simulated based on the standard k-ε closures and the LES approach with the WALE SGS model with three different grid schemes. The grid scheme with 2 million cells was found to be sufficient for both k-ε and LES simulations. The inlet boundary condition was set to a uniform velocity, with three corresponding inlet velocities of 0.176, 0.332, and 0.437 m/s. The turbulence intensity was estimated to be 2% for the inlet. The outlet was set as a pressure outlet and the other boundaries were set as no-slip wall boundaries.

3.4. Solver Settings

ANSYS Fluent was adopted as the solver for the simulations, and a pressure-based transient solver was used. In the present study, the LES approach with SGS models, including the Smagorinsky–Lilly (S-L), WALE, and DKE models, as well as the RANS approach with the standard, RNG and realizable k-ε, and standard and SST k-ω closures, was used for comparisons. The near-wall modeling approach was chosen for LES, and the log-law wall function was chosen for the k-ε closures and k-ω closures. The time step was adjusted using the Courant–Friedrichs–Lewy condition [33] so that the Courant number was less than one. The convergence criteria for the continuity, velocities, and k values were all set to 10⁻⁶.

To determine the necessary duration for time averaging of the LES simulation results, the velocity profiles at x = 100 mm obtained from various time durations are compared in Figure 4, where the results were obtained from the LES (results were time-averaged from t = 0 to the specific time) and standard k − ε closure predictions (results were extracted at the specific time). Even though the RANS model solved the flow in a Reynolds averaged manner, the sufficient simulation time still needed to be determined to avoid the effect of the initial condition. From the comparison of different time durations, the simulation time of 15 s was found to be fully sufficient for RANS predictions, and the duration of 25 s was sufficient for LES.

4. Results and Discussion

The mean velocity characteristics at the center plane were obtained from various turbulence closures and were compared with the present experimental data. The velocity profiles along the vertical direction (z axis) were first extracted from the center plane (y = 0) at various streamwise locations (x axis), including x = −75 mm, 50 mm, 75 mm, 100 mm, and 125 mm. The effect of the Reynolds numbers could be assessed by comparing the velocity profiles at the same streamwise location under various Reynolds numbers. The longitudinal velocity profiles along the centerline (y = z = 0) from x = −75 mm to x = 300 mm were also extracted to investigate the velocity variation along the centerline in this section.

4.1. Velocity Profiles in Upstream Region

It is critical to examine the steady flow profile in the upstream section before the cylinder. Figure 5 shows the vertical velocity profiles in the upstream at x = −75 mm with the three Reynolds numbers. From the experimental data, it can be observed that the velocity profile was slightly more uniform for higher Reynolds numbers. This suggests that the flow field might not be fully developed at the lowest Re, and it was probably due to the insufficient upstream length.

Despite the experimental imperfection, all the numerical results obtained from the different approaches collapsed into a single curve and were somewhat close to the experimental data, revealing the similar performance of these approaches in the prediction of steady upstream flow. A slight difference in the velocities at the centerline between the experimental data and the predicted results can be observed at the highest Re, with the variation around 2%, which was within the experimental uncertainty.

4.2. Velocity Profiles in Wake Region

It was concluded from Figure 5 that the different approaches all provided reasonably good predictions of the steady flow profile in the upstream section. However, behind the cylinder, there existed a wake region in which velocity fields were significantly affected by the cylinder, and the performance of conventional RANS closures was reported to yield poor accuracy [7] in this region. Therefore, the characteristics of the flow in the wake region behind the cylinder required further investigation. In this section, the mean velocity profiles were extracted from various turbulence closures as well as experimental results at different streamwise locations for comparison.

Figure 6 shows a comparison of the velocity profiles at x = 50 mm obtained from the various CFD approaches and experimental observations. From the figure, the velocity profile featured two peaks and one trough, with a local minimum along the centerline. The two peaks were due to flow acceleration with the presence of walls both at the top and bottom where the velocity decreased to zero. The presence of a trough located at the centerline in the wake region was consistent with previous studies [34,35] when the flow passed a cylinder. From the figure, both LES and RANS models with a k-ω closure were found to yield better agreement with the experiments, while RANS with a k-ε closure yielded much higher peaks and lower local minimum values. In terms of the trough, LES produced better predictions with a higher Re.

Further downstream, the velocity profiles at x = 75 mm were extracted and are plotted in Figure 7, in which the curves also generally had two peaks and one trough. Compared with Figure 6, the differences between the extreme values (peaks and troughs) became smaller, indicating the effect of the cylinder was weakened. From the figure, RANS with the three different k-ε closures still far over-predicted the peak values and under-predicted the minimum values. RANS with the three different k-ω closures gave much better predictions of the minimum values but over-predicted the peak values. In comparison, LES, especially the DKE and WALE SGS models, showed its capability of simulating both extreme values reasonably well.

Figure 8 shows the velocity profiles at x = 100 mm. Compared with Figure 6 and Figure 7, the flow profile tended to be much more uniform. In general, RANS models with k-ε closures were unable to reproduce this trend. Among the k-ω closures, the standard k-ω closure better predicted the minimum values at the centerline, but both the standard and SST k-ω closures over-predicted the peaks. In comparison, the LES predictions were consistent with the experimental results. Within LES, the WALE and DKE SGS models were found to be better than the S-L model, which tended to yield more uniform results at large and medium Re values.

In summary, RANS models with k-ε closures were not able to reproduce velocity variations in the wake region, and the k-ω closures provided better predictions but over-predicted the peaks when the flow went downstream from the cylinder. In general, LES, especially the WALE and DKE SGS models, provided the closest results to the experimental data in the wake region.

4.3. Velocity Profiles beyond Wake Region

With increasing distance from the cylinder to the downstream, the effect of the cylinder eventually disappeared, and the velocity profiles became uniform in general. Figure 9 shows the velocity profiles at x = 125 mm, and from the figure, the velocity became almost uniform, proving that the influence of the wake region can be negligible at this point.

From the figure, it can be observed that RANS models with k-ε closures were not able to simulate this transition and they tended to predict a much larger wake region. The k-ω closures generally captured this transition better, even though the standard k-ω closure tended to reproduce large wake regions at high Re values. However, both the standard and SST k-ω closures predicted much higher velocities compared with the experimental data. In comparison, LES was able to capture this transition well and the predicted velocities were closer to those in the experimental data. Among all the LES models, the LES S-L model was the best performer.

4.4. Longitudinal Velocity Profiles at Centerline

To better understand the flow characteristics, the magnitude of streamwise velocity along the centerline of the reactor was plotted from x = −75 to 300 mm, as shown in Figure 10. It was found that the centerline velocity started to drop from x = −50 mm until reaching zero at the surface of the cylinder. Once the flow passed the cylinder, the velocity increased from zero to a local peak, then decreased again to a local minimum located at around x = 42 mm and seemed invariant with Re. The minimum location was the location where the direction of streamwise velocity changed from negative to positive [35]. Beyond that, the velocity began to positively increase until reaching a steady state.

4.5. Discussion

The better performance of the LES model than any of the RANS models shown in the present study was anticipated and could be explained by the superior theory of LES to RANS. In turbulent flows, the transports of mass, momentum, energy, and other quantities are mainly affected by large eddies which are strongly determined by the boundary conditions, and thus are anisotropic [28]. On the contrary, small eddies are isotropic or tend to be isotropic and are rarely affected by boundaries. In the two types of RANS models, namely the k-ε-type and k-ω-type models, the Boussinesq hypothesis is adopted, which assumes that Reynolds stresses can be incorporated through turbulent viscosity or eddy viscosity (expressed by Equation (3)), to solve Navier–Stokes equations. The underlying assumption of the Boussinesq hypothesis is that the turbulent viscosity is assumed to be isotropic, which may not be true for large eddies, and would therefore lead to false results, which was revealed in the present study.

In contrast, eddies in LES are filtered into large and small parts based on the local grid sizes, and LES models resolve Navier–Stokes equations by computing the relatively large eddies using direct numerical simulation without any hypotheses and simulating the small eddies using appropriate assumptions, namely the SGS models. LES is a workable alternative because previous observation reveals that the scale of eddies ranges from the macro-scale to the micron-scale and most of the energy is contained in the large eddies [36]. Therefore, by direct simulations of large eddies, LES is expected to yield better predictions of the flow than RANS closures. For small eddies, the SGS models can give reasonable results even with relatively coarse grids.

In the present study, among the RANS closures, the k-ω closures yielded better predictions than the k-ε closures in the wake region. This phenomenon can be attributed to the better ability of k-ω closures, by substituting the dissipation rate in k-ε closures with the specific dissipation rate (in Appendix B.1 and Appendix B.2), in coping with the adverse pressure gradients, wall boundary layers, and shear flows [37,38] which exist in the wake region behind the cylindrical sleeve (UV lamp). However, such an ability of k-ω closures adversely becomes a weakness in the far field beyond the wake region, as shown in Figure 10. It is worth noting that SST yielded better predictions than the standard k-ω closure in the far field region, due to its combination of the standard k-ω in the near-wall region and the k-ε closure in the far field [39].

5. Conclusions

The hydrodynamics of the flow passing a cylinder within a closed conduit was investigated by both RANS and LES approaches with various turbulence closures and LES SGS models, respectively. Overall, RANS models with k-ε closures yielded the least accurate results with far under-estimations of the longitudinal velocity in the wake region, while the k-ω closures were in satisfactory agreement with the experimental data near to the cylinder but continuously over-predicted in the far field; LES provided the best fit with the experimental data. Overall, the results illustrated that LES is a better approach than RANS in predicting the flow field around complex boundaries, including possible flow separation, and is thus more suitable for flow field predictions within a UV reactor. In view of the need for reliable predictions of the flow field and inactivation, LES with suitable SGS models is therefore recommended for UV reactors with intensive wall interactions. It is recommended that the LES approach for both the mean flow and Reynolds stress distributions, integrated with a fluence rate model for UV light intensity and a kinetic model for microbial inactivation, should be further evaluated in the future for the design of UV disinfection reactors.