International Journal of Chemical Engineering

Volume 2015 (2015), Article ID 925639, 14 pages

http://dx.doi.org/10.1155/2015/925639

## Numerical Investigation of Vertical Plunging Jet Using a Hybrid Multifluid–VOF Multiphase CFD Solver

^{1}Department of Mathematical Sciences, Michigan Technological University, 1400 Townsend Drive, Houghton, MI 49931-1295, USA^{2}Chemical Sciences and Engineering Division, Argonne National Laboratory,
Argonne, IL 60439, USA

Received 9 April 2015; Revised 21 June 2015; Accepted 28 June 2015

Academic Editor: Jerzy Bałdyga

Copyright © 2015 Olabanji Y. Shonibare and Kent E. Wardle. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.