The impact of turbulence models on the simulation of residence time distribution in a multichamber ozone contactor

The ﬂ ow and tracer transport in an ozone contactor is simulated by using computational ﬂ uid dynamics. The standard k - 1 model, RNG k - 1 model, Realizable k - 1 model and SST k - v model are used to describe turbulence. A step change method is used to simulate the residence time distribution. The residence time and cumulative residence time are compared with labora-tory experiments. All turbulence models can capture the feature of the residence time distribution and cumulative residence time distribution. The residence time distribution in the initial period is sensitive to the used turbulence model, which can affect the value of MDI. The standard k - 1 model behaves better among four turbulence models. freestream SST v v k 1


INTRODUCTION
Ozone disinfection plays an important role in the advanced treatment of drinking water because of its strong disinfection effect, high disinfection efficiency, wide range of sterilization (Helen et al. 2003;Lee et al. 2015;Verma et al. 2016;Kolosov et al. 2018). However, bromate, which is a potential carcinogen (group 2B), may be produced when water contains bromide ions (Kim et al. 2007). Therefore, it is very important to control the formation of bromate in the ozone disinfection process. The amount of bromate produced in ozone disinfection process can be affected by ozone dosage, reaction time, ozone distribution, reaction temperature, and concentration of bromide ion. To control the bromate production and improve the disinfection efficiency, the hydrodynamics, raw water quality and operating conditions should be carefully determined during ozone disinfection.
Residence time distribution has been widely used to evaluate the hydrodynamics in ozone contactors (Wilson & Venayagamoorthy 2010). Commonly used values are t 10 , t 50 , t 90 , representing the ratio of the time for 10, 50 and 90% of inlet tracer flowing out of the ozone contactor, respectively (Ender & Mustafa 2018). In addition, the Morrill dispersion index (MDI), which is defined as the ratio of the characteristic residence time t 90 to t 10 , is used to evaluate the relationship between the mixing level of the fluid and the relative distribution of the tracer between t 10 and t 90 in the ozone contactor. Recently, a new parameter, namely recirculation factor, to characterize the contact tank, which is defined as the difference of unity and the ratio of the volume occupied by the recirculation zones and the overall volume of the reactor (Bruno et al. 2021). These parameters can be obtained by the tracer experiment (Martin et al. 1992;Wols et al. 2008;Kim et al. 2010a) and the computational fluid dynamics (CFD) (Choi et al. 2004;Zhang et al. 2007Zhang et al. , 2014aTalvy et al. 2011;Zou et al. 2017;Aparicio-Mauricio et al. 2020). The contact tank geometry has a significant on the mixing efficiency. Contact tank geometry has a significant effect on the mixing efficiency. Slot-baffled and perforated baffled contact tank have been developed to improve the mixing efficiency and optimize disinfection treatment (Aral & Demirel 2017;Demirel & Aral 2018a;Bruno et al. 2020Bruno et al. , 2021. Flow in ozone contactors are usually turbulent. The simulation of turbulence plays an important role in the prediction of residence time distribution. Large-eddy simulation (LES) was viewed as a powerful tool in the simulation of ozone contactors (Wang & Falconer 1998;Kim et al. 2010bKim et al. , 2013Zhang et al. 2013aZhang et al. , 2013bZhang et al. , 2014bDemirel & Aral 2016). Wang & Falconer (1998) investigated residence time distribution using k À 1 and largeeddy simulation (LES). The peak time was mainly influenced by the advection effect. Kim et al. (2010bKim et al. ( , 2013 conducted a three-dimensional numerical simulation of flow and tracer transport in a multichamber ozone contactor using LES. It was concluded that LES is a powerful tool in the simulation of ozone contactors. Zhang et al. (2013aZhang et al. ( , 2013bZhang et al. ( , 2014b simulated the flow and tracer transport in the ozone contactor using the open-source numerical library OpenFoam and both the k À 1 model and large-eddy simulation (LES). Good agreement between the numerical data and the experimental data were obtained for cumulative residence time distributions. Demirel & Aral (2016) simulated mixing efficiency in multi-chamber contact tanks based on vorticity field in combination with both k À 1 and LES. On the other hand, several researchers only used k À 1 model to simulate residence time distribution in ozone contactors with different baffling configurations (Wols et al. 2008;Demirel & Aral 2018b;Goodarzi et al. 2020;Kizilaslan et al. 2020).
Although LES is able to better simulate the disinfectant mixing with the flow, LES requires substantially finer meshes than those used for two-equation turbulence models and a very long flow-time for performing statistics of the flow. Thus, an expensive computational effort is required in LES. Considering a compromise in the accuracy and computational effort, the present study simulates the flow in the ozone contactor using the Reynolds-averaged governing equations. The effect of two-equation turbulence models on the residence time distribution is investigated in detail.

METHODOLOGY Flow equations
Ozone contactors are usually equipped with multiple baffles to direct flow. The simple plug flow model is not capable of describing this kind of flow. CFD is a powerful tool to simulate flow in ozone contactors. Considering the steady flow, the movement of water can be described by the continuity equation and the momentum equations as follows where r is the density, p is the pressure, x j is the coordinate component, u j is the velocity component, m is the dynamic viscosity, Àru 0 i u 0 j is the turbulent stress.

Eddy viscosity model
Following the Boussinesq assumption, the turbulent stress can be expressed as follows where m t is turbulent viscosity, k is the turbulent kinetic energy, 1 is the dissipation rate of the turbulent kinetic energy. m t should be determined by a suitable turbulence model.
Standard k À 1 model The Standard k À 1 model has been widely used to simulate the turbulent flow. The governing equations of the turbulent kinetic energy k and its dissipation rate of turbulence 1 are written as follows: where G k is the generation of turbulent kinetic energy due to the averaged velocity gradient, expressed as follows: The model constants C m , C 11 , C 21 , s k and s 1 take the values of 0.09, 1.44, 1.92, 1.0, 1.3, respectively.
RNG k À 1 model The RNG k À 1 model is derived based on the renormalization group theory. Similar to the Standard k À 1 model, the governing equations of the turbulent kinetic energy k and its dissipation rate of turbulence 1 are solved, which is written as follows: where a k and a 1 are the reciprocals of the effective turbulent Prandtl number of the kinetic energy k and its dissipation rate 1, R 1 is an additional term related to the rapidly strained flows. The model constants C m , C 11 , C 21 , a k and a 1 take the values of 0.0845, 1.42, 1.68, 1.393, 1.393, respectively.
Realizable k À 1 model The Realizable k À 1 model also solves the turbulent kinetic energy and its dissipation rate. However, the turbulent viscosity and the governing equation of the dissipation rate are modified to satisfy mathematical constraints on the normal Reynolds stresses. The governing equations of the turbulent kinetic energy and its dissipation rate of turbulence 1 are written as follows: The model constants C 1 , C 2 , s k and s 1 take the values of max 0:43, h h þ 5 , 1.9, 1.0, 1.2, respectively, where The standard k À v model can reproduce low-Reynolds number effects and shear flow spreading, however, it is of freestream sensitivity. The SST k À v model is a combination of the standard k À v model and k À 1 model. That is, the standard k À v model is applied near the wall whereas the k À 1 model is applied away from the wall. The transport equations of k and v are expressed as follows: with where Y k and Y v are dissipation term of k and v, respectively, D v is the cross-diffusion term, F 1 is the blending function. The model constants s k,1 , s k,2 , s v,1 and s v,2 take the values of 1.176, 1.0, 2.0, 1.168, respectively.

Governing equations for tracer
In order to analyze hydraulic efficiency, a numerical tracer study is conducted by solving the Reynolds averaged advection diffusion equation for passive tracer concentration as follows where C is the Reynolds-averaged tracer concentration.

PHYSICAL MODEL AND NUMERICAL SETTINGS
A three-dimensional ozone contactor is shown in Figure 1. The ozone contactor has a cross-section of 230 mm-Â240 mm and a length of 480 mm along the flow direction. Three baffles with the dimension of 180 mmÂ230 mm evenly distribute along the vertical direction with an interval of 120 mm, dividing the ozone contactor into four chambers with the same dimension of 230 mmÂ240 mmÂ120 mm. The inlet has a dimension of 230 mmÂ30 mm. Hexahedral mesh is used to discretize the computational domain, as shown in Figure 2. Three sets of grid are tested for grid sensitivity, namely grid 1 (1,953,640 cells), grid 2 (2,381,700 cells) and grid 3 (2,782,330 cells). Figure 3 shows the velocity and the turbulent viscosity computed using the standard k À 1 model along the line (x ¼ 0.3 m, y ¼ 0.115 m) for three sets of grid. These two parameters have a significant effect on the tracer transport such that they are selected for grid sensitivity. The nearly identical velocity profiles are obtained for Uncorrected Proof three sets of grid. Turbulent viscosity computed using grid 2 and grid 3 are also nearly the same. Thus, the present 2,381,700 cells can achieve satisfactory results. No-slip condition is imposed at all walls. Standard wall function is used for turbulence near walls, which writes for y Ã . 11:225 where U Ã ¼ u P C 1=4 m k 1=2 P t w =r is the dimensionless velocity at the wall-adjacent cell centroid P, y Ã ¼ ry P C 1=4 m k 1=2 P m is the dimensionless distance from the wall, k is the von Kármán constant. When y Ã , 11:225 exists, U Ã ¼ y Ã is assumed (Ansys 2020). The inlet velocity is 0.029 m s À1 , the turbulence intensity is 0.08, and the turbulent viscosity coefficient ratio is 0.1. Water surface is set as a symmetry. The pressure is specified at the outlet. At the inlet the concentration of tracer is set as C 0 ¼ 1. Thus, the variation of tracer concentration with time at the outlet is just the cumulative residence time distribution F(u). The residence time distribution ay the outlet can be obtained by ANSYS Fluent is used to solve governing equations. The pressure-velocity coupling is treated using the SIMPLE algorithm. All convection terms are discretized with second-order upwind differencing scheme. When the residuals for all equations except tracer concentration are less than 1:0 Â 10 À4 , the computations is assumed convergent. The residual criteria for tracer concentration is 1:0 Â 10 À6 . Compared to the standard k À 1 model, the computational efforts of RNG k À 1 model, Realizable k À 1 model and SST k À v model are 124, 140 and 160%, respectively. This is due to additional nonlinearities of these three models.  Figure 4. The experimental data from Kim et al. (2010a) is also included. All turbulence models capture three peaks with decreasing magnitudes on the curves of residence time distribution. The number of peaks are in accordance with that observed in the experiment. However, the magnitude of peaks for all turbulence models are greater than those experimental values, especially in the first peak. Accordingly, the computed residence time distribution is less than the experimental values after the normalized time is greater than 0.5. The initial peaks are 1.53, 2.01, 2.29, 2.54 for the Standard k À 1 model, the RNG k À 1 model, the Realizable k À 1 model and the SST k À v model, respectively, which is greater than the experimental value of 0.99. Among four turbulence models, the Standard k À 1 model obtains the lowest peak on the curve of residence time distribution. Correspondingly, the curve for the Standard k À 1 model is closer to the experimental curve than the other three turbulence models. Compared with the Standard k À 1 model, the first and second peaks of the other three turbulence models shift to the left. The phenomenon demonstrates the existence of a short stream in the ozone contactor. However, the predicted values of the valley in the RTD curve by the RNG k À 1 model, the Realizable k À 1 model and the SST k À v model are very close to the experimental value.
The cumulative residence time distributions for four turbulence models are shown in Figure 5. The cumulative residence time distribution curves are obtained by integrating the residence time distribution curves. All turbulence models predict a relatively bigger values of cumulative residence time than the experimental values. This discrepancy is obvious in the initial period of the experiment. As time proceeds, all predicted curves tends close to the experimental values. In general, all turbulence models reproduce the cumulative residence time. Among four turbulence models, the cumulative time distribution of the Standard k À 1 model is in better agreement with the experimental values. Therefore, the standard k À 1 model has the best simulation effect on the internal water flow in the ozone contactor.
From the RTD curve, the averaged mean residence time at the outlet can be calculated as follows The experimental value is 111.25 s while the computed values are 109.87 s, 109.39 s, 108.52 s and 107.42 s for the Standard k À 1 model, the RNG k À 1 model, the Realizable k À 1 model and the SST k À v model,  The comparisons of dimensionless residence time for four turbulence models are shown in Table 1. All values of t 10 , t 50 , t 90 , and MDI for the Standard k À 1 model have the smallest error, showing the good capability in the simulation of the internal flow of the ozone contactor.

Analysis of velocity field
The velocity field at the mid-plane (y ¼ 0:115m) is shown in Figure 6. Multiple re-circulation regions exist in the ozone contactor. In each compartment a big re-circulation region and a small one are formed. Water moves mainly along a curved flow passage, as indicated by green color in Figure 6. This curved passage determines the shape of residence time distribution at the outlet. All four turbulence models predicted the above feature. The difference lies the maximum velocity in the curved flow passage. The standard k À 1 model gives the maximum value of the maximum velocity, as shown in Figure 6(c). The patterns of the re-circulation region in the first compartment are similar for four turbulence models. The standard and RNG k À 1 models predict the similar  Uncorrected Proof pattern of the re-circulation region in other three compartments. However, the patters of the re-circulation region in other three compartments are totally different for Realizable k À 1 and SST k À v models. Especially, the flow velocity in the fourth compartment predicted by SST k À v model is less than other three models. Thus, SST k À v model gives the biggest magnitude of the first peak in Residence time distribution. In general, velocity range predicted by the standard k À 1 model are less than that by the other three turbulence models. That is, the tracer under the other three turbulence models can get to the outlet earlier, thus the first peaks in the RTD curve shift to left.

Analysis of turbulent viscosity
The turbulent viscosity at the mid-plane (y ¼ 0.115 m) is shown in Figure 7. The turbulent viscosity predicted by the Standard and RNG k À 1 models are similar in shape. However, the magnitude predicted by the Standard k À 1 models is the largest, so its hydraulic residence time is also the smallest. The patterns of the turbulent viscosity computed by Realizable k À 1 and SST k À v models are different from those computed by the Standard and RNG k À 1 models. In general, the magnitudes computed by RNG and Realizable k-1 models and SST k À v model are close.

CONCLUSION
Flow and tracer transport in the ozone contactor is conducted by using CFD. Residence time distribution is calculated based on the tracer concentration. All turbulence models capture the feature of residence time distribution and cumulative residence time distribution. The discrepancy between the experiment and all turbulence models is observed in the initial period. The standard k-1 model predicts the smallest peak and shifts the position of peak to right. The other three turbulence models obtain the similar position of peak whereas predict a more bigger values of peak. All turbulence models reproduce the mean residence time at the outlet. In terms of velocity and turbulent viscosity in the ozone contactor, multiple re-circulating region and dead zones exist in the reactor. Compared with the other three turbulence models, the standard k-1 model predicts