CFD-PBM Coupled modeling of bubble size distribution in a swirling-flow nanobubble generator

ABSTRACT Swirling-flow nanobubble generator is an efficient hydrodynamic cavitation method that can continuously produce enormous quantities of bulk nanobubbles. However, the development of hydrodynamic models of nanobubble generation remains challenging and is rarely found in the literature due to the associated modeling complexity. In this work, a hydrodynamic model was developed to predict bubble size distribution in a swirling-flow type nanobubble generator using a combination of computational fluid dynamics (CFD) and population balance method (PBM). The proposed model was evaluated by considering several combinations of bubble coalescence, breakage, and turbulence models. The results show that the combination of the turbulence coalescence and Luo breakage models predicted better than any other model combination. The selection of appropriate turbulence models could improve the modeling accuracy. The standard k-Ω model provided better predictions than other turbulence models for high flow rates, while the standard k-ε model was more appropriate for low flow rates. The bubble number density was successfully predicted for a 30 min generation time. The turbulence dissipation rate influenced the bubble number density and mass transfer, and the results of simulation considering this relation corresponded to the experimental results. Therefore, the coupled CFD-PBM model can aid in the design of nanobubble generators using virtual prototypes and reduce the development costs required for broader applications.


Introduction
Nanobubble technology has attracted the attention of researchers because it can benefit the discovery and development of new and advanced technologies in several fields. Nanobubbles are nanoscale bubbles dispersed in a liquid and collectively known as bulk nanobubbles. Nanobubbles are specified as bubbles having a diameter of 200 nm or less (Agarwal et al., 2011). Meanwhile, Temesgen et al. (2017) classify bulk nanobubbles as tiny bubbles with diameters less than 1 μm. The international standard for fine bubble technology has established the latter size classification (ISO, 2017). Nanobubbles are distinguished from regular bubbles by their longterm stability and high gas solubility in liquid (Ebina et al., 2013;Nirmalkar et al., 2018;Ohgaki et al., 2010;Weijs et al., 2012). The potential applications of nanobubbles have been extensively explored in the literature as a result of these properties. Nanobubbles have been employed in mineral processing and mining for flotation and particle separation (Ahmadi et al., 2014;Calgaroto et al., 2015;. Nanobubbles have also been used in water treatment, disinfection, and cleaning technologies (Agarwal et al., 2011;Imaizumi et al., 2018;Wu et al., 2008). Nanobubbles have been employed to fabricate high-performance batteries and hydrogen storage (Lim et al., 2019;Xia et al., 2018). Nanobubbles have been used in energy systems to enhance fuel characteristics, engine performance, and biodiesel generation (Ahmad et al., 2019;Nakatake et al., 2013;Oh et al., 2013). Nanobubbles can also be used in agriculture to boost plant and fish bioactivity and growth rates Mahasri et al., 2018).
Studies on nanobubble generation are frequently described in the literature, including patents and articles. However, most investigations are conducted on limited experimental sizes using non-continuous generating systems that may be challenging to scale up. Hydrodynamic cavitation is the most extensively utilized method for producing bulk nanobubbles. A swirling-flow nanobubble generator is a high-efficiency hydrodynamic cavitation technique for making large amounts of nanobubbles. Despite using a modest gas flow rate, the swirling-flow nanobubble generator has the highest mass transfer efficiency compared to venturi, ejector, and pressure dissolution techniques (H. Li et al., 2013;Terasaka et al., 2011). However, developing hydrodynamic models of nanobubble generation is still challenging and very rarely reported in the literature due to the associated modeling complexity. Nanobubble generator models based on numerical methods can explain the mechanisms that occur during the generation process, facilitating nanobubble generators with high efficiency, especially by optimizing the mass transfer efficiency, flow rate, and bubble size. The purpose of this study is to develop a hydrodynamic model that can predict bubble size distribution in a swirling flow type nanobubble generator; thus, the results can assist in the design of nanobubble generators and minimize the development costs for broader applications.
The two main frameworks used in two-phase flow simulation are Eulerian-Lagrange and Eulerian-Eulerian. Several studies have attempted to model the bubble generation process at the microscale and nanoscale using CFD (computational fluid dynamics) methods. Basso et al. (2018) modeled the generation of microbubbles from a venturi nozzle using CFD based on the Eulerian-Lagrange approach. Bubble dynamics can be accurately captured in this way, but in many cases, interactions between bubbles are generally neglected in these methods because their inclusion requires a high number of grids and high computational costs (Long et al., 2018;Maeda & Colonius, 2018). Compared to the Eulerian-Lagrange, the Eulerian-Eulerian method does not require solving the balance equations for individual bubbles, so computational resources are saved (Li, Yu, et al., 2019). Ren et al. (2019) attempted to model bubble generation at the nanoscopic level on a honeycomb structure through the three-dimensional (3D) CFD modeling using the Eulerian-Eulerian approach coupled with PBM. The CFD-PBM model could predict nanobubble size distributions. In CFD-PBM coupled simulation, the multiphase flow problem in CFD simulation was solved using the transport equations of mass, momentum, and energy. Meanwhile, PBM could estimate bubble evolution and dynamics by considering the effect of aggregation and breakage at the interface region between the continuous fluid and dispersed bubbles (Halfi et al., 2020;Lee, Jung, et al., 2020). The CFD-PBM coupled model can thus simulate fluid flow over the entire domain (Liao, 2020).
In this work, the nanobubble generator was modeled using CFD-PBM coupled simulation. The two-phase flow CFD simulation based on the Eulerian-Eulerian approach involves continuous-liquid fluid flow and bubble dispersion. The liquid phase was modeled based on the mass and momentum conservation equation followed by turbulence modeling, and the gas phase was modeled based on the Eulerian approach using a mixture model. The several bubble coalescence and breakage models and turbulence models were applied in the CFD-PBM simulation to obtain the nanobubble size distribution with the best accuracy.

Device design & experiments
Figure 1(a) shows a swirling-flow nanobubble generator mounted vertically. The apparatus's most innovative feature is a two-chambers nozzle in which the inlet and outlet are oriented on the same axis; thus, the outlet nozzle may be linked to other devices. A swirling flow will be generated in the inner chamber using guide holes of the nozzle. The fluid pressure achieves a minimum in the nozzle axis area; meanwhile, the tangential velocity of the liquid reaches a maximum, resulting in a very modest gas input pressure with low energy consumption. Bulk nanobubbles are produced through hydrodynamic cavitation at the orifice or nozzle outlet area. Bernoulli's equation accurately describes hydrodynamic cavitation as P + 1/2ρU 2 = C (constant) and U 2 + 2P/ρ = 2C/ρ (Munson et al., 2002). In the orifice region, the pressure P will be negative when the liquid flow rate exceeds 2C/ρ. In this area, cavitation inception will occur where small bubbles are formed containing cold vapor or other gases diffused in the liquid. The parameter that describes the condition of cavitation is the cavitation number C n , which is estimated using C n = P a − P v / 1 2 ρU 2 (Eisenberg, 1961). The cavitation number when cavitation inception is called the critical cavitation number, which is a function of geometry, and its value can vary widely. The pressure in the cavity in a two-phase flow with one flow component is solely determined by the vapor pressure. However, in the case of a nanobubble generator that uses multi-component flow, the pressure in the cavity is equal to the sum of the liquid's saturation pressure and the partial pressure of the gas entering the cavity.
This study generated bulk nanobubbles using the recirculation test bench comprised of a centrifugal pump, swirling-flow nozzle, and bubble column, as demonstrated in Figure 1(b). Deionized (DI) water was utilized in the experiment. To ensure the absence of contaminants, the AZ-86555 Benchtop meter was used to measure the DI water's quality with a conductivity of less than 3 ± 0.2 μs·cm −1 . The water temperature was kept at a constant of 20°C ± 3°C using a chiller and heat exchanger in the feed tank. The DI water was circulated throughout the system by a centrifugal pump, which the flow rate was regulated by a variable-frequency drive (VFD) controller. The oxygen was then combined with water in the swirling-flow nozzle. A QLY-5L oxygen generator with a greater than 95% product purity was employed to feed pure oxygen. The oxygen flow was measured with a gas flow meter 0−60 ± 1.8 L/min and a pressure sensor 0−1.2 ± 0.018 MPa. Bulk nanobubbles were formed under a saturated condition in the 7.4 L column. The bubble column was composed of translucent Plexiglas and had an exit pipe linked to the drain and feed tank to produce a cycling system. Experiments were performed using a fixed 1 L/min gas stream with two water-to-gas flow rate ratios (Q l /Q g ), 38 and 20, and three-generation times, 11 s, 15, and 30 min. At the Q l /Q g of 38, the time required to complete one generation cycle in a 7.4 L bubble column is 11 s. Prior studies revealed equilibrium that peak fragmentation in the gas phase occurred after 30 min (Michailidi et al., 2020). The nanobubble water's dissolved oxygen (DO) concentration was measured by the polarographic approach with a measuring range of 0-50 ± 0.35 mg/L using a DO meter (Lutron WA-2017SD). The equipment was cleaned before each set of tests by cycling it 10 min through a chromic acid solution, followed by three 15-min washes with DI water.
The Zetasizer Nano ZS (Malvern Instruments) based on dynamic light scattering (DLS) technique was used to analyze the nanobubble size distribution. DLS can characterize particles, emulsions, and molecules mixed in a liquid solution. Light is dispersed at varying intensities due to the Brownian motion of particles in a suspension (Malvern, 2021). The Stokes-Einstein equation calculated variations in light intensity induced by Brownian motion and converted them into particle sizes (Nirmalkar et al., 2018): According to Rayleigh's equation, the intensity of light scattering is proportional to particle diameter, L 6 . Small particles' motion is quicker, resulting in a more significant intensity of scattered light fluctuation. On the other hand, large particles' motion is sluggish and results in lesser fluctuation intensity (Wang et al., 2019). The Zetasizer Nano ZS was calibrated using the NIST monodisperse particle standard in our previous work resulting in measurement accuracy for nanoparticle size of ± 5.79% based on number intensity distribution (Alam, Sutikno, Soelaiman, & Sugiarto, 2021). In this work, the nanobubble size distribution for each generation time was based on average values of nine measurements, and each measure became the average value of 29 continuous runs. The average diameter of nanobubbles and standard error were obtained using the formulae for combining groups (Higgins & Green, 2008). If there are two sample groups, x i and y i , with sizes n and m respectively, the combined sample mean can be obtained bȳ The birth rate as a result of coalescence The death rate as a result of coalescence The birth rate as a result of breakage The death rate as a result of breakage Coalescence -Turbulent (Abrahamson, 1975;Saffman & Turner, 1956) Breakage -Luo (Luo & Svendsen, 1996) C 2 = 2.52,C 3 = 0.04, andC 4 = 0.01 As a result, the combined sample variance is

Numerical model
The numerical model utilized in this work was based on a two-phase flow study using the Eulerian-Eulerian approach using a mixture model, as described in our previous study (Alam et al., 2020). CFD-PBM equations were utilized to predict the size distribution of bulk nanobubbles by combining coalescence and breakage models with turbulence models. The equations used in PBM modeling to predict the nanobubble size distribution are shown in Table 1. The growth phenomenon was neglected because the nanobubble sizes were stable in the water due to the ion-shielding effect (Dressaire et al., 2008;Uchida et al., 2016). The bubble coalescence phenomenon was considered based on the variation between the turbulence (Abrahamson, 1975;Saffman & Turner, 1956) and Luo (Luo, 1993) models. The lowest number of eddies in the turbulence model is determined based on the Kolmogorov microscale, which relies on the turbulence dissipation rate and kinematic viscosity (Rodrigues et al., 2019). In the Luo model, bubble coalescence is similar to the turbulence model because the turbulent flow is considered. The coalescence probability is influenced by the bubble size, the density ratio of the two-phase fluid, and Weber's number.
The bubble breakage phenomenon was considered by several methods, including the Luo (Luo & Svendsen, 1996), Lakkonen (Laakkonen et al., 2006), and Lehr (Lehr et al., 2002) models. The Luo breakage model comprises the density likelihood function and breakage recurrence after bubble bursting. The model is based on the ideas of isotropic turbulence and probability. Furthermore, the Lehr breakage model is almost identical to the Luo model, but it uses different parameters based on Lehr's experiment. As in the Luo and Lehr models, the Laakkonen model can also model bubbly flow characteristics.

Simulation details
In this work, all simulations used the CFD software Ansys ® Fluent. The procedures of CFD-PBM coupled simulation followed the steps according to the modeling framework in our previous study (Alam et al., 2020).  nodes were used in the preliminary simulation. An inflation layer with a smooth transition option was placed on the 3D surface with a target wall y+ of 1.0 for the bubble column. Geometry meshing was performed using the Fluent solver preference with the program-controlled option, a high smoothing mesh, and a target skewness of 0.9.
In setting boundary conditions, the pressure-based, steady-state, and absolute velocity formulations were used in the simulation with a gravitational acceleration of 9.81 parallel to the nozzle axis (z-axis). The mixture model was selected using two eulerian phases (water and oxygen) with a dispersed interface modeling and implicit body force in the two-phase flow model. The oxygen stream rate was set at 1 L/min at inlet gas, with two cases of Q l /Q g ratios, 38 and 20. Therefore, the gas volume fractions in the vortex chamber for the Q l /Q g of 38 and 20 were 0.050 and 0.026, respectively. In PBM, the number of bubble bins used in nanobubble generation referenced previous research, in which the simulation focused on tracking bubbles smaller than 5,120 nm with a parent diameter of 1,280 nm in the gas inlet (Ren et al., 2018(Ren et al., , 2019. The reason is that the nanobubble size in the input area of the nanobubble generator is difficult to measure, so it remains unknown (Ren et al., 2019). In this study, bubble diameter (bin size) tracing in PBM included 19 bins, and the range and number of each bin are shown in Table 2. Because the Zetasizer Nano ZS has a measuring range of 1 nm to 10 μm, all of the simulation bins were included within the instrument measurement range, allowing for a simple comparison of simulation and experimental results. PBM employed the discrete method by combining several variations of bubble coalescence and breakage kernels, where bubble breakage formulation referred to the Hagesaether model (Hagesaether et al., 2002).
Based on the coalescence and breakage phenomena, the two-phase flow simulation was strongly influenced by turbulent flow. Therefore, the inclusion of a turbulence model in the CFD simulation will have significant effects, resulting in proportional simulation results and accurate predictions of bubble size distributions. In this study, several turbulence models were examined in the twoequation model to determine which produced the more precise nanobubble size distribution. The Reynolds Average Navier-Stokes (RANS) equation formed the basis for the two-equation model. This model incorporates the transport equation for turbulent kinetic energy, k, and turbulent energy dissipation, ε, or specific turbulent energy dissipation, (Nichols, 2008).
In multiphase CFD modeling coupled with PBM and cavitation model development, k-ε models, including the standard k-ε (Sutherland et al., 2019;Zhou et al., 2019), and RNG k-ε (Gemello et al., 2018) models, provide fast and accurate simulation results. The fluid flow is assumed to be under fully turbulent conditions in the standard k-ε model, in which the molecular viscosity effect is negligible (Launder & Spalding, 1972). Meanwhile, the RNG k-ε model is based on a statistical technique from the Renormalization Group (RNG), which improves the flow accuracy of the standard k-ε model (Yakhot & Orszag, 1986). In addition to considering the swirling effect, the RNG k-ε model also considers the derivative equation for the effective viscosity; thus, it can calculate turbulent flows with low Reynolds numbers.
In addition to the k-ε model, k-models are generally applied as the turbulence model in cavitation and CFD-PBM modeling, including the standard k- (Li, Bussonnière, et al., 2019) and the SST (shear stress transport) k- (Deng et al., 2019) models. The k-model incorporates the transport equation for turbulent kinetic energy (k) and the formula for specific turbulent energy dissipation ( ), where is the ratio of k to ε . The standard k-model was derived from the Wilcox turbulence model, which considers the effects of shear flow and compressibility at low Reynolds numbers (Wilcox, 2006). The SST k-model improves the calculation of the turbulent viscosity due to considering the influence of turbulent shear stress, overcoming weaknesses in previous equation models (ANSYS, 2013;Menter, 1994). The four abovementioned turbulence models were compared to yield the highest accuracy modeling results in this study. The RNG k-ε model was used in the preliminary simulation with Q l /Q g ratios of 38 and 20, and the results were analyzed in the inlet and outlet areas of the bubble column.

Mesh independence study
The mesh independence study referred to the procedures in previous studies that effectively determine the best grid resolution based on the grid convergence index (GCI) . GCI is a simple approach for universally reporting grid-convergence investigations that are not confined to integer refinement. The GCI is based on extended Richardson's extrapolation, which compares discrete solutions with various grid spacings (Roache, 1997). GCI can be calculated based on Richardson's extrapolation with : A grid refinement ratio larger than 1.3 was suggested to separate the discrete error from other faults . Therefore, the grid refinement ratio (r) can be calculated as follows : The grid independence test was carried out in this work for three grids, namely fine, medium, and coarse grids with a grid refinement ratio larger than 1.3. The F s value of 2.5 was recommended to compare three or more grids (Hassanien et al., 2017). The order of convergence, p can be calculated using The numerical solution's relative error for each grid condition can be computed using The GCI was used to calculate the numerical solution error based on the mean bubble diameter of each CFD model. The GCI between fine and medium grids, as well as between medium and coarse grids, is determined using Eq. (28). A higher GCI value indicates increased grid error, implying that a more refined grid design is necessary. The GCI of coarse and medium grids, as well as medium and fine grids, were compared to assess grid adequacy.

Results and discussion
The nanobubble size distributions obtained using the Zetasizer Nano ZS for the Q l /Q g ratios of 38 and 20 with three-generation times of 11 s, 15, and 30 min are shown in Figure 3. The error bars in Figure 3(a−d) represent the standard deviation of the nine measurements, each of which comprises 29 runs. The size distributions based on the intensity at the Q l /Q g ratios of 38 and 20 are shown in Figure 3(a,b), respectively. Both Q l /Q g ratios produced a bimodal distribution with a bubble size of less than 1 μm, indicating that the bubbles could be categorized as bulk nanobubbles or ultrafine bubbles. After converting signal intensities to bubble numbers using the Mie theory (Malvern, 2021), the majority of the bubble sizes were less than 200 nm at Q l /Q g ratios of 38 and 20, as shown in Figure 3(c,d), respectively. As mentioned in Section 2, the signals of large particles are more potent than those of tiny particles because the intensity of scattered light is proportional to L 6 . In this work, the numerical simulations predicted the nanobubble sizes based on the bubble number density, thus, bubble size prediction was compared with experimental results using the number distribution. The error bars of 11s, 15min, and 30min overlap each other, indicating that the generation time has no significant effect on the size of the nanobubbles due to a low degree of bubble coalescence. However, it has an effect on increasing the bubble concentration as reported in the previous study (Etchepare, Oliveira, Nicknig, et al., 2017). If the measurement results for the three-generation times are combined, the bubble size distributions for the Q l /Q g ratios of 38 and 20 are shown in Figure 3(e,f), respectively. Using Eq. (2) and (3), the mean bubble diameters based on the sum of the three-generation times at Q l /Q g 38 and 20 are 82.77±58.47 nm and 84.73±59.64 nm, respectively. The greater the Q l /Q g ratio, the higher will be the inlet water pressure of the swirling-flow nozzle. If we define the Q l /Q g ratio as x and the input water pressure as y, then the relationship between the inlet water pressure and Q l /Q g is y = 11.219x -91.663, with the coefficient of determination, R 2 = 0.9868. Variation in the Q l /Q g ratio affects the solubility of the gas in the liquid, which can be proven by examining the DO concentration. Figure 4 shows the DO concentrations after a 30-minute generation time with Q l /Q g ratios of 38 and 20. At a Q l /Q g of 38, the highest DO concentration was 32.4 mg/L, while at a Q l /Q g of 20, it was 27.14 mg/L. The DO concentrations were comparable to those achieved in a recent work, yielding a 31.7 mg/L DO concentration (Ebina et al., 2013). After 30 min of generation, there was no more increase in DO concentration. The optimum generation time is predicted to be approximately 30 min when the highest fragmentation of the gas phase occurs, as previously reported in a study (Michailidi et al., 2020). The solubility of a gas in water satisfies Henry's law, based on the dissolved gas's partial pressure (Takahashi & Miyahara, 1979). The DO concentration is proportional to the interfacial area of the bubble, which is determined by the bubble size and concentration. The higher the DO concentration, the smaller the bubble size and the higher the bubble concentration.
The preliminary simulation results of PBM are shown in Figure 5. The simulation included the turbulent RNG k-ε, the Luo breakage model, and the turbulence coalescence model at a Q l /Q g of 38. Figure 5(a) compares the initial number density distributions at Q l /Q g ratios of 38 and 20 at the gas inlet region. The size was consistent with the boundary conditions, in which the initial bubble diameter was set as 1,280 nm, or bin 15. The number densities at Q l /Q g ratios of 38 and 20 were equivalent due to the constant gas flow rate of 1 L/min. Figure 5(b) compares the number density at the outlet column for all bins. For both Q l /Q g ratios, the number densities of bubble size larger than 2,000 nm were higher, essentially at the column outlet. Hence, most of the bubbles tended to coalesce into larger bubbles at the column outlet. The large bubbles quickly rose to the water's surface and collapsed;  thus, only nanobubbles remained in the water for DLS detection. Upon limiting the number density at the outlet area to bins with sizes smaller than 1,000 nm (bins 1-14), most of the remaining bubble sizes were smaller than 250 nm, as shown in Figure 5(c). These results corresponded with the DLS results. The bubble concentration was higher at a Q l /Q g of 38 than at a Q l /Q g of 20. This result was consistent with the DO test results, in which the DO concentration was higher at Q l /Q g 38 than at Q l /Q g 20.
Furthermore, mesh independence was studied at Q l /Q g of 38 and 20 using the same boundary conditions and model combinations as the preliminary simulation of PBM. The three grid resolutions used in the test for both cases of Q l /Q g are fine, medium, and coarse grids, with the number of grids being 454,951, 296,159, and 191,598, respectively. The ratio of grid refinement between medium and fine grids (r 2,1 ) and between the coarse and medium grids (r 3,2 ) calculated according to Eq. (29) are both equal to 1.5, subscript 1 denotes the most refined grid. Since the study is focused on predicting bubbles' size, the grid convergence index (GCI) calculation is based on the mean bubble diameter obtained from simulations using the three grid resolutions. Figure 6 shows the results of the grid independence test on the Q l /Q g ratios of 38 and 20. The number density distributions were converted into percentages for bin sizes < 1,000 nm because the experiments did not determine the concentration of bubbles per unit volume. In Figure 6(a), the relative error of the mean bubble diameter at Q l /Q g 38 between the medium and fine grids (e 2,1 ) and between the coarse and medium grids (e 3,2 ) calculated according to Eq. (31) are 0.002% and 0.016%, respectively. Meanwhile in Figure 6(b), e 2,1 and e 3,2 at Q l /Q g 20 are 2.832% and 0.478%, respectively. The order of convergence for both grids ratio p 2,1 and p 3,2 at Q l /Q g 38 calculated according to Eq. (30) are both equal to 1.2. Meanwhile, p 2,1 dan p 3,2 at Q l /Q g 20 are both equal to −0.3. By using Eq. (28) and F s value of 2.5, GCI 2,1 and GCI 3,2 at Q l /Q g 38 are 0.01% and 0.06%, respectively. Meanwhile, GCI 2,1 and GCI 3,2 at Q l /Q g 20 are −61.88% and −10.38%, respectively. In both cases of Q l /Q g , the GCI between medium and coarse grids (GCI 3,2 ) is greater than the GCI between fine and medium grids (GCI 2,1 ), indicating the need for a precise grid design. The relative error between the medium and fine grids (e 2,1 ) at Q l /Q g 38 produces the smallest value. However, the solution of the medium grid at Q l /Q g 20 provides a significant difference compared to the fine grid. The results of the grid independence test for each bin at 38 and 20 Q l /Q g ratios are shown in Figure 6(c,d), respectively. The larger the bin size, the greater the relative error (e) and the grid convergence index (GCI) for both Q l /Q g ratios. However, the number intensities of the bubble size distribution are primarily in the range of bin size less than 200 nm. At a high Q l /Q g ratio (Q l /Q g of 38), the GCI 2,1 and GCI 3,2 did not show a significant difference, with values for both GCIs being less than 3% at bin sizes smaller than 200 nm. On the other hand, the GCI value at a low Q l /Q g ratio (Q l /Q g of 20) showed a significant difference, with the most prominent value being GCI 2,1 , suggesting the need for a higher grid resolution. Since the study aims to produce bulk nanobubbles with a high flow rate, the PBM simulations for the following cases are determined using the most refined grid resolution of 454,951 grids to overcome significant errors in the result of low flow rates. Figure 7 shows the results of PBM modeling were compared by experiments using combined bubble coalescence and breakage models. Based on comparing the PBM simulation and experimental results at Q l /Q g of 38 (Figure 7(a)) and Q l /Q g of 20 (Figure 7(b)), the Turbulence and Luo models combination gives better predictions than other combination models. However, there is a large discrepancy between the simulation and experimental results. Therefore, the modeling needs to be continued by trying several variations of the turbulence model combining Turbulence coalescence and Luo breakage models. The models combined to form each pair had a noticeable influence on the turbulent flow.
In the turbulence coalescence model, the lowest number of eddies is based on the Kolmogorov microscale η (Rodrigues et al., 2019). Therefore, bubble coalescence can be divided into two sub-range mechanisms-viscous and inertial. The bubble size is smaller than the Kolmogorov microscale in the viscous mechanism. This mechanism will affect bubble coalescence by the local shear stress in the vortex, as demonstrated in Eq. (9). In the inertial mechanism, the bubble size exceeds the Kolmogorov microscale. In this mechanism, bubbles will be dragged away by velocity fluctuations in the flow, as demonstrated in Eq. (10). In addition to the influence of the hydrodynamic interactions of the viscous force, bubble coalescence is affected by attractive forces through van der Waals forces, as demonstrated in Eq. (11) and Eq. (12). This study selected the Hamaker constant as 5.31×10 −24 J based on empirical data previously collected for microbubbles smaller than 100 μm (Vakarelski et al., 2010).
Furthermore, the effects of the turbulent eddies on the breakup phenomenon are crucial. In a turbulent field, the arrival or bombardment of eddies causes changes in relative velocity on the bubble surface area, which has amounts of energy for the fragmentation of bubbles. However, only eddies with a length scale less than or equal to the bubble diameter can cause bubble oscillations, while the large eddies mainly cause bubble translatory motion (Luo & Svendsen, 1996). In the Luo model, the breakup rate is calculated using Eq. (19) with the lower limit of the integral, ξ min (ξ min = λ min /L). The incomplete gamma functions can express the integrand, making computation more straightforward. The smallest size of eddies in isotropic turbulence's inertia subrange, λ min , was set as 21.4 times the Kolmogorov microscale, η (η = (μ l /ρ l ) 3/4 /ε 1/4 ) (Tennekes & Lumley, 1972). The higher the turbulent energy dissipation rate, the lower the Kolmogorov microscale and the smaller the size of eddies. Therefore, the smaller eddies have higher breakage frequencies resulting in smaller bubble sizes and higher number densities.
The RANS turbulence models were used to estimate turbulent scales in this study. The selection of appropriate turbulence models is expected to improve the PBM modeling accuracy. Figure 8 shows the bubble size distribution results based on the influence of various turbulence models at Q l /Q g of 38 and 20. The turbulence model significantly affects the size distribution in PBM simulations Figure 6. The results of the grid independence test: (a) and (b) bubble size distribution using the three grid resolutions for 38 and 20 Q l /Q g ratios, respectively, (c) and (d) the relative error and the grid convergence index for each bin at 38 and 20 Q l /Q g ratios, respectively. Subscript 1 denotes the most refined grid.   at Q l /Q g of 38, as shown in Figure 8(a). The size distributions of SST k-, standard k-, and standard k-ε are shifted farther to the right than that of the preliminary turbulence model (RNG k-ε). A similar result can be found at Q l /Q g of 20 that the size distribution of an initial simulation using the RNG k-ε model can be improved by selecting the appropriate turbulence model as shown in Figure 8(b). However, based on the bubble number density aspect, only a few bin sizes from the simulation results are still within the margin of error from the experimental results. Table 3 compares mean bubble sizes between experimental and PBM simulation results using various turbulence models at Q l /Q g ratios of 38 and 20. The experimental results for the three-generation times are averaged to compare with the simulation results, producing mean bubble diameters of 84.73±59.64 nm and 82.77±58.47 nm at the Q l /Q g ratios of 38 and 20, respectively. At a Q l /Q g of 38, the standard k-model provides better predictions than other turbulence models, while at a Q l /Q g of 20, the standard k-ε model can predict better than other turbulence models. Therefore, the standard kmodel is more appropriate for high flow rates, while the standard k-ε model is more suitable for low flow rates.
The bubble number density distributions at Q l /Q g ratios of 38 and 20 when the standard k-model was applied later in PBM modeling using a combination of the turbulence and Luo models can be seen in Figure 9. The simulation used the bin size closest to the experimental results (bin 7, or an 80 nm diameter). Nanobubbles with a diameter of 80 nm at both Q l /Q g ratios began formation with high concentrations in the orifice area or nozzle outlet, and the bubble number at the column outlet subsequently decreased, indicating that some bubbles coalesced into larger bubbles and burst on the surface. At a Q l /Q g of 38 with a nozzle inlet pressure of 334.66 kPa, the maximum bubble number density in the column outlet area for a bubble size of 80 nm was 1.79×10 10 m −3 , or 1.79×10 4 mL −1 , in the 7.4 L bubble column with one generation cycle of 11 s. According to Tanaka et al. (2020), the concentration of nanobubbles will continue to increase linearly with the generation time. According to Etchepare, Oliveira, Nicknig, et al. (2017), nanobubbles resist shear forces from pump impellers as well as pressures of up to 5 bar during multiple generations. Therefore, the bubble concentration with a total generation time of 30 min was estimated as 2.93×10 6 mL −1 , corresponding to the results presented in earlier studies (Ebina et al., 2013;Ren et al., 2019).
The bubble concentration was lower at a Q l /Q g of 20 than at a Q l /Q g of 38 because the turbulent energy dissipation rate was also lower, as demonstrated in Figure 10.  The strength of the turbulent flow affects mass transfer, which is determined by the turbulent energy dissipation (Lamont & Scott, 1970;Prasher & Wills, 1973). Likewise, the generation or exposure time is influenced by vortexes or turbulent flow at the microscopic scale (Garcia-ochoa & Gomez, 2009). These characteristics prove that the solubility is higher at a Q l /Q g of 38 than at a Q l /Q g of 20, corresponding to the measured DO concentrations.

Conclusions
CFD-PBM coupled simulations can predict bubble size distributions on a nanoscale. CFD-PBM simulations evaluated included several coalescence and breakage model combinations. The turbulence coalescence and Luo breakage models can predict better than all other combinations of models. The turbulence model considers hydrodynamic interactions of viscous and van der Waals forces based on the mechanisms of viscosity and inertia. Meanwhile, the Luo model is affected by turbulence fluctuation and bubble surface instability. PBM simulations using various turbulence models were evaluated to improve the simulation accuracy. The results demonstrated that the standard k-model could predict better than other turbulence models for a high flow rate at Q l /Q g of 38, while the standard k-ε model was more appropriate for a low flow rate at Q l /Q g of 20. The bubble number density of 2.93×10 6 mL −1 was successfully predicted at the highest Q l /Q g for a 30 min generation time. The turbulence dissipation rate significantly affects the bubble number density and gas-liquid mass transfer. The higher the turbulence dissipation rate, the higher the bubble concentration and the gas-liquid mass transfer rate. This result corresponded with the experimental results obtained for the dissolved gas concentration. Therefore, the coupled CFD-PBM results can aid in the design of nanobubble generators using virtual prototypes and reduce the development costs required for broader applications. Since the DLS measurements did not determine the concentration of bubbles per unit volume, more research is needed to verify nanobubble concentration using a nanoparticle tracking analyzer (NTA) and improve modeling accuracy based on bubble number density. Further studies are also required to assess the effect of altering geometry and process parameters on the simulation results in producing an optimal design.