Abstract
The flows in stage-wise liquid-liquid extraction devices include both phase segregated and dispersed flow regimes. As a additional layer of complexity, for extraction equipment such as the annular centrifugal contactor, free-surface flows also play a critical role in both the mixing and separation regions of the device and cannot be neglected. Traditionally, computional fluid dynamics (CFD) of multiphase systems is regime dependent—different methods are used for segregated and dispersed flows. A hybrid multiphase method based on the combination of an Eulerian multifluid solution framework (per-phase momentum equations) and sharp interface capturing using Volume of Fluid (VOF) on selected phase pairs has been developed using the open-source CFD toolkit OpenFOAM. Demonstration of the solver capability is presented through various examples relevant to liquid-liquid extraction device flows including three-phase, liquid-liquid-air simulations in which a sharp interface is maintained between each liquid and air, but dispersed phase modeling is used for the liquid-liquid interactions.
1. Introduction
While multiphase flows present unique challenges for computational fluid dynamics (CFD) simulation, a host of solution methods exist for simulation of well categorized flows. For “dispersed” flows in which one phase is continuous and the other is distributed in fine droplets, one can use Lagrangian particle tracking for small phase fractions (less than ~10%) in which each individual fluid particle is followed through the fluid in response to local flow conditions. For high phase fraction dispersed flows, a multi-fluid Eulerian-Eulerian solution method with interphase mass and momentum transfer can be applied. For stratified flows in which the fluid phases have a clearly defined phase interface, free-surface capturing methods such as Volume of Fluid (VOF) can be employed. Real flows, such as those encountered in liquid-liquid extraction devices, are not so easily categorized and can span multiple flow regimes (both spatially and temporally). In theory, interface capturing methods could be used for direct simulation of dispersed flows given that a mesh spacing of ~10x smaller than the smallest droplet can be used; however, accurate physical capturing of droplet-droplet interactions requires yet finer mesh resolution or droplet coalescence is severely overpredicted. In practice such meshing—and the small timesteps required (on the order of s) for stable solution—is not feasible for realistic turbulent multiphase flows and will not be in the foreseeable future even on large computers unless CFD algorithm developments are made which allow significant timestep acceleration or time parallelization (spatial decomposition is the only option currently for CFD parallelization). This timestep limit is due to the fact that interface capturing methods require explicit solution and are limited by Courant number: Thus, the timestep is directly proportional to the mesh spacing ( is the flow velocity)—that is, if the mesh spacing is cut in half, the timestep must essentially be decreased by the same margin. Consequently, for complex multiphase flows in which both dispersed flow and segregated flow regions are present one would like to couple these two methods into a single solver. In such a method, interface capturing would be used in regions where meshing is sufficient to resolve large droplets and bulk fluid-fluid interfaces or for phase pairs where interdispersion can be neglected; dispersed flow models would be used in regions where droplet characteristics move into the “subgrid” scale. As an example, for complex multiphase flows such as those found in liquid-liquid extraction devices where two immiscible liquids are mixed and air can also be present, one could employ sharp interface methods for certain phase pairs (e.g., liquid-air) and at the same time use dispersed modeling for others (e.g., liquid-liquid).
The idea of coupling these two methods for solution of such flows was explored by Cerne et al. [1]. They employed a simplified switching routine based on the gradient of the volume fraction across neighboring cells to flag cells as either VOF or two-fluid and solved the appropriate number of equations in each cell—resulting in complicated numerical issues due to solving models with different numbers of equations across the same domain [2]. To avoid such issues, this same research team has shifted toward multi-fluid-VOF coupling via the addition of interface capturing on top of an Eulerian multi-fluid solver. In this way, the multi-fluid formulation with momentum equations for each phase is applied across the entire domain and an interface sharpening algorithm (in their case a conservative level-set method which is similar in concept to the interface compression method described later) is applied for sharp interface regions [2, 3]. Strubelj and Tiselj [2] give a good overview of methods that have been employed for this coupling along with details regarding difficulties in coupling the phase momentum equations at the sharp interface (where the phase velocities should be equal). Again, the simple switching function of Cerne et al. [1] has been used by these authors who acknowledge its somewhat arbitrary nature and identify this as a key area of work to make the coupling more physically based.
A recent study by Yan and Che [4] has also attempted to couple multi-fluid and VOF methods with a “switching” mechanism based on the identification and transport of “large” (grid-resolved) and “small” (sub-grid) interface structures with models for exchange between the two. While this provides a means of phenomenological switching between the grid-resolved and subgrid scales and enables the use of accurate interface reconstruction methods (piece-wise linear interface construction (PLIC) was used), it requires the solution of two transport equations for each dispersable phase and is not easily generalizable to multiple phases beyond two.
The current study is part of a broader research effort to deliver computational tools for detailed simulation of liquid-liquid extraction (solvent extraction) unit operations with the aim of providing a pathway for prediction of key operational performance measures (e.g., stage efficiency, and other-phase carryover (backmixing)) for any given set of conditions. Such predictive capability will help inform process-level modeling tools and deliver insight into unit design and operation. To accomplish this, methods are required which can adequately predict liquid-liquid mixing and interfacial area generation as well as the formation and transport of small droplets. Liquid-liquid contacting equipment used in solvent extraction processes has the dual purpose of mixing and separating two immiscible fluids. Consequently, such devices inherently encompass a wide variety of multiphase flow regimes from segregated to dispersed flow types. From a simulation perspective one might perform separate simulations for the mixing [5] and separation [6] functions of stage-wise extraction devices such as centrifugal contactors or mixer settlers, but in this case assumptions must be made regarding the coupling between the two regions. This is problematic, especially in the case of the annular centrifugal contactor (Figure 1). In this device, mixing occurs in a narrow annulus between the stationary outer housing and the inner rotating cylinder (which itself is the centrifuge). The connection between the two regions occurs at the bottom of the rotor where the dispersion of the two fluids enters the separating region in the hollow rotor. These two regions are not necessarily independent and the coupling is poorly understood. These goals and simulation challenges have led us towards the development of a hybrid method building upon the work of the others mentioned above and extended for application to the multiphase flow and operation of these devices.
Of the equipment types generally used for solvent extraction processing of used nuclear fuel, centrifugal contactors have the largest relative knowledge gap and at the same time the greatest opportunity for significant benefits to a future fuel cycle facility. Thus, the present focus has been on these devices for which this multi-fluid-VOF coupling methodology is of particular importance to capture both the liquid-liquid dispersion flow as well as the complex, dynamic fluid-rotor interaction. A thorough review of CFD modeling efforts for annular centrifugal contactors (also called annular centrifugal extractors) has recently been published by Vedantam et al. [7]. Three-phase liquid-liquid-air simulation of the flow in an annular centrifugal contactor using a VOF-based method with single momentum equation for the mixture has been performed previously by Wardle [8]. While simulations of flow in centrifugal contactor-related geometries from other authors are limited, there is one available study in which a multi-fluid solution method is used to simulate the separation of two liquids in a simplified rotor geometry [9]. Recent work by Li et al. [10] reports three-phase simulations in a coupled mixer/separator centrifugal contactor model using a pure multi-fluid approach. Curiously, the authors found very little mixing between the two liquids. This appears to have been due to the assumption of a large dispersed phase diameter—however, it may also have been related to smearing of the fluid-rotor contact resulting from diffuse liquid-air interfaces. Coupling of the multi-fluid method with interface capturing as developed in this current effort allows simulation of annular centrifugal contactor mixing zone and rotor flows in which a liquid-air free surface is captured.
2. Computational Methodology and Implementation
The implementation of the methodology described here has been done using the OpenFOAM CFD package (version 2.1) which provides a collection of libraries and utilities for constructing custom CFD solvers for a wide variety of applications. A working version of the solver named has been included as part of the release of version 2.1 of the toolkit and is available for download at http://www.openfoam.org/.
3. Multifluid Momentum Equations and Implementation of Interphase Drag Models
The multi-fluid 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 are the inter-phase momentum transfer or drag force and the surface tension force . For the examples shown here, only the drag force is accounted for—though surface tension capability (based on the continuum surface force model of Brackbill [12]) and surface contact angle effects have also been subsequently added to the solver. The drag term is given by where the subscripts and denote the continuous and dispersed phase values and where is As implemented in the solver, the drag calculation is generic such that the model must simply return the value of . A variety of correlations could be used for calculation of the drag coefficient (in (5)) and several common models are available in OpenFOAM. As an example, the commonly used model of Schiller and Naumann [13] was used for the test cases here. In this model, the drag coefficient is a function of the Reynolds number according to where in which is the continuous phase kinematic viscosity. In this solver, calculation of the drag coefficient can be done by specifying a dispersed phase or by independent calculation with each phase as the “dispersed phase” and the overall drag coefficient applied to the momentum equations taken as the volume fraction weighted average of the two values. This is a so-called blended scheme which is a useful approximation for flows with regions in which either phase is the primary phase.
A constant droplet diameter size (independently defined for each phase) is used in the work reported here, but models for variable droplet size are compatible with this flexible framework. Indeed models based on a reduced population balance method have been developed and will be reported separately.
3.1. Interface Capturing
The interface sharpening methodology employed here is the interface compression method of Weller [14]. Gopala and van Wachem [15] give a useful comparison of this type of method with several other interface sharpening and reconstruction algorithms employed in a variety of CFD codes. The interface compression scheme of Weller [14] adds an additional “artificial” compression term to the LHS of the volume fraction transport equation (2) as where the velocity is applied normally to the interface to compress the volume fraction field and maintain a sharp interface. The term ensures the term is only active in the interface region. In addition, a bounded differencing scheme is employed for discretization of (8). The value for the artificial interface “compression velocity” is given by The term gives the interface unit normal vector for the direction of the applied compression velocity. The magnitude of the velocity is used since 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. While can mathematically be any value 0, if we restrict , (9) reduces to and is then 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 multi-fluid 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 implementation of the solver is configured such that the interface compression coefficient is 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 general, the interface compression method for interface capturing is not as physically accurate as interface reconstruction methods such as Piecewise Linear Interface Construction (PLIC). However, it is much easier to implement and performs faster, and most importantly, unlike PLIC it is mass conservative [15]. One problem that this method can suffer from is parasitic wavy currents at the interface which are particularly problematic for small-scale, surface tension driven flows (e.g., capillary rise), but become less important for advective flows such as is the current target of this work. The false interfacial currents can also be minimized by maintaining a low number through subtimestepping and restricting as was done here. For a discussion of other limitations common to VOF methods see [16].
3.1.1. Dynamic Switching
With the computational framework as outlined above, it is possible to imagine a unified method which allows simulation of complex flows with any combination of regimes ranging from fully dispersed to fully segregated. However, as noted by the various researchers who have attempted coupling between multi-fluid and VOF methods [1, 2, 4], the key area of uncertainty is the method by which switching occurs between the dispersed and segregated models. As outlined above, the current implementation could be made to allow for dynamic switching of the interface sharpening only in regions where the flow is segregated through implementation of a spatially nonuniform field(s). Others have suggested that one could choose to switch according to some predetermined flow regime map [2]. Such an approach, however, would likely be useful only for simple geometry flows (e.g., pipe flow) and a more general physics-based approach is needed for application to a broader class of flows such as those under investigation here. From a more general physical perspective, one would like to use the sharp interface method in regions of the flow. From a more general physical perspective, one would like to use the sharp interface method in regions of the flow where the actual droplet size and corresponding local mesh resolution is sufficient for accurate capture of droplet curvature. Conversely, where the droplet size falls below the mesh-resolvable scale, inter-phase dispersion modeling would be employed (sharpening deactivated). Such a scheme would require prediction of local droplet size based on a population balance or similar method and comparison with local mesh spacing using some cutoff criteria.
For the present implementation limited testing was done for a two-phase version of the solver using a switching function based on the work by Cerne et al. [1] and later used by others [3]. This technique applies a switch according to the magnitude of the gradient of the volume fraction itself—with the assumption being that when the gradient of the volume fraction becomes smaller than some cutoff value it is an indication of actual phase dispersion and interface sharpening is turned off. In this case, the normalized magnitude of the gradient of the volume fraction () is used as defined by Thus, when the value of is set to 1. A cutoff value of 0.4 is recommended by Cerne et al. [1]; however, the formulation used here has been normalized to 1 so the corresponding value may be somewhat different. One drawback of this mechanism for switching is that it is somewhat of a “self-fulfilling prophecy.” Where the interface is already sharp (steep gradient in the volume fraction ), you apply interface compression to keep it sharp, and where it is dispersed, you let it stay dispersed. Even so, it gives a reasonable preliminary model for the coupling and has been implemented and tested for a two-phase version of the solver only. An example case is described in Section 4.1.
For the majority of examples reported here, which are three-phase liquid-liquid-air cases, it was assumed that dispersion of air into the liquids can be neglected and thus the parameter was fixed for each phase pair. A value of 1 (interface capturing on) was used for any air-liquid interfaces and interface capturing was turned off ( = 0) for the liquid-liquid interfaces. In this case, no entrained air will be captured and the air-liquid interface will be everywhere sharp while the liquid phases will be allowed to interdisperse and no sharp interfaces will be imposed. This assumes that the entrainment of air has a negligible effect on the flow and mixing of the two liquids.
3.2. Solution Procedure for Multifluid-VOF Coupling
The general solution procedure for the hybrid solver using the equations above is as follows: (1)update timestep according to Courant number limit (ratio of timestep to interface transit time in cell);(2)solve coupled set of volume fraction equations with interface sharpening for selected phase pairs ((8) with multiple subtimesteps); (3)compute drag coefficients;(4)construct equation set for phase velocities and solve for preliminary values; (5)solve pressure-velocity coupling according to Pressure Implicit with Splitting of Operators (PISO) algorithm:(a)compute mass fluxes at cell faces;(b)define and solve pressure equation (repeat multiple times for non-orthogonal mesh corrector steps); (c)correct fluxes; (d)correct velocities and apply BCs; (e)repeat for number of PISO corrector steps; (6)compute turbulence and correct velocities; (7)repeat from 1 for next timestep.
3.2.1. Numerical Considerations for Stability of Momentum Coupling and Phase Conservation
In the limit of a sharp interface, the velocities on either side of the interface must be equal to meet the so-called no-slip interface condition. This is an inherent feature of traditional VOF simulations as all phases share a single momentum equation and thus the phase velocities are the same everywhere. Imposition of a sharp interface through the addition of interface compression “on top” of a multi-fluid formulation in which each phase has its own momentum equation requires that an additional artificial drag is imposed to equalize the velocities at the interface. In the work of Strubelj and Tiselj [2] in which multi-fluid-VOF coupling was performed, an arbitrary function proportional to the inverse of the time step divided by 100 (resulting in a large value) was imposed to force large interphase drag coefficients at the interface. In this case, rather than devise some arbitrary formulation, small “residual drag” and “residual phase fraction” constants were added for each phase pair (typically both equal to ) to stabilize the phase momentum coupling. These added residual values were only used in calculation of the drag for momentum coupling stability and did not affect actual phase fractions or overall phase conservation.
In order to ensure phase conservation for the coupled phase fractions with added interface sharpening, it was necessary to incorporate limiters on the phase fraction as well as on the sum of the phase fractions prior to the explicit solution of the phase fraction equation system. These additional limiters have been incorporated in a new multiphase implementation of the Multidimensional Universal Limiter with Explicit Solution (MULES) solver framework within OpenFOAM. This multiphase MULES implementation used in is also leveraged for enhancing the phase conservation performance of the n-phase VOF-only solver . As in the case of the standard VOF solver, solution of the volume fraction transport equations was done using subtimestepping over several subintervals of the overall time step to maintain solution stability according to the Courant number limit (1) while maximizing overall timestep to minimize time to solution for the transient solver. It was found that an overall Cr number limit of 1.5 (based on velocities near sharp interface) with 5 subtimesteps could deliver stable results.
4. Results: Example Cases
The following example cases demonstrate the capability of the simulation methods to capture, on a per-phase-pair basis, both dispersed and segregated flows. Only the first case, on the breaking of a dam considers treatment of the sharpening coefficient as a volumetric field with local, dynamic switching based on the gradient of the volume fraction as mentioned above. In all other cases, the value of is generally set to 1 (imposed interface sharpening) for liquid-air interfaces and to 0 for liquid-liquid phase interactions. The properties used are representative of water, oil, and air at room temperature conditions. As done in the work of Padial-Collins et al. [9], a constant droplet size of 150 microns was assumed for the liquids. A value of 1 mm was used for air. Interphase drag was treated via the blended method mentioned earlier. Visualizations were done using ParaView version 3.12.
4.1. Liquid-Liquid “Column” Collapse
This example is a modification to the classic collapsing liquid column 2D test case in this case for a liquid-liquid system with two initial regions of dispersed phase volume fraction as shown in Figure 2. The domain is 58.4 cm square, with the short barrier on the bottom surface having a width of 2.4 cm and a height of 4.8 cm. Figure 3 shows a comparison of successive time snapshots of simulations having set the interface sharpening coefficient to 0 (VOF behavior) and 1 (multi-fluid behavior). The behavior is as expected for the two cases—that is, with , immediately upon startup droplets with a characteristic size similar to the mesh size are formed, and despite this, overall segregation of the phases is relatively slow. When interface sharpening is not imposed () and the multi-fluid behavior is governed by the interphase drag correlation, separation of the phases appears to be more physical and occurs on a faster time scale.
(a)
(b)
A variation of the above case was performed starting from the same initial state but in which the value of was allowed to vary locally (as 0 or 1 only) according to (11). Successive time snapshots of phase fraction and the sharpening coefficient field are shown in Figure 4. Note that the interphase drag has been increased in this case (through a larger assumed droplet size) resulting in slightly slower separation dynamics as compared to the earlier case shown in Figure 3.
(a)
(b)
It was observed that around s the first region of active interface sharpening appears, but left-right interface motion is sufficient that the sharpened region is not maintained. Between 10 and 12 s, a stable, sharpened interface appears and grows until it covers the length of the phase interface around 16-17 s.
This example demonstrates the functionality of dynamic interface sharpening switching based on the volume fraction gradient for a simple test case. As noted earlier, a switching function of this type is somewhat arbitrary. Ideally, one would like a physical basis for governing switching according to a comparison of the local predicted droplet size and the local mesh spacing—where mesh resolution is sufficient to resolve droplets adequately, sharpening is activated, and where droplet size falls into the subgrid scale, sharpening is deactivated. One could imagine a very flexible model which could simulate multiple flow regimes in this way in a multi-scale manner with the multifluid method capturing sub-grid phase particle transport analogous in principle to the idea of the sub-grid scale modeling done in Large Eddy Simulation (LES).
4.2. Horizontal Settler
The next test case simulates the separation of a 50 : 50 liquid-liquid dispersion entering into a 2D, horizontal settler as shown in Figure 5 where the initial state is stratified layers of oil (red) and water (blue). Gravity acts in the vertical direction. Flows of each phase exit from the corresponding surfaces in the upper right of the domain. The length of the main body of the domain is 10 cm and the overall height is 2.25 cm. The development of a “dispersion band” between the regions of separated phases was observed. The dispersion band was not static, but was found to be disturbed by longitudinal waves generated by a periodic vortex at the back of the so-called dispersion disk (wall just upstream of the inlet). The dispersion disk is placed just upstream of the inlet to direct the entering dispersion toward the central vertical position of the separating region and try to minimize downstream disturbances. Simulations were done with no dispersion disk in order to verify the effectiveness of this feature. It was indeed observed that the overall width of the dispersion band is greater without the dispersion disk as inlet disturbances propagate farther downstream and disrupt the separation of the two phases leading to more entrainment in the exit streams.
(a)
(b)
(c)
(d)
(e)
4.3. Annular Mixer
As noted above, a principle application requiring the capability of this hybrid solution method is the liquid-liquid-air flow in an annular centrifugal contactor. This device (Figure 1) consists of an annular region with a rotating inner cylinder and stationary outer cylinder in which the two immiscible liquids are mixed in the presence of a free surface. This complicates the physics significantly, requiring sharp interface capturing to accurately predict the intermittent liquid-rotor contact [17]. At the same time, dispersed phase modeling is needed to predict the mixing and flow of the liquid-liquid dispersion.
To demonstrate the capability of the solver to capture such flow dynamics, simulations were conducted in an idealized annular mixer both for a 2D, axisymmetric case and a fully 3D case. In both cases, the inner radius is 2.54 cm and the outer 3.17 cm (annular gap of 0.63 cm) and the height of the annulus is 7 cm. The top surface is open to air at constant atmospheric pressure and the bottom surface is treated as a wall. Unless otherwise stated, the rotation rate of the inner cylinder is 3600 RPM (377 rad/s) resulting in a surface velocity of 9.56 m/s. Turbulence was treated using Large Eddy Simulation (LES) with the Smagorinsky sub-grid model. A uniform quadrilateral mesh was used for the 2D model with spacing of 0.2 mm (32 cells across the annular gap). In order to explore mesh dependency of the new solver, additional simulations were done with mesh spacings of 0.4 mm (coarse) and 0.1 mm (fine). The relative mesh spacings for the three sizes can be seen in Figure 6. For the 3D model, the base case simulation was done with a mesh spacing ~0.4 mm (15 hexahedral cells across the annular gap, 675 K cells total) though for comparison an additional run was also done for a finer mesh (~0.25 mm, 2.4 M cells).
(a)
(b)
(c)
4.3.1. , Axisymmetric Model
Figure 7 gives a time series of snapshots showing volume fractions for water (blue), oil (red), and air (cyan) from startup through s for simulation on the medium mesh refinement (Figure 6(M)). It is clear that even for this very turbulent flow, the hybrid solver is able to maintain a sharp interface for the liquid-air free surface while at the same time allowing phase inter-dispersion for the two liquids. Unlike the simplifications often used by other CFD studies of this type of flow (e.g., [18]), the presence of air and the existence of the free surface has a significant impact on the characteristics of the flow and breaks down any Taylor-Couette vortices that would be present in the liquid-liquid only case. It was observed that there was one relatively stable vortex at the bottom which was characterized by a light-phase rich region at the center and rotation in the clockwise direction (flow inward along the bottom surface). The companion vortex (counterclockwise rotation) just above this lower one was found to be periodically formed and then break away and travel upward as the liquids are spun off the rotor and move toward a maximum height on the outer wall.
Sharp interface capturing methods such as the VOF method used here are inherently mesh dependent; the finer the mesh, the finer the interface features that can be captured. In order to explore the mesh dependency that the hybrid solver inherits from the VOF method, simulations for the 2D axisymmetric annular mixer model were performed on two additional meshesmone twice as coarse and one twice as fine as the base case. Snapshots of the phase fractions at s after startup (a) and for a time average over the period from s to s (b) are shown in Figure 8. Even for the relatively short time after start-up shown in Figure 8(a) the flows observed for the simulations on the three meshes have clear differences. Additionally, the differences are apparent not only in the shape of the liquid-air free surface, which is to be expected, but also in the multi-fluid dispersion behavior between the two liquids. Comparison of the time-averaged behavior (Figure 8(b)), however, shows better general consistency in the overall liquid height. The Taylor-Couette vortex near the bottom mentioned previously was found to be most prominent in the medium mesh case.
(a)
(b)
Though it is not readily apparent from the snapshots in Figure 7, as has been observed previously for flow in this configuration at such conditions [17], the overall height of the liquid on the outer wall exhibited an oscillatory behavior. Figure 9 shows a plot of the liquid height on the outer wall versus time for the 2D axisymmetric simulation for the three different meshes. As noted previously, there is significant variation in the temporal evolution of the flow on the three mesh densities while the general behavior is similar. As will be shown later, the oscillations in liquid height for the 3D simulations exhibited a much more clear periodicity; for the 2D simulations, the oscillation magnitude in liquid height was largest for the most coarse mesh. The minimum height of the oscillation of the liquid corresponds with a maximum contact area between the liquid and the rotor after which the liquid is accelerated and spun out and up the housing wall leading to a maximum liquid height corresponding to a minimum in fluid-rotor contact. The amount of overall contact between the liquid and the rotor has an impact on the level of mixing that occurs between the two liquid phases. Figure 10 shows a plot of the time-averaged liquid contact on the rotor side and integrated values corresponding to the fractional liquid “coverage” in each case. While not unexpected, there are complex mesh dependencies for the hybrid solver which require additional investigation.
In order to provide a point of reference for the additional computational cost of the multi-fluid hybrid scheme versus an all-VOF simulation (single, shared momentum equation and sharp interfaces everywhere), a simulation was done using OpenFOAM’s multiphase-capable VOF solver for the fine mesh case (44,800 hexahedral cells) of the 2D annular mixer problem described here. The simulation was done using the exact same solver settings (discretization schemes, number of subtimesteps on volume fraction solutions, etc.) and the same number of parallel processors (12 cores were used in this case). The simulations were compared out to 1 second of flow and it was found that the hybrid multi-fluid/VOF solver is only 39% more costly per timestep than the comparable case with the all-VOF solver (8.69 CPU seconds/step versus 6.23 CPU seconds/step for all-VOF solver). As the bulk of the computational time is spent in the solution of the volume fractions and more so in the pressure-velocity coupling (Items #2 and #5, resp., in the solution procedure given in Section 3.2), only a modest computational cost is incurred due to the additional phase momentum equations and momentum coupling in the multi-fluid formulation that is not required in the VOF-only solver.
4.3.2. Model
Figure 11 shows snapshots of the liquid phase fractions in the 3D annular mixer model soon after start-up (~0.25 s) and ~3 s after startup at 3600 RPM. In contrast with the 2D axisymmetric approximation, there is significant azimuthal variation in the flow and liquid-rotor contact. The liquid-liquid dispersion also exhibited a steady height oscillation due to being periodically thrown off the rotor and up the outer wall. A plot of the liquid height on the outer wall (azimuthal average) for two rotor speeds (2400 and 3600 RPM) are shown in Figure 12 along with the 3600 RPM value for the fine mesh for comparison. At the lower rotor speed, the liquid exhibits only minimal height variation with some periodicity while at the high rotor speed a steady oscillation develops ~1 s after startup. Table 1 gives a comparison of the liquid height oscillation frequency (calculated from trough-trough times for the ~6 cycles between 1.5–3 s) for the 3D model with different variations: the base case blended drag model with the droplet diameters for both liquid phases at 150 microns, aqueous dispersed with 150 micron droplet size, and aqueous dispersed with a 50 micron droplet size. The results for the fine mesh (blended case only) are also shown for comparison. There is not a very strong dependence on oscillation frequency for blended versus aqueous dispersed at the same droplet diameter. For the smaller droplet diameter, the oscillation frequency was found to decrease slightly. In terms of mesh dependency, despite oscillation “phase” offset due to differences at early times, the fully developed oscillation frequency is very comparable for the two mesh densities.
A companion experimental effort has also been initiated to provide means of validating the advanced simulation capability presented by this new solver for flows in actual annular centrifugal contactor equipment. While the simulations presented here are for a simple annular mixer, it was found that certain characteristics of the flow showed good agreement with preliminary experimental observations and a brief comparison will be made here. Experiments were conducted using a CINC-V2 centrifugal contactor modified with a quartz outer cylinder as reported by Wardle et al. [17]. The tests included here were done using a liquid-liquid system consisting of 1 M nitric acid with 1 M aluminum nitrate as the aqueous phase () and 40% (by volume) tributyl phosphate in dodecane (). The aqueous phase was dyed with methylene blue to aid in visual phase discrimination. While tests were done for a variety of housing vane types and inlet flow rates, conditions were selected for comparison here which gave a comparable annular liquid height in the mixing zone. For additional details regarding the experiments, please see Wardle 2012 [11].
It was observed that while the 3D annular mixer geometry has a fixed volume with no inlet/outlet flow, the general features of the flow compared favorably with observations in experiment. Figure 13 shows a comparison snapshot of the liquid-liquid flow as viewed from the housing side from the simulation (a) and a high-speed image from experiment (b). Both images are near a minimum liquid height and show similar banding and tendrils of the dispersed, heavy phase (blue). In both simulation and experiment, it was observed that as the dispersion approached a minimum liquid height, aqueous phase striations such as these appeared. Additionally, contact with the rotor for the collapsing liquid resulted in a sharp “spurting” of air across a line near the top of the main body of liquid followed by the rise of the liquid up the outer wall. In addition to such qualitative comparisons, the oscillation frequency of the liquid surface was also found to be quite comparable quantitatively under conditions of similar liquid height. As reported previously in Table 1, the mean oscillation frequency from simulation was found to be 4.6 Hz. From high-speed video analysis, the frequency observed in experiments was slightly higher at 4.9–5.3 Hz. Values from two different cases (4 straight mixing vanes (4V) and 8 curved mixing vanes (CV)) with similar liquid height are given in Table 1. Note that this is also quite comparable to what was found previously for similar conditions under single-phase flow [19] which ranged from 4.5 to 5.4 Hz at 3600 RPM depending on the total feed flow rate—with high frequency being observed at lower flow rates (lower overall liquid height). This seems to demonstrate that the phenomenon of liquid height oscillation (and the corresponding periodic liquid–rotor contact) is more a function of the annular geometry and rotor speed than it is of the fluid properties, flow configuration (whether net axial flow is present or not), or geometry beneath the rotor—provided that the conditions and configuration result in comparable liquid volume (height) in the annular region.
(a)
(b)
5. Conclusions
A hybrid multiphase CFD solver has been developed which combines the Euler-Euler multi-fluid methodology with VOF-type sharp interface capturing on selected phase-pair interfaces. A variety of examples of cases have been presented here which are relevant to liquid-liquid extraction and demonstrate the functionality and flexibility of the new solver. The multiphase flow simulation capability described here has application to a variety of complex flows which span multiple regimes from fully segregated to fully dispersed. While the target application is flow in liquid-liquid extraction devices, this methodology could be useful to the simulation of other multiphase flows which are currently restricted to a single flow regime. For example, coupling these methods with heat transfer and phase change could provide a tool for simulations of gas-liquid channel flows such as those seen in nuclear reactors spanning bubbly, slug, churn-turbulent, and annular flows.
The primary goal of the overall research effort of which this work is a part is the prediction of mass-transfer efficiency in stage-wise liquid-liquid extraction devices. This requires prediction of liquid-liquid interfacial area and consequently capturing of the dispersed phase droplet size distribution in some manner. The solver presented here provides a flexible foundation for building in the necessary models from any of the available methods. In the solver as part of the OpenFOAM release, the droplet diameter of each phase is a field variable and has been implemented as a callable library such that additional droplet diameter models can be easily implemented and selected at runtime. Extension of the solver to include variable droplet size using the reduced population balance method of Attarakih et al. [20] as implemented in [21] has recently been performed and exploration of droplet breakup and coalescence models is underway with very promising results.
Disclosure
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 US Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in the said paper to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the government.
Acknowledgments
The authors gratefully acknowledge the use of the Fusion Linux cluster at Argonne National Laboratory and the Fission Linux cluster at Idaho National Laboratory for computational resources. K. E. Wardle was supported by the US DOE Office of Nuclear Energy’s Fuel Cycle R&D Program in the area of Separations. K. E. Wardle would also like to thank Jing Gao of the University of Illinois at Chicago for experimental contributions as a summer research aide.