Incompressible Two-phase Flow Research Articles

High-Order Bound-Preserving Local Discontinuous Galerkin Methods for Incompressible and Immiscible Two-Phase Flows in Porous Media

High-Order Bound-Preserving Local Discontinuous Galerkin Methods for Incompressible and Immiscible Two-Phase Flows in Porous Media

A Compressible Formulation of the One-Fluid Model for Two-Phase Flows

In this paper, we introduce a compressible formulation for dealing with 2D/3D compressible interfacial flows. It integrates a monolithic solver to achieve robust velocity–pressure coupling, ensuring precision and stability across diverse fluid flow conditions, including incompressible and compressible single-phase and two-phase flows. Validation of the model is conducted through various test scenarios, including Sod’s shock tube problem, isothermal viscous two-phase flows without capillary effects, and the impact of drops on viscous liquid films. The results highlight the ability of the scheme to handle compressible flow situations with capillary effects, which are important in computational fluid dynamics (CFD).

Immiscible Two-Phase Flow in Porous Media: Effective Rheology in the Continuum Limit

AbstractWe consider steady-state immiscible and incompressible two-phase flow in porous media. It is becoming increasingly clear that there is a flow regime where the volumetric flow rate depends on the pressure gradient as a power law with an exponent larger than one. This occurs when the capillary forces and viscous forces compete. At higher flow rates, where the viscous forces dominate, the volumetric flow rate depends linearly on the pressure gradient. This means that there is a crossover pressure gradient that separates these two flow regimes. At small enough pressure gradient, the capillary forces dominate. If one or both of the immiscible fluids percolate, the volumetric flow rate will then depend linearly on the pressure gradient as the interfaces will not move. If none of the fluids percolate, there will be a minimum pressure gradient threshold to mobilize the interfaces and thereby get the fluids moving. We now imagine a core sample of a given size. The question we pose is what happens to the crossover pressure gradient that separates the power-law regime from the high-flow rate linear regime and the threshold pressure gradient that blocks the flow at low pressure gradients when the size of the core sample is increased. Based on analytical calculations using the capillary bundle model and on numerical simulations using a dynamical pore-network model, we find that the crossover pressure gradient and the threshold pressure gradient decrease with two distinct power laws in the size. This means that the power-law regime disappears in the continuum limit where the pores are infinitely small compared to the sample size.

Fully higher-order coupling of finite element and level set methods for two-phase flow with a new explicit projection method

This work presents the development of a higher-order finite element method for modeling unsteady, incompressible and inviscid two-phase flows using the level set method. A new explicit projection method is introduced to solve the incompressible Euler equation, where the pressure and the velocity fields are solved separately. To implement high order temporal schemes, the Hodge decomposition in this projection method is not applied to the velocity but to its time derivative. A discontinuous Galerkin discretization is used to solve both the level set and the nonlinear advection terms present in the Euler equations, while a Galerkin discretization is used to solve the pressure Poisson equation. The present method is a fully high order numerical approach, both in space and temporal dimensions. Numerical results provided by several 2D benchmark test cases in naval hydrodynamics (i.e hydrostatic equilibrium, sloshing and dam break) validate the computational code, illustrating the power of the proposed method. Compared to the CFD codes based on low order VOF/finite-volume methods, the present approach demonstrates greater accuracy on coarse meshes, leading to reduced CPU time.

The subdivision-based IGA-EIEQ numerical scheme for the Navier–Stokes equations coupled with Cahn–Hilliard phase-field model of two-phase incompressible flow on complex curved surfaces

We develop an accurate and robust numerical scheme for solving the incompressible hydrodynamically coupled Cahn–Hilliard system of the two-phase fluid flow system on complex surfaces. Our algorithm leverages a number of efficient techniques, including the subdivision-based isogeometric analysis (IGA) method for spatial discretization, the explicit Invariant Energy Quadratization (EIEQ) method for linearizing nonlinear potentials, the Zero-Energy-Contribution (ZEC) method for decoupling, and the projection method for the Navier–Stokes equation to facilitate fully decoupled type implementations. The integration of these methodologies results in a fully discrete scheme with desired properties such as linearity, second-order temporal accuracy, full decoupling, and unconditional energy stability. The implementation of the scheme is straightforward, requiring the solution of a few elliptic equations with constant coefficients at each time step. The rigorous stability proof of unconditional energy stability and the implementation procedure are given in detail. Numerous numerical simulations on complex curved surfaces are carried out to verify the effectiveness of the proposed numerical scheme.

A lattice Boltzmann model for incompressible gas and liquid two-phase flows combined with free-surface method

A novel lattice Boltzmann (LB) model is proposed to study the gas and liquid two-phase flows with large density and viscosity ratios. In the model, both the gas and liquid phases are considered as viscous incompressible fluids, which are governed separately by the two-relaxation-time LB equations. They are coupled by a momentum exchange method at the interface. The interaction between the gas and liquid phases is explicitly described and naturally involved in the model. The interfacial conditions in the model are validated by the benchmarks of the layered Poiseuille flow and the Laplace law. The feasibility of combining this model with the bubble model and the wetting scheme is proven through transient flow problems of single bubble rising and capillary intrusion. The validity of this model is confirmed by more complex flows including solid–liquid–gas coupling and droplet breaking problems by simulating shearing a droplet on a substrate and a droplet falling on a liquid film. The results demonstrate that the present model can be used to describe both the gas and the liquid flows. This work provides a solution to model the simulation of the dynamical behaviors of multi-phase flows.

A conservative sharp interface method for incompressible viscous two-phase flows with topology changes and large density ratio

Topology change of interfaces often leads to unresolved interface structures such as tiny droplets/bubbles, which requires the regularization of the order parameter and consequently causes artificial modification of volume fraction of the fluids. This could give rise to unphysical oscillation of pressure and cause the failure of computation, in particular for two-phase flows involving large density ratio. We develop a robust conservative sharp-interface method for incompressible viscous two-phase flows with topology changes and large density ratio. The method consists of a cut-cell meshing approach that dynamically generates unstructured meshes in the vicinity of the interface, a second-order finite volume scheme in the Arbitrary Lagrangian–Eulerian framework, and a projection method. The jump conditions at the interfaces are incorporated into the calculation of fluxes at the cell faces of the unstructured meshes that coincide with the reconstructed interfaces, to ensure that they are strictly satisfied at the interface. Unlike a pseudo-compressible formulation used in the previous study, the projection method prevents the oscillation of pressure in the presence of topology change. Special treatments are proposed to further improve the robustness of the method when dealing with unresolved interface structures, including removal of tiny droplets, adjustment of spatial discretization and restriction of the maximum curvature calculated. In such a way, two-phase flows involving topology change and large density ratio can be numerically resolved in a sharp-interface, efficient and robust manner. The performance of the proposed method is systematically examined by a series of numerical tests, including spurious currents, oscillation of a suspended droplet, droplet deformation in shear flows, breakup of a liquid thread, rising bubbles and partial coalescence of a droplet with a pool. The numerical results are compared against analytical solution, benchmark results or experimental data, and good agreement has been achieved.

Discretization schemes for the two simplified global double porosity models of immiscible incompressible two-phase flow

We present the discretization schemes for the two simplified homogenized models of immiscible incompressible two-phase flow in double porosity media with thin fractures. The two models were derived previously by the authors by different linearizations of the nonlinear local problem called the imbibition equation which appears in the homogenized model after passage to the limit as ε → 0. The models are fully homogenized with the matrix-fracture source terms expressed as a convolution.

Numerical simulation for the conserved Allen–Cahn phase field model of two-phase incompressible flows by an efficient dimension splitting method

A second-order, stable, and efficient numerical method is proposed for the simulation of high-dimensional two-phase incompressible flow. The model considers the coupled system of incompressible Navier–Stokes equations with variable density and the conserved Allen–Cahn equations. To save storage and computing costs, a dimension splitting scheme is proposed, which combines a second-order pressure correction method with a splitting effect and an alternating direction implicit method to split a two- or three-dimensional problem into a series of one-dimensional sub-problems that can be solved in parallel. By reconstructing the time derivative terms, the variable coefficient velocity equations are transformed into constant coefficient linear equations. In addition, a second-order stabilization term is added to enhance stability. Finally, a post-processing method is proposed which can preserve both the maximum bound principle and the mass conservation. For space discretization, the weighted essentially non-oscillatory scheme and central difference scheme are used on the marker and cell grid, which can effectively avoid numerical oscillation caused by high Reynolds numbers. Numerical experiments verify the convergence, stability, parallel efficiency, and effectiveness of post-processing, and demonstrate some simulations of two-phase incompressible flows.

Lattice Boltzmann method for studying the dynamics of a single rising bubble in shear-thickening power-law fluids

Bubble motion in non-Newtonian fluids is widespread in various industrial processes such as crude oil extraction, enhancement of boiling heat transfer, CO<sub>2</sub> sequestration and wastewater treatments. Systems containing non-Newtonian liquids, as opposed to Newtonian liquids, have shear-dependent viscosity, which can alter the hydrodynamic characteristics of the bubbles, such as their size, deformation, instability, terminal velocity, shear rate, and ultimately affect the bubble rising behaviors. In this work, The dynamic behavior of bubble rising in a shear-thickened fluid was studied using an incompressible lattice Boltzmann non-Newtonian gas-liquid two-phase flow model. The effects of the rheological exponent <i>n</i>, the Eötvös number (<i>Eo</i>), and the Galilei number (<i>Ga</i>) on the bubble deformation, terminal velocity, and the shear rate were investigated. The numerical results show that the degree of bubble deformation increases as <i>Eo</i> grows, and the effect of <i>n</i> on bubble deformation degree is correlated with <i>Ga</i>. On the other hand, the terminal velocity of the bubbles increases monotonically and nonlinearly with <i>Ga</i> for a given <i>Eo</i> and <i>n</i>, and the effect of <i>n</i> on the terminal velocity of the bubbles gets stronger as <i>Ga</i> increases. When <i>Ga</i> is fixed and small, the terminal velocity of the bubble increases and then decreases with the increase of <i>n</i> at small <i>Eo</i>, and increases with the increase of <i>n</i> when <i>Eo</i> is large; but when <i>Ga</i> is fixed and large, the terminal velocity of the bubble increases with the increase of <i>n</i> in a more uniform manner. In addition, regions with high shear rates can be found near the left and right ends of the bubble. The size of these regions increases as both <i>Eo</i> and <i>Ga</i>, and it first increases and then decreases as <i>n</i> increases. Finally, the orthogonal experimental method was used to obtain the impact of the aforementioned three factors on the shear rate and terminal velocity. For the shear rate, the order of influence is <i>n</i>, <i>Ga</i> and <i>Eo</i> in descending order. For the terminal velocity, <i>Ga</i> has the greatest influence, followed by <i>n</i>, and <i>Eo</i> has the least influence.

The extended Discontinuous Galerkin method for two-phase flows with evaporation

In this work, we present an extended Discontinuous Galerkin method for simulating transient, incompressible two-phase flows, which include heat transfer and thermally driven evaporation at the interface of single-component systems. This expands our previous work to include the consideration of non-material interfaces and a coupling between velocities and temperature gradients at the interface. The phase boundary is represented by the zero-set of a level set function, while effects due to surface tension are treated by the Laplace-Beltrami formulation. This sharp interface model allows for a sub-grid accurate representation of the solution fields. By using compactly supported polynomial solutions, discontinuities at the interface can be sharply represented without employing additional reconstruction schemes. The approach is validated through well-known evaporation test cases. This includes two 1D test cases, known as Stefan and Sucking problem, a 2D film boiling and finally the 3D growth of a vapor bubble, known as Scriven test case.

A conservative sharp interface method for two-dimensional incompressible two-phase flows with phase change

A conservative sharp interface method is proposed in this work to simulate two-dimensional/axisymmetric incompressible two-phase flows with phase change. In this method, we use the cut cell method to generate unstructured meshes near the interface, of which the cell edges overlap with the interface at each time step. On such mesh, the mass and heat transfer during phase change and all the jump conditions can be incorporated into the calculation of fluxes at the cell edges, to ensure that they are strictly satisfied at the interface in a sharp manner. The governing equations, including the incompressible Navier–Stokes equations, heat equation, and vapor mass fraction equation, are discretized by a second-order finite volume method in the arbitrary Lagrangian–Eulerian framework. To well couple the mass, heat, momentum, and interface evolution, the solution procedure is carefully designed and performed with several techniques. In such a way, the sharp discontinuity of the velocity, stress, temperature gradient, and vapor fraction, caused by the mass/heat transfer during phase change, can be simulated accurately and robustly. The performance of this method is systematically examined by cases of phase change at or below the saturated temperature, including vapor bubble in superheated liquid, film boiling, droplet evaporation at different relative humidity conditions, droplet evaporation under gravity, and droplet evaporation under forced convection. The applicability of the present method for incompressible two-phase flows with phase change is well demonstrated by comparing the numerical results with the benchmark, theoretical or experimental ones.

Quantitative comparison of dam break and bubble flows based on CF-VOF and IR-VOF schemes

In this study, an incompressible and immiscible two-phase flow solver was employed to assess the accuracy of phase-interface capturing methods in simulating two-phase flows. Three benchmark investigations were conducted involving two dam break flows with deformable free surfaces and the upward motion of a gas bubble through a liquid column. Four phase-interface capturing schemes, namely, MULES, HRIC and CICSAM from the Color Function Volume of Fluid (CF-VoF), and Interface Reconstruction Volume of Fluid (IR-VoF), were utilized. Comparative analysis of the phase interfaces indicated that CICSAM and IsoAdvector consistently provided sharp interfaces across all cases when compared to MULES and HRIC. The mass conservation performance of all schemes excelled in the dam-break case with relatively simple free-surface development and the rising bubble case. An examination of the local pressure distribution over time revealed that methods yielding diffusive interfaces produced inaccurate results. Thus, the importance of employing denser grids or schemes that yield sharper interfaces, such as CICSAM and IsoAdvector, to enhance simulation accuracy was underscored. Overall, the results of this analysis confirm that the IsoAdvector scheme, which provides sharp interfaces, is a reliable option for simulating incompressible immiscible two-phase flows.

TwoPhaseFlow: A Framework for Developing Two Phase Flow Solvers in OpenFOAM

We present a new OpenFOAM based open-source framework, TwoPhaseFlow, enabling fast implementation and testing of new phase change and surface tension force models for two-phase flows including interfacial heat and mass transfer. Capitalizing on the runtime-selection mechanism in OpenFOAM, the new models can easily be selected and benchmarked against analytical solutions and existing models. The framework currently includes the following three interface curvature calculation methods for surface tension: 1) the height function method, 2) the parabolic fit method and 3) the reconstructed distance function method. As for phase change, two models are available: 1) Interface heat resistance and 2) direct heat flux. These can be combined in three solvers: 1) InterFlow for isothermal, incompressible two-phase flow, 2) compressibleInterFlow for compressible, non-isothermal two-phase flow and 3) multiRegionPhaseChangeFlow for compressible, non-isothermal two-phase flow with conjugated heat transfer. By design, addition of new models and solvers is straightforward and users are encouraged to contribute their specific models, solvers, and validation cases to the library.

Some remarks on simplified double porosity model of immiscible incompressible two-phase flow

The paper is devoted to the derivation, by linearization, of simplified homogenized models of an immiscible incompressible two-phase flow in double porosity media in the case of thin fissures. In a simplified double porosity model derived previously by the authors the matrix-fracture source term is approximated by a convolution type source term. This approach enables to exclude the cell problem, in form of the imbibition equation, from the global double porosity model. In this paper we propose a new linear version of the imbibition equation which leads to a new simplified double porosity model. We also present numerical simulations which show that the matrix-fracture exchange term based on this new linearization procedure gives a better approximation of the exact one than the corresponding exchange term obtained earlier by the authors.

Investigations on the shallow water wave attenuation over continuous porous oyster reef-like structure

Observations over the decades have found that oyster reefs in shallow estuarine environments have the potential to protect shorelines or offshore structures by dissipating wave energy. Laboratory experiments were conducted to observe the shallow water wave attenuation in different terrains of the continuous porous oyster reef-like structures. To better reproduce this process, a numerical wave flume that is within the OpenFOAM scheme was developed in this study. The volume-averaged Reynolds-Averaged Navier–Stokes (VARANS) equations were solved for two-phase incompressible flow around the porous structures and the VOF method was used for tracking the free surface. The model was first validated by the present laboratory experiments. Subsequently, dimensional analysis was conducted to address the factors influencing wave attenuation over continuous porous reef-like structures. The model was then applied to evaluate the impacts of wave factors and morphological factors of porous oyster reef-like structures on the wave attenuation. By the end, an empirical formula was proposed to predict the wave transmission coefficient based on the numerical results and the experimental data.

Inverse asymptotic treatment: Capturing discontinuities in fluid flows via equation modification

A major challenge in developing accurate and robust numerical solutions to multi-physics problems is to correctly model evolving discontinuities in field quantities. These manifest themselves as interfaces between different phases in multi-phase flows, or as shock and contact discontinuities in (either single- or multi-phase) compressible flows. A plethora of bespoke discretization schemes have been developed to capture both types of discontinuities in physics-based simulations. However, when a quick response is required to rapidly emerging challenges such as the need to design novel hypersonic vehicles, the complexity of such implementations impedes a swift transition from problem formulation to computation. This is exacerbated by the need to compose multiple interacting physics, which may cause specialized schemes to break unexpectedly as governing equations need to be changed and/or added. We introduce “inverse asymptotic treatment” (IAT) as a unified framework for capturing discontinuities in fluid flows that enables building directly computable models based on standard, off-the-shelf numerics. By capturing discontinuities through modifications at the level of the governing equations, rather than via reliance on specialized discretization schemes, IAT can seemlessly handle additional physics and thus makes it easier for novice end users to quickly obtain numerical results for a variety of multi-physics scenarios. This strategy also facilitates the use of a “multi-physics” compiler that automates the conversion of the modified PDEs to numerical source code readable by popular computational frameworks like OpenFOAM. We outline IAT in the context of phase-field modeling of two-phase incompressible flows, and then demonstrate its generality by showing how localized artificial diffusivity (LAD) methods for single-phase compressible flows can be viewed as instances of IAT. Through the real-world example of a laminar hypersonic compression corner, we illustrate IAT’s ability to, within a span of just a few months, generate a directly computable model whose wall metrics predictions for sufficiently small corner angles come close to that of NASA’s state-of-the-art VULCAN-CFD solver. Finally, we propose a novel LAD approach via “reverse-engineered” PDE modifications, inspired by total variation diminishing (TVD) flux limiters, to eliminate the problem-dependent parameter tuning that plagues traditional LAD. Through canonical numerical tests, we demonstrate that this “limiter-inspired” LAD approach, when combined with second-order central differencing, can robustly and accurately model compressible flows.

On a Navier–Stokes–Cahn–Hilliard system for viscous incompressible two-phase flows with chemotaxis, active transport and reaction

On a Navier–Stokes–Cahn–Hilliard system for viscous incompressible two-phase flows with chemotaxis, active transport and reaction

Fully discrete, decoupled and energy-stable Fourier-Spectral numerical scheme for the nonlocal Cahn–Hilliard equation coupled with Navier–Stokes/Darcy flow regime of two-phase incompressible flows

In this paper, we introduce fully discrete Fourier-Spectral numerical scheme to solve the nonlocal Cahn–Hilliard equation coupled with Navier–Stokes/Darcy equations, which represent a phase-field model for two-phase incompressible flow in either the free flow regime or a Hele-Shaw cell. The proposed scheme achieves full decoupling while maintaining linearity and energy stability through a combination of the Scalar Auxiliary Variable (SAV) method, which discretizes the nonlinear potential, and the “Zero-Energy-Contribution” (ZEC) method, which handles the coupled nonlinear advective/surface tension terms. The efficiency of this scheme is attributed to its linear decoupling structure and the fact that it requires only a few elliptic equations with constant coefficients to be solved at each time step. We rigorously establish the scheme’s unconditional energy stability. Further, some numerical simulations are provided in both 2D and 3D to show its effectiveness, including its accuracy, stability, and efficiency.

Flow and mixing dynamics in face-to-face and rear-end collisions of pairs of equal-sized droplets

In this work, the behaviors of pairs of equal-sized droplets in rear-end and face-to-face collisions were simulated using the improved lattice Boltzmann method for incompressible two-phase flows. First, the time evolution of the droplet shape was investigated by tracing colored particles, and this was compared between the rear-end and face-to-face collisions. For collinear collisions, the droplet shapes in the rear-end collisions were found to be similar to those in the face-to-face collisions. However, the behaviors of the tracer particles were different: the droplets in the rear-end collisions mixed more easily than those in the face-to-face collisions. For offset collisions, it was found that the rolling motion of the coalesced droplet accelerates the mixing inside it in both face-to-face and rear-end collisions. A new index—the total mixing intensity—was introduced, and the droplet mixing can be quantitatively evaluated by calculating its value. The results indicate that the droplet mixing process of a collinear collision can be characterized by the velocity ratio, which is defined as the ratio of the center-of-mass velocity to the relative impact velocity.