Simulation Of Multi-fluid Flows Research Articles

A Consistent and Implicit Rhie–Chow Interpolation for Drag Forces in Coupled Multiphase Solvers

The use of coupled algorithms for single fluid flow simulation has proven its superiority as opposed to segregated algorithms, especially in terms of robustness and performance. In this paper, the coupled approach is extended for the simulation of multi-fluid flows, using a collocated and pressure-based finite volume discretization technique with a Eulerian–Eulerian model. In this context a key ingredient in this method is extending the Rhie–Chow interpolation technique to account for the unique flow coupling that arises from inter-phase drag. The treatment of this inter-fluid coupling and the fashion in which it interacts with the velocity-pressure solution algorithm is presented in detail and its effect on robustness and accuracy is demonstrated using 2D dilute gas–solid flow test case. The results achieved with this technique show substantial improvement in accuracy and performance when compared to a leading commercial code for a transonic nozzle configuration.

A five-equation model based global ale method for compressible multifluid and multiphase flows

In this work, a mixing cell closure model was derived in a Lagrangian formulation based on the hypothesis of isentropic compression or expansion. The proposed closure model is used to construct a five-equation model for the simulation of multifluid and multiphase flows. Thereafter a global ALE(Arbitrary Lagrangian-Eulerian) method was developed for the five-equation model. The five-equation based ALE method is a two-stage ALE method, including a Lagrangian phase and a rezone-remap phase. In the first phase, the model equation was discretized with a cell-centered Lagrangian scheme, and a multifluid Riemann solver was proposed by extending traditional Riemann solver to mixed cells with different materials. The five-equation based ALE method can be used to simulate multi-material flows with large deformation, which is a challenge for many traditional ALE methods. Moreover, it can be used to simulate multiphase flows. A series of multi-material and multiphase test problems were simulated with the five-equation based ALE method, and numerical results agree well with exact solutions and reference results.

Numerical simulation of compressible multifluid flows using an adaptive positivity‐preserving RKDG‐GFM approach

SummaryThe Runge‐Kutta discontinuous Galerkin method together with a refined real‐ghost fluid method is incorporated into an adaptive mesh refinement environment for solving compressible multifluid flows, where the level set method is used to capture the moving material interface. To ensure that the Riemann problem is exactly along the normal direction of the material interface, a simple and efficient modification is introduced into the original real‐ghost fluid method for constructing the interfacial Riemann problem, and the initial conditions of the Riemann problem are obtained directly from the solution polynomials of the discontinuous Galerkin finite element space. In addition, a positivity‐preserving limiter is introduced into the Runge‐Kutta discontinuous Galerkin method to suppress the failure of preserving positivity of density or pressure for the problems involving strong shock wave or shock interaction with material interface. For interfacial cells in adaptive mesh refinement, the data transfer between different grid levels is achieved by using a L2 projection approach along with the least squares fitting. Various numerical cases, including multifluid shock tubes, underwater explosions, and shock‐induced collapse of a underwater air bubble, are computed to assess the capability of the present adaptive positivity‐preserving RKDG‐GFM approach, and the simulated results show that the present approach is quite robust and can provide relatively reasonable results across a wide variety of flow regimes, even for problems involving strong shock wave or shock wave impacting high acoustic impedance mismatch material interface.

The simulation of compressible multi-fluid flows by a GRP-based energy-splitting method

This paper serves to simulate compressible multi-fluid flows using an energy-splitting method based on the generalized Riemann problem (GRP) solver. The emphasis is put on two-dimensional problems to address the fundamental role of transversal effect in the construction of numerical fluxes, which exhibits the genuine multi-dimensionality of the underlying numerical schemes. In contrast with the popularly-used methods of line which use exact or approximate Riemann solvers, the GRP solver, as the discontinuous version of the Lax–Wendroff method, takes the temporal-spatial coupling process to include the transversal effect. Combined with the energy-splitting technique developed earlier in Lei and Li [24], the resolution of flow fields around material interfaces is improved significantly.

Contact preserving Riemann solver for incompressible two-phase flows

In numerical simulation of multi-fluid flows, it is very crucial to accurately capture the interface separating the fluids, which has the properties of a contact wave. Since HLLC (Harten Lax van-Leer with contact) formulation has the ability to model contact waves exactly, a novel HLLC-type Riemann solver is proposed for calculation of convective flux in a system where the incompressible-Navier–Stokes equations are tightly coupled with the volume of fluid (VOF) equation. The efficacy of the new formulation is evaluated by solving various standard test problems involving two-fluid flows.

High-Order Finite-Element Framework for the Efficient Simulation of Multifluid Flows

In this paper, we present a comprehensive framework for the simulation of Multifluid flows based on the implicit level-set representation of interfaces and on an efficient solving strategy of the Navier-Stokes equations. The mathematical framework relies on a modular coupling approach between the level-set advection and the fluid equations. The space discretization is performed with possibly high-order stable finite elements while the time discretization features implicit Backward Differentation Formulae of arbitrary order. This framework has been implemented within the Feel++ library, and features seamless distributed parallelism with fast assembly procedures for the algebraic systems and efficient preconditioning strategies for their resolution. We also present simulation results for a three-dimensional Multifluid benchmark, and highlight the importance of using high-order finite elements for the level-set discretization for problems involving the geometry of the interfaces, such as the curvature or its derivatives.

A locally extended finite element method for the simulation of multi-fluid flows using the Particle Level Set method

Abstract The simulation of immiscible two-phase flows on Eulerian meshes requires the use of special techniques to guarantee a sharp definition of the evolving fluid interface. This work describes the combination of two distinct technologies with the goal of improving the accuracy of the target simulations. First of all, a spatial enrichment is employed to improve the approximation properties of the Eulerian mesh. This is done by injecting into the solution space new features to make it able to correctly resolve the solution in the vicinity of the moving interface. Then, the Lagrangian Particle Level Set (PLS) method is employed to keep trace of the evolving solution and to improve the mass conservation properties of the resulting method. While the local enrichment can be understood in the general context of the XFEM, we employ an element-local variant, which allows preserving the matrix graph, and hence highly improving the computational efficiency.

Parallel load balancing strategy for Volume-of-Fluid methods on 3-D unstructured meshes

Volume-of-Fluid (VOF) is one of the methods of choice to reproduce the interface motion in the simulation of multi-fluid flows. One of its main strengths is its accuracy in capturing sharp interface geometries, although requiring for it a number of geometric calculations. Under these circumstances, achieving parallel performance on current supercomputers is a must. The main obstacle for the parallelization is that the computing costs are concentrated only in the discrete elements that lie on the interface between fluids. Consequently, if the interface is not homogeneously distributed throughout the domain, standard domain decomposition (DD) strategies lead to imbalanced workload distributions. In this paper, we present a new parallelization strategy for general unstructured VOF solvers, based on a dynamic load balancing process complementary to the underlying DD. Its parallel efficiency has been analyzed and compared to the DD one using up to 1024 CPU-cores on an Intel SandyBridge based supercomputer. The results obtained on the solution of several artificially generated test cases show a speedup of up to ∼12× with respect to the standard DD, depending on the interface size, the initial distribution and the number of parallel processes engaged. Moreover, the new parallelization strategy presented is of general purpose, therefore, it could be used to parallelize any VOF solver without requiring changes on the coupled flow solver. Finally, note that although designed for the VOF method, our approach could be easily adapted to other interface-capturing methods, such as the Level-Set, which may present similar workload imbalances.

Parallelization Strategy for the Volume-of-fluid Method on Unstructured Meshes

The Volume-of-Fluid (VOF) is one of the most widely used methods for interface tracking in the simulation of multi-fluid flows. The interface between different fluids is generated from the volume fraction scalar fields, which account for the ratio of volume of each fluid in each control volume. Then, an advection equation is solved to obtain the new distribution of the fluids after momentum is applied. Since this is a time-consuming process, parallelization techniques play an essential role. In the VOF approaches most of computing cost of the algorithm is concentrated in operations with the cells that form the interface, i.e. the cells in which coexist different fluids. When the interface is not homogeneously distributed throughout the domain, the standard domain decomposition strategy results in an unbalanced partition. A possible strategy to overcome this limitation is to adapt the domain decomposition to the interface distribution, however, this approach presents a number of drawbacks mainly related to the dynamic location of the interface. In this paper a new strategy, based in a load balancing process complementary to the domain decomposition, is presented with the aim to overcome the limitations of standard domain decomposition based approaches

Sensor and Numerical Simulator Evaluation for Porous Medium Desiccation and Rewetting at the Intermediate Laboratory Scale

Soil desiccation, in conjunction with surface infiltration control, is considered at the Hanford Site as a potential technology to limit the flux of technetium and other contaminants in the vadose zone to the groundwater. An intermediate‐scale experiment was conducted to test the response of a series of instruments to desiccation and subsequent rewetting of porous media. The instruments include thermistors, thermocouple psychrometers, dual‐probe heat pulse sensors, heat dissipation units, and humidity probes. The experiment was simulated with the multifluid flow simulator STOMP, using independently obtained hydraulic and thermal porous medium properties. All instrument types used for this experiment were able to indicate when the desiccation front passed a certain location. In most cases the changes were sharp, indicating rapid changes in moisture content, water potential, or humidity. However, a response to the changing conditions was recorded only when the drying front was very close to a sensor. Of the tested instruments, only the heat dissipation unit and humidity probes were able to detect rewetting. The numerical simulation results reasonably match the experimental data, indicating that the simulator captures the pertinent gas flow and transport processes related to desiccation and rewetting and may be useful in the design and analysis of field tests.

A simple incompressible flux splitting for sharp free surface capturing

SUMMARY This paper first applies a flux vector-type splitting method based on the numerical speed of sound for computing incompressible single and multifluid flows. Here, a preconditioning matrix based on Chorin's artificial compressibility concept is used to modify the incompressible multifluid Navier–Stokes equations to be hyperbolic and density or volume fraction-independent. The current approach can reduce eigenvalues disparity induced from density or volume fraction ratios and enhance numerical stability. Also, a simple convection-pressure flux-splitting method with high-order essentially nonoscillatory-type primitive variable extrapolations coupled with monotone upstream-centered schemes for conservation laws-type volume fraction recompressed reconstruction is used to maintain the preservation of sharp interface evolutions in multifluid flow simulations. Benchmark tests including a solid rotation test of a notched two-dimensional cylinder, the evolution of spiral and rotational shapes of deformable circles, a dam breaking problem, and the Rayleigh–Taylor instability were chosen to validate the current incompressible multifluid methodology. An incompressible driven cavity was also chosen to check the robustness of the proposed method on the computation of single fluid incompressible flow problems. Copyright © 2011 John Wiley & Sons, Ltd.

Desiccation of Unsaturated Porous Media: Intermediate‐Scale Experiments and Numerical Simulation

Soil desiccation (drying), involving water evaporation induced by air injection and extraction, is a potentially robust vadose zone remediation process to limit migration of inorganic or radionuclide contaminants through the vadose zone. Desiccation also has the potential to improve gas‐phase‐based treatments by reducing water saturation and therefore increasing sediment gas‐phase permeability. Before this technology can be deployed in the field, concerns related to energy limitations, osmotic effects, and potential contaminant remobilization after rewetting must be addressed. A series of detailed, intermediate‐scale laboratory experiments, using unsaturated homogeneous and heterogeneous systems, was conducted to improve our understanding of energy balance issues related to soil desiccation. The experiments were subsequently simulated with the multifluid flow simulator STOMP, using independently obtained hydraulic and thermal porous medium properties. In all experiments, the injection of dry air proved to be an effective means for removing essentially all moisture from the test media. Observed evaporative cooling generally decreased with increasing distance from the gas inlet chamber. The fine‐grained sand embedded in the medium‐grained sand of the heterogeneous system showed two local temperature minima associated with the cooling. The first one occurred because of evaporation in the adjacent medium‐grained sand, whereas the second minimum was attributed to evaporative cooling in the fine‐grained sand itself. Results of the laboratory tests were simulated accurately only if the thermal properties of the flow cell walls and insulation material were taken into account, indicating that the appropriate physics were incorporated into the simulator.

Direct numerical simulation of complex multi-fluid flows using a combined front tracking and immersed boundary method

In this paper a simulation model is presented for the Direct Numerical Simulation (DNS) of complex multi-fluid flows in which simultaneously (moving) deformable (drops or bubbles) and non-deformable (moving) elements (particles) are present, possibly with the additional presence of free surfaces. Our model combines a Front Tracking (FT) model developed by van Sint Annaland et al. (2008. Numerical simulation of dense gas–solid fluidized beds: a multiscale modeling strategy. Ann. Rev. Fluid Mech. 40, 47–70.) and an Immersed Boundary (IB) model developed by van der Hoef et al. (2008. Numerical simulation of dense gas–solid fluidized beds: a multiscale modeling strategy. Ann. Rev. Fluid Mech. 40, 47–70.) The FT part circumvents the explicit computation of the interface curvature. The IB part incorporates both particle–fluid and particle–particle interaction via a direct forcing method and a hard sphere Discrete Particle (DP) approach. In our model a fixed (Eulerian) grid is utilised to solve the Navier–Stokes equations for the entire computational domain. The no-slip condition at the surface of the moving particles is enforced via a momentum source term that only acts in the vicinity of the particle surface. For the enforcement of the no-slip condition Lagrangian force points are used, which are distributed evenly over the surface of the particle. Dissipative particle–particle and/or particle–wall collisions are accounted via a hard sphere DP approach using a three-parameter particle–particle interaction model accounting for normal and tangential restitution and tangential friction. The capabilities of the hybrid FT-IB model are demonstrated with a number of examples in which complex topological changes in the interface are encountered.

Scalable Modeling of Carbon Tetrachloride Migration at the Hanford Site Using the STOMP Simulator

Numerical simulation has been applied in support of the U.S. Department of Energy's efforts to characterize the nature and distribution of CCl4 contamination in the deep vadose zone at the Hanford site, near Richland, WA. Three‐dimensional simulations were executed using layered and heterogeneous distributions of soil properties to investigate the vertical and lateral distribution of CCl4 beneath its release point (216‐Z‐9 trench) and the effects of soil vapor extraction. The complexity of the modeled physical processes—namely, the nonlinearities associated with multifluid subsurface flow, including phase transitions and hysteresis in the relative permeability–saturation–capillary pressure functions—limits the grid resolution when executed using single‐processor computers. To achieve higher grid resolutions and acceptable detail in the simulation results for subsurface distribution and remediation of CCl4, execution on multiple processors was required. We have developed a scalable implementation of a multifluid subsurface flow and transport simulator, STOMP, with capabilities for volatile organic compounds. Developing scientific software for execution on parallel computers has unique challenges. The guiding objectives for developing this scalable code were to keep the source coding readable and modifiable by subsurface scientists, to allow for both sequential and scalable processing, and to depend on domain scientists for code parallelization and scalable linear system solvers.

An object-oriented and quadrilateral-mesh based solution adaptive algorithm for compressible multi-fluid flows

In this paper, an object-oriented and quadrilateral-mesh based solution adaptive algorithm for the simulation of compressible multi-fluid flows is presented. The HLLC scheme (Harten, Lax and van Leer approximate Riemann solver with the Contact wave restored) is extended to adaptively solve the compressible multi-fluid flows under complex geometry on unstructured mesh. It is also extended to the second-order of accuracy by using MUSCL extrapolation. The node, edge and cell are arranged in such an object-oriented manner that each of them inherits from a basic object. A home-made double link list is designed to manage these objects so that the inserting of new objects and removing of the existing objects (nodes, edges and cells) are independent of the number of objects and only of the complexity of O ( 1 ). In addition, the cells with different levels are further stored in different lists. This avoids the recursive calculation of solution of mother (non-leaf) cells. Thus, high efficiency is obtained due to these features. Besides, as compared to other cell-edge adaptive methods, the separation of nodes would reduce the memory requirement of redundant nodes, especially in the cases where the level number is large or the space dimension is three. Five two-dimensional examples are used to examine its performance. These examples include vortex evolution problem, interface only problem under structured mesh and unstructured mesh, bubble explosion under the water, bubble-shock interaction, and shock-interface interaction inside the cylindrical vessel. Numerical results indicate that there is no oscillation of pressure and velocity across the interface and it is feasible to apply it to solve compressible multi-fluid flows with large density ratio (1000) and strong shock wave (the pressure ratio is 10,000) interaction with the interface.

Saturation‐Dependent Hydraulic Conductivity Anisotropy for Multifluid Systems in Porous Media

The hydraulic conductivity of unsaturated anisotropic soils has recently been described with a tensorial connectivity–tortuosity (TCT) concept. We extend this concept to unsaturated porous media with two or three immiscible fluids. Mathematical expressions to describe the conductivity of each fluid in anisotropic porous media under unsaturated condition are derived in the form of symmetric second order tensors. The theory is applicable to the generalized hydraulic conductivity model and compatible types of saturation–pressure formulation. The extended model shows that the anisotropic coefficient of any one of the fluids is independent of the saturation of other fluids. Synthetic Miller‐similar soils generated to be anisotropic were defined by allowing the saturated hydraulic conductivity to have different correlation ranges for different directions of flow. The extended TCT concept was tested using synthetic soils with four levels of heterogeneity and four levels of anisotropy. Numerical experiments of infiltration of two liquid phases, water and the nonaqueous phase liquid (NAPL) carbon tetrachloride, were performed to test the extended model. The results show that, similar to water in a two‐fluid (air–water) system, NAPL retention curves in a three‐fluid (air–NAPL–water) system were independent of flow direction but dependent on soil heterogeneity, whereas the connectivity–tortuosity coefficients are functions of both soil heterogeneity and anisotropy. The extended TCT model accurately describes unsaturated hydraulic functions of anisotropic soils and can be combined into commonly used relative permeability functions for use in multifluid flow and transport numerical simulations.

Behavior of a Viscous LNAPL Under Variable Water Table Conditions

An intermediate-scale experiment in a 1.02-m-long, 0.75-m-high, and 0.05-m-wide flow cell was conducted to investigate the behavior of a viscous LNAPL under variable water table conditions. Two viscous LNAPL volumes (0.4 L) were released, one week apart, from a small source zone on top of the flow cell into a partly saturated, homogenously packed porous medium. Following a redistribution period of 30 days after the second release, the water table was increased 0.5 m in 50 minutes. After the water table rise, viscous LNAPL behavior was monitored for an additional 45 days. Fluid saturation scans were obtained periodically with a fully automated dual-energy gamma radiation system. Results show that both spills follow similar paths downwards. Within two hours after the first LNAPL arrival, the capillary fringe was reduced across the cell by approximately 0.04 m (22%). This reduction is directly related to the decrease in the air-water surface tension from 0.072 to 0.057 N/m. LNAPL drainage from the unsaturated zone was relatively slow and a considerable residual LNAPL saturation was observed after 30 days of drainage. Most of the mobile LNAPL moved into the capillary fringe during this period. After a rapid 0.5 m water table rise, the LNAPL moved up in a delayed fashion. The LNAPL used the same path upwards as it used coming down during the infiltration phase. After 45 days, the LNAPL had moved up only approximately 0.2 m. Since the LNAPL had only moved up a limited amount, nonwetting fluid entrapment was limited. The experiment was simulated using the STOMP multifluid flow simulator, which includes entrapped and residual LNAPL saturation formation. A comparison indicates that the simulator is able to predict the observed phenomena well, including residual saturation formation in the vadose zone, and limited upward LNAPL movement after the water table rise. The results of this experiment show that viscous mobile LNAPL, subject to variable water table conditions, does not necessarily float on the water table and may not appear in an observation well.

Multi-physics treatment in the vicinity of arbitrarily deformable gas–liquid interfaces

A novel three-dimensional problem formulation is introduced for the simulation of turbulent interfacial multi-fluid flows. The strategy is built around the large eddy simulation (LES) concept, and can be employed for interfacial heat and mass transfer problems in which use can be made of either scalar transfer correlations, or exact mass/energy jump conditions. This multi-physics treatment capability at arbitrarily deformable interfaces translates into two main features: (i) a reconstructed distance function (RDF) is introduced to define a level-set interface-normal length scale, and (ii) an interfacial shear velocity is defined on the distance function support for further use in near-interface transport models. The solution algorithm uses VOF with piecewise planar interface reconstructions on a twice-as-fine mesh, and infers the convective mass fluxes from the interface solution to promote momentum conservation. The interfacial shear velocity defined on the distance function support is introduced to accommodate the asymptotic behaviour of turbulence approaching the interface in a proximity-dependent manner. Provided with highly accurate distance function data, the scheme generates near-interface damping functions that are second-order accurate and independent of interface orientation. The damping was found to be affected by errors introduced into shear velocity estimates by the high-frequency errors in the RDF scheme near the interface. The methodology has been applied for the simulation of a wave breaking scenario featuring multiple modes and interfacial length scales.

A Comparison of Models Describing Residual NAPL Formation in the Vadose Zone

A major shortcoming of multifluid flow simulators is the inability to predict the retention of nonaqueous phase liquid (NAPL) in the vadose zone after long drainage periods. Recently, three theoretical models—the Wipfler and van der Zee (WVDZ) model; the Van Geel and Roy (VGR) model; and the Lenhard, Oostrom, and Dane (LOD) model—have been proposed for describing residual NAPL formation. The WVDZ model assumes a critical total liquid saturation below which all NAPL becomes residual. The VGR and LOD models are extensions of an existing hysteretic relative permeability–saturation–capillary pressure model and assume formation of residual NAPL during NAPL drainage and imbibition, respectively. In this study, we compared model predictions against results of a series of static pressure cell experiments. We found no experimental evidence supporting the WVDZ concept of a critical total liquid saturation. The other two models yielded reasonable predictions. The VGR and LOD models were then incorporated into a multifluid flow simulator, and simulations of two transient column experiments were conducted. Both models performed considerably better than simulations without considering the formation of residual NAPL, underwriting the importance of incorporating this process in simulators. Although the VGR and LOD models are based on different conceptual models, no clear performance differences could be observed when simulation results were compared against the transient experimental data.

A front-tracking/front-capturing method for the simulation of 3D multi-fluid flows with free surfaces

A method for simulating incompressible, imiscible, unsteady, Newtonian, multi-fluid flows with free surfaces is described. A sharp interface separates fluids of different density and viscosity. Surface and interfacial tensions are also considered and the required curvature is geometrically approximated at the fronts by a least squares quadratic fitting. To remove small undulations at the fronts, a mass-conserving filter is employed. The numerical method employed to solve the Navier–Stokes equations is based on the GENSMAC-3D front-tracking method. The velocity field is computed using a finite-difference scheme on an Eulerian grid. The free-surface and the interfaces are represented by an unstructured Lagrangian grid moving through an Eulerian grid. The method was validated by comparing the numerical results with analytical results for a number of simple problems. Complex numerical simulations show the capability and emphasize the robustness of this new method.