Abstract
A novel hybrid multiphase flow solver has been used to conduct simulations of a vertical plunging liquid jet. This solver combines a multifluid methodology with selective interface sharpening to enable simulation of both the initial jet impingement and the long-time entrained bubble plume phenomena. Models are implemented for variable bubble size capturing and dynamic switching of interface sharpened regions to capture transitions between the initially fully segregated flow types into the dispersed bubbly flow regime. It was found that the solver was able to capture the salient features of the flow phenomena under study and areas for quantitative improvement have been explored and identified. In particular, a population balance approach is employed and detailed calibration of the underlying models with experimental data is required to enable quantitative prediction of bubble size and distribution to capture the transition between segregated and dispersed flow types with greater fidelity.
1. Introduction
While computational fluid dynamics (CFD) methods have advanced in many areas such that predictive modeling is possible (e.g., external aerodynamics), multiphase flows remain an area where prediction is yet out of reach for the majority of applications. Multiphase flows present a unique challenge due to diversity of flow regimes characterized by a broad range of scales—dispersed bubbles/droplets at the microscale up to macroscale free surface flows. A compounding complication is the regime dependency limitation of existing multiphase CFD methodologies; specific simulation methods are applicable only to a specific class of flows. For example, free surface flows in which the dynamics of a fluid-fluid interface is important to the overall physics under investigation are typically modeled with a shared momentum equation and a phase fraction field using a sharp interface capturing method such as Volume of Fluid and level set or the like where fluid interfaces are resolved on or by, in the case of interface tracking methods, the computational mesh. However, such methods are infeasible for simulation of dispersed flows characterized by many small fluid particles (bubbles/droplets) which cannot be fully resolved on a computation mesh. In such cases, a so-called multifluid or Euler-Euler approach is taken in which phases are treated as interpenetrating continua, each governed by its own momentum equation and having exchange terms to account for interphase momentum transfer (drag, lift, virtual mass, etc.). This methodological segregation of flow regimes is acceptable for many classes of flows but presents severe limitations for problems in which both segregated and dispersed flow features are present or where transitions between flow types are critical to the phenomena of interest. Recently, a hybrid CFD solver has been developed which aims to overcome the issue of regime dependency by combining the Euler-Euler multifluid method with sharp interface capturing in a regime flexible framework [1].
Many real multiphase flows are highly turbulent and characterized by the generation of very small dispersed droplets/bubbles in a second continuous carrier phase. For such flows, the fluid particle (droplet/bubble) size and its distribution are a critical factor for important phenomena such as interphase heat, mass, and momentum transport processes. Population balance modeling is the most common approach tackling this problem [2]. In such a method, the distribution of fluid particle sizes is typically captured through transport equations for size classes or moments [3] with exchange terms for transfer between classes due to breakage and coalescence. Simplified methods for capturing only a single mean diameter are also available [4]. While the underlying models for breakup and coalescence rates required by population balance based methods are largely empirical and a general, transferable model set is not available [5–7], parametric tuning versus experimental data can produce useful insights [8]. It is also possible to capture the local maximum stable fluid particle size through hydrodynamic theory of breakup based on the Weber number [9].
With the addition of fluid particle size capturing using a number density transport approach, the hybrid multifluid–VOF CFD method is given even greater flexibility [10]. Along with the potential for improved fidelity in the prediction of polydispersed flow fields, it is also possible to use the local prediction of droplet size to guide the switching between dispersed flow modeled regions and interface sharpened regions. Such a method which makes a more direct connection between the local drop size and the ability of the mesh to resolve the interfacial features should have advantages over alternate methods which have been employed elsewhere [11, 12]. The work presented here builds upon that reported previously in Wardle and Weller [1] through the addition of such a novel dynamic algorithm to switch between VOF-resolved and multifluid phase capturing. This new capability adds greater flexibility to the original hybrid method and allows for capturing of flow regime transitions.
The goal of this present publication is to provide an evaluation of this new regime flexible multiphase simulation framework through application to a specific flow phenomenon—the vertical plunging liquid jet—which encompasses both dispersed and segregated flow features. The vertical plunging liquid jet case, despite a simple configuration in many respects, demonstrates the deceptive complexity of multiphase flows. At short times as the jet plunges into the quiescent liquid pool, the flow features are dominated by free surface flows and large cavity structures which can realistically be resolved using sharp interface methods. At longer times, however, the entrained gas cavities break up and a bubbly plume forms [13]. The vertical plunging jet configuration also has industrial application for bottle and mold filling processes where careful consideration of entrainment and gas holdup is required.
Given the relevance of this flow configuration, a number of researchers have studied vertical plunging jets by both experimental and computational means. Zidouni Kendil et al. [14] performed PIV and numerical studies looking at the characteristics of the velocity field below the free surface. For numerical simulations, they treated the plunging liquid jet as a two-phase problem and neglected flow phenomena above the pool surface. Qu et al. [13] computationally studied the influence of jet velocity and variations of jet falling length on the jet penetration depth using two different approaches: multiphase mixture model approach and level set method. On comparing the predictions from both models with experiment (high-speed video), the level set method was more accurate in the areas of the deformation of the water surface and gas entrainment. High resolution large eddy simulations (LES) have also been conducted using a Volume-of-Fluid style solver within the open-source CFD toolkit OpenFOAM (interFoam) [15]. Excellent agreement between simulation and experiment was found for the dynamics and pinch-off time for the initial cavity following impingement of the jet into the pool. The simulation behavior at longer times was not reported. Schmidtke and Lucas [16] investigated the role of different drag models on the gas void fraction below free surface. They considered two approaches in their work. In the first, water was treated as a continuous phase and gas as a dispersed phase everywhere in the domain while, in the second, the different morphologies of the gas, that is, continuous above and bubbly below the free surface, were taken into consideration and the problem was treated as a three-phase flow.
Of particular relevance to the work reported here, Hänsch et al. [17] extended the inhomogeneous Multiple Size Group (MUSIG) model [18] in ANSYS CFX software with the addition of a continuous gas phase in an attempt to combine a dispersed and a continuous gas phase in one computational domain. Their approach was based on the implementation of an additional interface sharpening algorithm within the Eulerian multifield framework and was also validated on a plunging liquid jet case. Similar to the work of Yan and Che [19], the dispersed gas and continuous gas were treated as separate fields with unique momentum equations within the framework. In order to transition between dispersed and continuous gas fields, a “clustering method” was introduced which is aimed to counteract dispersion of the sharp interface; this is conceptually similar to the compression velocity employed in the OpenFOAM interFoam solver scheme [20] except that in the case of the work of Hänsch et al. the force was employed on the gas phase momentum equation and induced physical clustering of dispersed gas until a critical phase fraction of 0.3 was reached and a transition to complete coalescence to continuous gas was enforced. The physical basis for this clustering force is not clear.
In the current work, we employ the multiphaseEulerFoam solver described in Wardle and Weller [1] along with a local resolution approach to transitions between dispersed and segregated regimes. The advantages and limitations of this methodology will be evaluated in relation to the vertical plunging liquid jet phenomena.
2. Computational Methods
Only a general overview of the hybrid multiphase flow simulation methodology will be given here as the details are reported elsewhere [1]. At a conceptual level, on top of the framework of an -phase Eulerian multifluid solution framework (one momentum equation for each phase with interphase momentum coupling terms), a Volume-of-Fluid (VOF) sharp interface capturing algorithm is applied for desired phase pairs. This is done through addition of interface compression to the volume fraction transport equations and the use of limiters to maintain boundedness of phase fractions and their sum. This solver has been developed using the open-source CFD toolkit OpenFOAM and has been included in the release of OpenFOAM as multiphaseEulerFoam. Version 2.2.1 of OpenFOAM and the corresponding solver is used in this work. The base solver will be described in Section 2.1 and the expanded solver will be described in Section 2.2.
2.1. Hybrid Multifluid–VOF Solver
The multifluid model equations for incompressible, isothermal flow are given by sets of mass and momentum equations for each phase :where the density, phase fraction, and velocity for phase are given by , , and , respectively, and is gravity. The two interfacial forces included are the interphase momentum transfer or drag force and the surface tension force . As described elsewhere [1], the drag force coefficient is taken from the model of Schiller and Naumann [21] and the continuum surface force model of Brackbill et al. [22] is used for surface tension.
The volume fraction transport equation for phase modified with the interface compression term is given bywhere the velocity is applied in a direction normal to the interface to compress the volume fraction field and maintain a sharp interface. The term ensures the modification is only active in the interface region. In addition, a bounded differencing scheme is employed for discretization of (2). The value for the artificial interface “compression velocity” is given byThe term gives the interface unit normal vector for the direction of the applied compression velocity. The magnitude of the velocity is used as dispersion of the interface, which is being counteracted by the compression velocity, can only occur as fast as the magnitude of the local velocity in the worst case. The coefficient is the primary means of controlling the interfacial compression and can be considered simply a binary coefficient which switches interface sharpening ON or OFF (0). With set to 0 for a given phase pair, there is no imposed interface compression resulting in phase dispersion according to the multifluid model. Conversely when it is set to 1, sharp interface capturing is applied and VOF-style phase fraction capturing occurs (forcing interface resolution on the mesh). The interface compression coefficient can be defined and applied independently for all phase pairs. Thus, a sharp interface can be maintained at all interfaces between specific phases (e.g., air–water and air–oil) and dispersed phase modeling with no interface compression can be used for other phase pairs (e.g., water–oil). In this current work, the coefficient has been extended as a spatially varying volumetric scalar field such that it can be used to control the model behavior locally as described in the following section.
A conceptually similar hybrid multiphase approach has recently been reported by Marschall and Hinrichsen [23]. Their work focuses on gas-liquid two-phase flows for which they provide a rigorous derivation for the method including additional interaction terms and models. While we present here an example for gas-liquid flow, the method used here is intended for more general classes of multiphase flows beyond just two phases.
2.2. Dynamic Switching Hybrid Solver
As part of its ongoing in-house development and in order to facilitate this work, a number of enhancements have been made to increase the flexibility of the base solver and enable its applicability to flow regime transitions (dispersed to segregate and vice versa).
While the original solver allows for independent assignment of for each phase pair, it remains a constant value in both space and time. In the expanded solver the value for each phase pair has been implemented as a spatially and temporally varying volumetric field variable. Additionally, in order to provide consistency in the application of interphase forces, the drag force has been scaled by such that it is only applied in regions where the = 0 and dispersed flow exists. The surface tension force has also similarly been scaled by such that it is applied under interface sharpened conditions only.
2.2.1. Reduced Population Balance
As described in Wardle [10], a reduced population balance model (also referred to as a number density transport approach) has also been implemented according to the method of Attarakih et al. [4] and Drumm et al. [24]. This method considers only a single moment of the droplet distribution and thus reduces the population balance down to two quantities: the total number (or number density) and the total volume. Since the volume fraction equation for the dispersed phase is already solved only one additional transport equation is required thus limiting the added computational burden. The particle number density is related to the particle mean mass diameter according toThe transport equation for is given bywhere the source term is a straightforward function of the droplet size-dependent breakage () and aggregation () rates:where is the number of daughter droplets which for the typical assumption of binary breakup is given a value of 2. While a number of breakup and coalescence models have been implemented in the solver, a mixed breakup/coalescence model set will be used for the work reported here (following the work of Drumm et al. [24]) with the breakup rate kernel of Martínez-Bazán et al. [25] and the coalescence rate kernel of Prince and Blanch [26]. Since no experimental data for bubble size are available for the test case simulated here, these correlations are chosen as representative models which are sufficient to demonstrate the methodology.
In addition to the number density transport method, an additional algebraic diameter model has also been implemented for comparison. This model is based on the well-known hydrodynamic analyses of Hinze [9] in which the maximum stable droplet size can be determined by a critical Weber number defined aswhere is the continuous phase density, is the interfacial tension and the fluctuating velocity for eddies of size , and , where , are directly related to the turbulence dissipation rate . As done by Sano et al. [27], the first definition is used and the instantaneous slip velocity between the dispersed and continuous phases () was used as an approximation of . The value for the turbulence dissipation rate is the local value computed by the subgrid scale turbulence model as described later on. From a fundamental perspective, it is not clear whether the use of the critical Weber number model in this way is rigorously consistent as it is intended to describe only how the maximum diameter scales with dissipation rate. Nonetheless, it provides a simpler method with only a single adjustable parameter and is presented here for comparison.
The inception of spurious interfacial currents is a well-known drawback of sharp interface methods such as the interface compression scheme used here [28]. The currents arise from imprecise calculation of the interface normal and consequently the interface curvature . Improved methods for calculation of the interface curvature have been reported in the literature [29]. It has also been reported that simply smoothing of the calculated interface normal field can significantly reduce these spurious currents. We have applied an iterative smoothing approach on the interface normal calculation as outlined by Raeini et al. [30]. Additional manipulation of the field itself has also been done as noted in the following section.
2.2.2. Switching Algorithm(s)
Two different switching algorithms have been implemented for local modification of the coefficient field. The first is based on the normalized gradient of the volume fraction similar to what was proposed by Cerne et al. [11] and later used by others [12, 31, 32]. In this formulation, interface sharpening is switched based on comparison of the gradient of the volume fraction with some specified cut-off value . Thus when the interface begins to become unresolved, the sharp interface is deactivated and when it naturally becomes more sharp compression is activated to maintain it in this state. This method can be seen as a sort of “self-fulfilling prophecy” in that where the interface is sharp you keep it sharp and when it is dispersed you let it disperse and it is included only for comparison. Rather than simply using the gradient itself—the magnitude of which would be sensitive to mesh dimensions—we normalize the local gradient by the global maximum value and apply a cutoff value of 0.4 which seems consistent with the original concept presented by Cerne et al. [11].
With prediction of the local mean droplet diameter from the reduced population balance implementation, it is possible to employ a more flexible algorithm which involves comparison of the local mean diameter and the local mesh size. The mesh size can be estimated using the cube-root of the cell volume ) and set to 1 (ON) whenwhere is a user-specified multiplier giving the number mesh cells of size required to resolve a droplet of size on the mesh. The value of must be greater than 1 as it takes several cells to adequately resolve the interface of a droplet/bubble. Some have suggested this value can be as little as 6 [33] while elsewhere significantly greater values up to 32 are reportedly required to accurately resolve interface dynamics [34]. Note that the compressive interface capturing scheme used here is mass conservative unlike the piecewise linear interface construction (PLIC) methods used in both of the aforementioned papers where mass loss due to inaccurate interface representation is of primary concern. Thus, it is expected that values in the lower end of the range could be sufficient.
For unswitching from a sharpened state (from = 1 to = 0) actuation based on the drop size predicted from the reduced population balance has less meaning since interface sharpening seeks to resolve all droplets on the mesh. Consequently, a different switching criterion should be considered. One such method that has been proposed by Hecht et al. [35] is switching based on the calculated interface curvature. The interface curvature of a spherical fluid particle is defined as the inverse of the sphere radius (). If one defines an interface resolution quality (IRQ) factor such thatit can be shown by substituting from (8) that the factor is equivalent to —that is, the number of mesh spacings spanning a given particle having a radius equivalent to the local interface curvature . Conceptually, this means that as the mesh size is decreased for a given particle diameter, the resolution of this feature on the mesh improves and is reflected in an increase in IRQ. Conversely, as the resolved particle size decreases (or local curvature increases), the ability of the mesh to resolve this feature is decreased and IRQ reflects this. Thus, the proper criterion for switching from = 1 to = 0 is given byIn the simulations presented here for this switching algorithm, values of = 3 and = 2 were used for activation and deactivation of interface sharpening, respectively. While these values represent a very aggressive treatment of interface sharpening, they are sufficient to demonstrate the capability of the algorithm.
Given that IRQ-based deactivation of interface sharpening relies on the accurate calculation of the interface curvature , care was taken to ensure that this quantity was treated appropriately in the region away from the interface to avoid spurious switching due to noise in . It was found that scaling by a factor of (which has a value of 1.0 at the interface and tapers to 0.0 elsewhere) was sufficient to ensure that only curvature in the interface regions was considered. Inaccuracies in the curvature calculation in the standard interFoam family of solvers (which includes multiphaseEulerFoam) and the potential for smoothing to improve accuracy have been noted by others [36].
2.3. Computational Parameters and Flow Conditions
The problem to be investigated here involves a jet of water issued from a nozzle with diameter, 6 mm, into a stationary pool of water. The setup is similar to the experiment of Qu et al. [13] but with smaller rectangular dimension as used by Deshpande et al. [15]. The 3D computational domain chosen is a rectangular box with length, 0.1 m, breadth, 0.1 m, and a depth of 0.3 m. The box contains water filled to a depth of 0.2 m so that the nozzle is positioned 0.1 m above the water surface. The setup of the computational domain is illustrated in Figure 1.
The computational domain was discretized using CUBIT meshing software. Two grids defined as standard mesh (S) and refined mesh (R), shown in Figure 2, were used in this study. To obtain the mesh, the geometry was partitioned into two, each with different meshing scheme specification. The submap scheme and pave scheme were specified for the outer and inner block, respectively. The inner block of the standard mesh was further refined (within OpenFOAM) to obtain the so-called “refined” mesh so that the “standard” mesh had a uniform grid size of 1 mm while the refined mesh had a minimum grid size of 0.5 mm in the jet region.
At the start of the simulation, the domain was filled with water up to a depth of 0.2 m. The velocity fields for air and water were each set to 0.0 m/s in all directions in the domain. The pressure was also initialized to zero.
A fully developed velocity profile with a mean axial velocity of 1.5 m/s was specified at the inlet for the liquid velocity and zero gradient for pressure. At the outlet, a mass flow rate equal to the mass flow at the inlet was specified to maintain a constant water level in the domain. A no-slip boundary condition was specified at all walls and an “atmospheric” boundary condition (0.0 gauge pressure) was specified for pressure on the upper surface of the box (see Figure 1).
The fluid properties employed in the simulation are summarized in Table 1. The surface tension coefficient was taken to be 0.072 N/m. In some cases such as those employing the standard multiphaseEulerFoam solver, a constant dispersed bubble diameter, mm or m, was assumed. As noted above, the enhanced solver allows for variable dispersed bubble size based on either the reduced population balance model or the critical Weber number method. In both cases, the diameters were assumed to vary over the range from 0.1 mm (100 m) to 25 mm. The drag coefficient was calculated based on the correlation by Schiller and Naumann [21]. While many correlations are reported in the literature [37], this simple, well-known model is sufficient to demonstrate the new computational method presented here. Additional corrections for bubble swarming and other effects were not considered.
Turbulence was neglected only for the simulations carried out with interFoam. For cases using the hybrid multifluid–VOF solver(s), a large eddy simulation (LES) treatment of turbulence was used with the Smagorinsky subgrid model [38]. While the role of the subgrid scale turbulence model and quantitative prediction of the turbulence dissipation rate is acknowledged to be of great importance to the overall accuracy of the simulation, the primary purpose of this work is to demonstrate the feasibility of the flow regime flexible tool used here. This framework is compatible with other more complex subgrid scale turbulence modeling schemes beyond the one used here.
The governing equations were discretized in time using the backward Euler method. The resulting discretized systems were solved using the Geometric Algebraic Multigrid method for pressure with Gauss-Seidel as the choice for smoothing and Preconditioned Biconjugate Gradient (PBiCG) method with Diagonal Incomplete-LU (DILU) preconditioning for velocity. The absolute tolerance was set to and for the equations of pressure and velocity, respectively.
3. Results and Discussion
As a baseline for comparison, simulations were first conducted using the standard OpenFOAM 2.2.1 release versions of the interFoam and multiphaseEulerFoam solvers followed by simulations using the in-house enhanced dynamic switching version of the solver. Where possible, qualitative comparisons with experimental observations of Qu et al. [13] and others will be made along with comparisons with simulations previously reported for similar configurations in the literature as summarized earlier.
Simulations using the standard solver were conducted for two different values of the dispersed phase (air) bubble size: 2 mm and 500 micron each for both sharpened ( = 1) and unsharpened ( = 0) settings giving a total of four configurations. Simulations using the enhanced solver were conducted for four different test configurations. In each case, simulations were run for 1 s following startup and then results were time averaged over an additional 1 s of flow.
The interFoam results reported here are from simulations on the “refined” mesh while the majority of those for the other two solvers with the exception of the cases reported in Section 3.3.1 were done on the “standard” mesh.
3.1. VOF Solver
InterFoam is a two-phase incompressible flow solver in the OpenFOAM CFD toolkit based on a modified volume of fluid approach in which an “artificial” compression term is added to the volume fraction transport equation as in (2) to reduce the effect of the numerical diffusion of the interface. A detailed description of the interFoam solver as well as its implementation can be found in Deshpande et al.’s paper [36].
Using the simulation configuration described above, the test case was run with the interFoam solver. Snapshots of the initial transient jet impact and entrainment process were recorded to determine the accuracy of the simulation using side-by-side comparison of the 1.5 m/s jet velocity experiment of Qu et al. [13] with the simulation. Results are consistent with those reported by Deshpande et al. [15]. The images shown in Figure 3 show the evolution of the 0.5 volume fraction field iso-surface colored by velocity. For this solver only, the results were obtained neglecting the turbulence of the dispersed phase. As seen in Figure 3, the entire process of impingement and penetration of the initial air pockets is well captured in the simulation.
3.2. Hybrid Multifluid–VOF Solver
Simulations using the hybrid solver, multiphaseEulerFoam, in its standard form as part of OpenFOAM-2.2(.x) library have been conducted for a number of cases to explore the impact of various parameters on the simulation behavior for the plunging jet phenomena. To explore the role of interface sharpening, simulations were done both with sharpening turned on, that is, (VOF behavior), and without resolving the interface, that is, (multifluid behavior). In addition, the influence of the dispersed phase bubble size—which for the standard solver is a required input parameter rather than a predicted quantity—was also explored using two different constant diameters: mm and μm.
3.2.1. Role of Dispersed Phase Droplet Size
Figure 4 shows selected time snapshots of the result of the simulation with imposed interface sharpening ( = 1) for two different diameters: mm and μm at two representative time steps. It can be observed from Figure 4 that the case using a smaller uniform diameter (labeled II in the figure) resulted in more entrained bubbles compared to the case with a larger specified bubble diameter.
This difference can be further observed through plots of the void fraction and axial velocity profiles along a line perpendicular to the jet axis at a vertical height at the midheight of the initial liquid pool as shown in Figures 5(a) and 5(b). As with other plots reported here, these results were time averaged over the period from = 1 s to 2 s. The smaller bubble diameter resulted in greater gas entrainment and a slightly increased centerline velocity (Figure 5(b)) with a plume that is more narrow as can also be seen in Figure 4 for the = 1.2 s conditions. This difference in gas entrainment fraction can also be seen in Figure 6 which gives a plot of this quantity along the vertical direction at the jet centerline.
3.2.2. Role of Interface Sharpening
Results for the case with no interface sharpening () are shown in Figure 7. From the figure, it can be seen that the case with a smaller assumed diameter ( μm) has a longer penetration depth than for the larger diameter. This result is not unexpected as the drag force imposed on bubbles is directly proportional to the bubble diameter. However, in both cases, the overall amount of entrainment and penetration was significantly less than for the cases with interface sharpening.
3.3. Dynamic Switching Hybrid Multifluid–VOF Solver
With the enhanced solver that includes the additional capabilities described in Section 2.2 including variable droplet size and dynamic, local switching interface sharpening, the potential exists for a more complete capturing of the phenomena encountered in the vertical plunging liquid jet configuration. In order to explore the range of options for bubble size prediction and dynamic switching, four different setups were simulated. Two cases used bubble diameter prediction using the reduced population balance and two were done using local diameter prediction from the simple critical Weber number model. For the cases using the reduced population balance model, simulations were conducted using each of the two different switching algorithms: the first based on the normalized gradient of the volume fraction field (“normGrad”) and the second based on the comparison of the local bubble size and the local mesh size (“dropSize”) as described in Section 2.2.2. Preliminary simulations showed the prediction of very small bubbles using the selected breakup and coalescence models (see Section 2.2.1). As a comparison of many and varied breakup and coalescence models available in the literature was deemed to be beyond the scope of this work, simple scaling of the rates was applied for each term. A breakup rate scaling factor of 0.025 and a coalescence scaling factor of 10 were used in all cases. As noted later, selection and tuning of breakup and coalescence models parameters using reliable experimental data are required for predictive capability of these tools. In addition to simulations using the reduced population balance method, simulations were conducted using the critical We number model using two different values 1.2 and 2.5 to capture the range of values reported in the literature. For these cases, switching was done using the dropSize model. In all cases, the simulations were initialized with = 1 everywhere.
Figure 8 shows snapshots of the instantaneous resolved bubble field and dispersed gas distribution on a cross-sectional plane for the time s for the four different cases. In addition, the region of interface sharpening is shown in the inset at the lower right of each image. Time-averages of the gas fraction on a cross-sectional plane for each case are shown in Figure 9. In comparing the two different switching algorithms, it was apparent that the normGrad method (Figure 8(a)) readily transitioned to the unsharpened state and was not reactivated to any significant degree. This resulted in a large spreading of the plume region extending all the way to the bottom of the domain but having a very dispersed gas fraction. Conversely, the dropSize switched cases (Figures 8(b)–8(d)) were found to maintain the sharpened state over the length of the incoming jet and the liquid pool surface. Compared to the rPBE cases which resulted in very small bubble sizes (see Figure 10(a)), the relatively large diameters predicted by the critical We models resulted in smaller regions of unresolved dispersed phase (). As shown in (7), the critical Weber number gives an estimate of the maximum stable droplet size rather than the mean. Future work might consider whether scaling of the predicted diameter from the critical Weber number model by a corresponding ratio of the mean to max diameter would be more appropriate.
(a)
(b)
(c)
(d)
Based on simple estimative analysis of images reported elsewhere [13], the average bubble size in the plume region appears to be approximately 2 mm. Thus, the PBE based models significantly underpredicted the bubble size; note that in the plots for bubble diameter in Figure 10(a) for these models the corresponding values have been scaled by a factor of 10. It is not clear whether this result is due to inaccuracies in the breakup and coalescence models or in the prediction of the subgrid turbulence dissipation rate which is critical to accurate prediction of these rates. Since detailed measurements of bubble size distribution were not available for these test conditions, no attempt was made to perform calibration of the underlying PBE models though such an effort will clearly be necessary for quantitative predictions to be expected. For the critical Weber number diameter models, the bubble diameter was somewhat overpredicted with the smallest bubbles being in the center of the jet and having diameters from 5 to 10 mm. Again, calibration of the We model parameter with experimental bubble size data would be necessary for quantitative predictions.
The jet penetration can be seen by the axial profile of the entrained gas fraction as shown in Figure 11. The experiment reported by Qu et al. [13] for these conditions was found to have a plume penetration depth of approximately 14.1 cm. In that work the plume depth was obtained from image analysis with an extent defined by an arbitrary normalized greyscale value of 0.15 in the time-averaged images. It is not clear to what value of gas fraction in Figure 11 this corresponds. If one assumes the penetration depth occurs at the point where the mean gas fraction drops below 5% then the approximate values for each of the four cases are as follows: rPBE, normGrad 9 cm, rPBE, dropSize 17 cm, We = 2.5, dropSize 16 cm, and We = 1.2, dropSize 16 cm. The dropSize cases all show similar behavior—though slightly greater in penetration than that reported in the experiment—while the case which was switched using the normGrad method of Cerne et al. [11] gave much reduced gas entrainment and penetration depth due to more substantial reliance on multifluid dispersed modeling regions coupled with the excessively small bubble sizes predicted by the rPBE model. A combination of the normGrad switching method with bubble diameter prediction using the critical Weber number model (i.e., We, normGrad) was not done but would be expected to have slightly improved results compared to rPBE, normGrad given that larger diameters would be predicted and therefore some resolved bubbles may be possible.
3.3.1. Refined Mesh
Simulations were run on the refined mesh for two of the cases above—the first being for the reduced PBE case with dropSize switching. Given the computational cost of this simulation (approximately per 120 K CPU hrs per 1 s of flow time) this simulation was only run to s. The second case, which was mapped from the solution of the first case at .9 s, was done using the critical We diameter model with a value of 1.0. This simulation was then run from s out to a final time of 1.25 s.
A series of snapshots showing the evolution of the incoming jet and entrained cavity and bubbles are shown in Figure 12. It was observed that, during the initial impingement phase, the incoming liquid stream was slightly broken up due to a switching off of interface sharpening on the tip of the jet where curvature is very high. Due to the excessively small bubble sizes predicted by the reduced PBE models used here, sharpening was not immediately reactivated. This resulted in some breakup of the initial cavity and deviation from the behavior seen in the experiment and when using pure VOF modeling (Figure 3); that is, the pinch-off of the initial cavity was not observed or occurred very soon after the initial impact. Similar behavior was seen for the standard mesh simulations using the reduced PBE though the amount of breakup was more significant given the reduced resolution of the interface giving a smaller value of the IRQ used for unswitching (9). In order to accurately capture and maintain the sharp interface during this initial transient time, improvements will need to be made to the breakup and coalescence models to give more quantitative prediction of the bubble sizes found in the real physical system.
The simulation results including an inset showing the field at the end time for both bubble diameter treatments are shown in Figure 13. Given the larger bubble size predicted from the critical Weber number model and the refinement of the mesh, more of the domain was consistently being sharpened. Though quantitative comparison of the two results to compare the size and number of the resolved bubbles in the two cases was not done, they were found to be fairly similar as would be expected if both cases approached a fully resolved VOF condition.
(a)
(b)
4. Conclusions
The flow phenomena observed for the impact and entrainment of a vertical plunging liquid jet offers a challenge to traditional regime-dependent CFD simulations and thus is a good test of the capabilities of the enhanced hybrid multiphase solver. It was observed that the results were sensitive to the details of the bubble diameter models and switching algorithms employed. When deactivation of interface sharpening was more substantial and predicted bubble size for the dispersed phase remained small (rPBE, normGrad), the simulation tended toward multifluid physics with a large dispersed plume and low entrainment. Conversely, when interface sharpening remained activated due to better mesh resolution or larger bubble size (e.g., critical We diameter model) the simulation result overall appeared more like a VOF model with most of the bubbles being resolved; however, even in this case the hybrid nature of the solver allowed for regions of high interface curvature (low interface resolution quality) to be unsharpened and revert to multifluid physics for modeling the transport of the bubbles on the subgrid scale according to the bubble diameter and interphase momentum transport terms (i.e., interphase drag) within the multifluid formulation. More work is needed to refine the transition between resolved and unresolved scales and determine the switching parameters which are appropriate for various classes of multiphase flows. At present no explicit transfer of fluid particle size information was made when transitioning between VOF-captured and multifluid modeled regions. Additionally, efforts at developing the methodology for prediction of experimentally observed fluid particle sizes require careful consideration of the underlying models used in the source term of the particle number density transport equation.
As has been noted elsewhere [10] and as evidenced by the large number of breakup and coalescence models available in the literature and the wide discrepancy among the corresponding model parameters reported in the literature, it is overly optimistic to consider any PBE-based model as predictive unless the underlying submodels for breakup and coalescence have been carefully calibrated to specific data for the fluid system and physics to be considered. This is a challenge yet to be conquered in the area of predictive population balance modeling. As these models also rely on accurate prediction of the local turbulence dissipation rate, significant research is also needed in the area of turbulence modeling and physically rigorous application of large-eddy simulation techniques to multiphase systems.
An alternative approach that can be facilitated by the hybrid multiphase simulation methodology presented here is to consider multiphase flows with the understanding that they are inherently multiscale in nature and a multiscale simulation approach is required. One such promising strategy has been proposed by Van den Akker [39] in which multiphase direct numerical simulations (DNS) are performed to predict droplet/bubble dynamics on the subgrid scale and these results are then coupled in some way to a macroscale simulation of the overall flow being explored. Given the disparity in computational cost at the macro- and microscales, rather than a direct coupling, information translation through an intermediary method such as in situ adaptive tabulation [40] offers one possible strategy. Lattice Boltzmann Methods (LBM) due to their excellent parallel scaling performance and high-order accuracy offer an attractive means for these multiphase DNS simulations.
One can envision a multiscale framework in which the hybrid multiphase solver employed here relies on simulation information for fluid particle dynamics (e.g., breakup/coalescence rate and size distribution) from dynamic LBM simulations being conducted for the specific flow conditions (e.g., fluid properties, phase fraction, and turbulence conditions) found in the concurrent macroscale simulation. Such a multiscale, multiregime, multiphase flow simulation strategy could have application to a wide variety of complex multiphase flows where current methods are unable to capture the diverse range of physics encountered. One example familiar to the authors and that has been the main target for the development of the methods demonstrated here is liquid–liquid extraction in which both dispersed and fully segregated multiphase flow regimes for three phases (typically two liquids and air) are found in a single flow device Wardle [41].
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors gratefully acknowledge the computation resources of the Fission Linux Cluster at Idaho National Laboratory. Visualization was done with ParaView 4.1 (http://www.paraview.org/). This work was supported by the U.S. Department of Energy, Office of Nuclear Energy, under Contract DE-AC02-06CH11357. Government License Notice: the submitted paper has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science Laboratory, is operated under Contract no. DE-AC02-06CH11357. The U.S. government retains for itself and others acting on its behalf a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce and prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the government.