Incompressible Two-phase Flow Research Articles

A novel fully-decoupled, linear, and unconditionally energy-stable scheme of the conserved Allen–Cahn phase-field model of a two-phase incompressible flow system with variable density and viscosity

In this article, we first use the nonlocal conserved Allen–Cahn dynamics to reformulate the phase-field model with variable density and viscosity for the two-phase incompressible flow system, and then construct a novel fully-decoupled, linear, and unconditionally energy stable scheme to solve the model. The development of the scheme is based on a combination of a novel decoupling technique, the penalty method of the Navier–Stokes equation, the quadratization method of linearizing the double-well potential, and the operator Strang-splitting method. The scheme is highly efficient, and it only needs to solve a series of fully-decoupled linear elliptic equations at each time step, in which the Allen–Cahn equation and the pressure Poisson equation only have constant coefficients. We rigorously prove the unconditional energy stability of the scheme and conduct numerical simulations in 2D and 3D to demonstrate the accuracy, stability, and effectiveness of the scheme numerically.

An energy stable linear numerical method for thermodynamically consistent modeling of two-phase incompressible flow in porous media

In this paper, we consider numerical approximation of a thermodynamically consistent model of two-phase flow in porous media, which obeys an intrinsic energy dissipation law. The model under consideration is newly-developed, so there is no energy stable numerical scheme proposed for it at present. This model consists of two nonlinear degenerate parabolic equations and a saturation constraint, but lacking an independent equation for the pressure, so for the purpose of designing efficient numerical scheme, we reformulate the model forms as well as the free energy function, and further prove the corresponding energy dissipation inequality. Based on the alternative reformulations, using the invariant energy quadratization approach and subtle semi-implicit treatments for the pressure and saturation, we for the first time propose a linear and energy stable numerical method for this model. The fully discrete scheme is devised combining the upwind approach for the phase mobilities and the cell-centered finite difference method. The unique solvability of numerical solutions and unconditional energy stability are rigorously proved for both the semi-discrete time scheme and the fully discrete scheme. Moreover, the scheme can guarantee the local mass conservation for both phases. We also show that the upwind mobility approach plays an essential role in preserving energy stability of the fully discrete scheme. Numerical results are presented to demonstrate the performance of the proposed scheme.

Numerical modelling of wave transformation and runup over rough fringing reefs using VARANS equations

Computational Fluid Dynamics (CFD) models are gaining popularity recently in efforts to understand the complex dynamics of wave interaction with coral reef profiles. To provide an alternative approach to the current CFD models to account for the reef surface with large bottom roughness, a 3D numerical wave tank based on the CFD tool OpenFOAM® is developed in this study. The VARANS equations are solved for two-phase incompressible flow with the k-ω SST model for the turbulence closure and VOF method for tracking the free surface. The reef surface with high friction is modeled by using a porous media model in the VARANS equations. The numerical model is first validated with laboratory experiments. The model is then applied to evaluate the impacts of forcing condition (wave and water level), reef morphology (reef-flat width, fore-reef slope, and back-reef slope), and reef roughness on wave transformation and wave runup. Moreover, an empirical formula is derived from the numerical results to predict wave runup. Subsequently, the simulated cross-reef variations of energy dissipation, wave spectrum, wave shape parameters (skewness and asymmetry) and nonlinearity parameter (Ursell number) are discussed for both the smooth and rough reefs. Finally, the cross-reef variations of the TKE and its dissipation rate are also examined for both reefs via the numerical simulations.

A consistent and conservative Phase-Field model for thermo-gas-liquid-solid flows including liquid-solid phase change

In the present study, a consistent and conservative Phase-Field model is developed to study thermo-gas-liquid-solid flows with liquid-solid phase change. The proposed model is derived with the help of the consistency conditions and exactly reduces to the consistent and conservative Phase-Field method for incompressible two-phase flows, the fictitious domain Brinkman penalization (FD/BP) method for fluid-structure interactions, and the Phase-Field model of solidification of pure material. It honors the mass conservation, defines the volume fractions of individual phases unambiguously, and therefore captures the volume change due to phase change. The momentum is conserved when the solid phase is absent, but it changes when the solid phase appears due to the no-slip condition at the solid boundary. The proposed model also conserves the energy, preserves the temperature equilibrium, and is Galilean invariant. A novel continuous surface tension force to confine its contribution at the gas-liquid interface and a drag force modified from the Carman-Kozeny equation to reduce solid velocity to zero are proposed. The issue of initiating phase change in the original Phase-Field model of solidification is addressed by physically modifying the interpolation function. The corresponding consistent scheme is developed to solve the model, and the numerical results agree well with the analytical solutions and the existing experimental and numerical data. Two challenging problems having a wide range of material properties and complex dynamics are conducted to demonstrate the capability of the proposed model.

Application of a high-order MP scheme to computation of multi-phase flows with heat and mass transfer

Reliable prediction of unsteady multiphase flows with mass changes, large density ratios, and complex flow features is very challenging due to the need to accurately resolve both complex interface deformations and steep gradients associated with material interfaces. This paper describes an application of our recently developed high-order monotonicity preserving (MP) scheme – a reliable and robust scheme for long-term computations of flows contain complex smooth structures interspersed with discontinuities (Ha and Lee, 2020) − to the solution of unsteady incompressible two-phase flows in the presence of heat and mass changes. The flow is governed by the multiphase Navier–Stokes equations, which comprises of the mixture continuity equation, the mixture momentum equations, the mixture energy equation, and the continuity equation for the individual phase. The numerical discretization is implemented on generalized curvilinear grids for facilitating flow simulations of the complex geometries. Our high-order MP scheme is applied to the calculations of the convective and surface tension terms. To obtain an effectively simultaneous solution of all equations, a novel dual-time stepping approach in conjunction with an implicit Runge–Kutta three-stage scheme is proposed. The reliability of the proposed numerical algorithm is then examined for a 1D Stefan condensation problem and two condensing flows including single bubble and steam-jet condensation problems. Particularly, for the single bubble condensation test, the bubble characteristics at different time instants are evaluated both qualitatively and quantitatively. For the steam-jet condensation test, the capability for capturing the evolution of interface deformations under different thermal and hydrodynamic flow conditions is highlighted. Finally, the effects of subcooling degree and steam mass flow rate on heat and mass transfer characteristics are investigated and discussed.

A decoupled, stable, and linear FEM for a phase-field model of variable density two-phase incompressible surface flow

The paper considers a thermodynamically consistent phase-field model of a two-phase flow of incompressible viscous fluids. The model allows for a non-linear dependence of the fluid density on the phase-field order parameter. Driven by applications in biomembrane studies, the model is written for tangential flows of fluids constrained to a surface and consists of (surface) Navier–Stokes–Cahn–Hilliard type equations. We apply an unfitted finite element method to discretize the system and introduce a fully discrete time-stepping scheme with the following properties: (i) the scheme decouples the fluid and phase-field equation solvers at each time step, (ii) the resulting two algebraic systems are linear, and (iii) the numerical solution satisfies the same stability bound as the solution of the original system under some restrictions on the discretization parameters. Numerical examples are provided to demonstrate the stability, accuracy, and overall efficiency of the approach. Our computational study of several two-phase surface flows reveals some interesting dependencies of flow statistics on the geometry.

New efficient time-stepping schemes for the Navier–Stokes–Cahn–Hilliard equations

In this paper, we construct and analyze a class of efficient time-stepping schemes for the Navier–Stokes–Cahn–Hilliard equation system stemming from phase field modeling of two-phase incompressible flows. The proposed schemes are based on an auxiliary variable approach for the Cahn–Hilliard equation and a delicate treatment of the terms coupling the Navier–Stokes equation and the Cahn–Hilliard equation. The former has been found to be an efficient tool for numerical integration of the Cahn–Hilliard equation in recent years. A detailed comparison with existing schemes is given, and the advantage of the new schemes are emphasized. In the theoretical aspect we rigorously prove that the designed schemes are unconditionally stable in the sense that some kind of energy remains bounded during the time stepping. In the implementation, we show how the schemes can be reformulated into six linear second-order elliptic equations with constant coefficients. The efficiency of the proposed algorithm is verified through several numerical examples.

A linear second-order in time unconditionally energy stable finite element scheme for a Cahn–Hilliard phase-field model for two-phase incompressible flow of variable densities

We propose a novel second-order BDF time stepping method of variable time step sizes combined with a classical residual-based stabilized finite element spatial discretization using the Streameline-Upwind Petrov–Galerkin (SUPG)/pressure stabilization Petrov–Galerkin (PSPG)/grad-div stabilization for solving the phase–field model for two-phase incompressible flow of different densities and viscosities in the advection dominated regime. In the case of uniform time step size and without extra stabilization, the scheme is shown to satisfy a discrete energy law. Benchmark test of the Rayleigh–Taylor instability under high Reynolds number and Péclect number demonstrates that the scheme captures details of the instability comparable to results in the literature by schemes based on sharp-interface models.

A parallel cell-centered adaptive level set framework for efficient simulation of two-phase flows with subcycling and non-subcycling

We develop a unified adaptive level set (LS) framework using the multi-level collocated grid for incompressible two-phase flows. This framework allows us to advance all variables level by level using either the subcycling or the non-subcycling method such that the data advancement on each level is fully decoupled. A series of synchronization operations are designed to keep the momentum and mass conserved across all levels. A multi-level re-initialization method for the LS function is also proposed, which greatly improves the mass conservation of the two-phase flows. The collocated grid allows the use of a single set of differential schemes and interpolation operations for all variables, which greatly simplifies the numerical implementation. The capability and robustness of the computational framework are validated by a variety of canonical problems, including the inviscid shear layer, gravity wave, rising bubble, and Rayleigh-Taylor instability. It is shown that the present multi-level scheme can accurately resolve the interfaces of the two-phase flows with gravitational and surface tension effects while having good momentum and energy conservation. At last, a three-dimensional dam breaking problem is simulated to show the efficiency and significant speedup of the proposed framework.

Global well-posedness of a Navier–Stokes–Cahn–Hilliard system with chemotaxis and singular potential in 2D

We study a diffuse interface model that describes the dynamics of incompressible two-phase flows with chemotaxis effects. This model also takes into account some significant mechanisms such as active transport and nonlocal interactions of Oono's type. The system under investigation couples the Navier–Stokes equations for the fluid velocity, a convective Cahn–Hilliard equation with physically relevant singular potential for the phase-field variable and an advection-diffusion-reaction equation for the nutrient density. For the initial boundary value problem in a smooth bounded domain Ω ⊂ R 2 , we first prove the existence and uniqueness of global strong solutions that are strictly separated from the pure states ±1 over time. Then we prove the continuous dependence with respect to the initial data and source term for the strong solution in energy norms. Finally, we show the propagation of regularity for global weak solutions.

Magnetic bubbles dynamics and heat transfer characteristics under influence of non-uniform magnetic fields

The computational study of transient immiscible and incompressible two-phase flows is one of the most common and desirable way for investigation of engineering phenomena and physics science. In the previous studies, generally bubbles current have been used as an active method for increasing heat transfer, however, due to existence of hydraulic boundary layers, the bubbles were not able to cross over this layer to thinning the thermal boundary layer and consequently the efficiency of this method was not very considerable. In this study, by considering potential of magnetic field, the effect of co-applying of external non uniform magnetic field and magnetic bubbles in enhancing the heat transfer efficiency in a 3-D tube has been investigated. The computational model consisted of the Navier–Stokes equation for liquid phase and VOF model for interface tracking are carried out by OpenFOAM. The external magnetic field has been considered non-uniform and time dependent. The results predicted that magnetic bubbles and external magnetic field due to their effect on thermal boundary layer increased significantly heat transfer and Nusselt number. Furthermore, results indicated magnetic bubbles can act as an active torbulators in the flow field and can be applied for increasing recirculation and secondary flow in the flow field. The average temperature and magnetic field over times for different cases have been discussed in the results.

A simple augmented IIM for 3D incompressible two-phase Stokes flows with interfaces and singular forces

In this paper, a Stokes solver on a MAC staggered grid for three-dimensional incompressible two-phase flows with interfaces and singular forces is presented. The main purpose is to develop the simplified IIM to ensure the accuracy near the interface represented by a level-set function. The idea of the simplified IIM is to apply Taylor expansion only along the normal direction and incorporate the jump conditions up to the second normal derivatives into the finite difference scheme. A Newton iteration method is applied to approximate the orthogonal projections for points near the interface. The augmented variables are introduced to decouple the coupled jump conditions of velocity and pressure. The determination of the augmented variables is finally reduced to a linear system, which is solved by the GMRES method. A CG-Uzawa type method is used to solve the discetized Stokes equation with some correction terms. Numerical results show that the proposed algorithm is efficient and can reach global second-order accuracy. The present algorithm can also preserve the area conservation and the discrete divergence free condition very well.

Eliminating spurious currents in phase-field-theory-based lattice Boltzmann equation for two-phase flows

Spurious currents are frequently observed near an interface in the equilibrium multiphase flow system by lattice Boltzmann equation (LBE). These unphysical phenomena are the result of force imbalance of LBE at a discrete level. In this paper, we develop a well-balanced Cahn–Hilliard equation-based LBE for incompressible two-phase flows. The effects of small initial perturbation of order parameter or dynamic pressure and nonisotropic discretization of gradient in force term on eliminating the spurious currents are investigated systematically. Numerical simulations including flat interface and stationary droplet problems are carried out to show the capability of present LBE for eliminating the spurious currents and its accuracy. The results predicted by the present LBE are compared with those by mixed isotropic discretizations scheme (MIDS) frequently used in the LBE community. Numerical results show that the initial perturbation of order parameter or dynamic pressure and nonisotropic discretization of gradient term has no significant effect on eliminating the spurious currents by present LBE, while the MIDS is sensitive to them.

Fully-discrete finite element numerical scheme with decoupling structure and energy stability for the Cahn–Hilliard phase-field model of two-phase incompressible flow system with variable density and viscosity

We construct a fully-discrete finite element numerical scheme for the Cahn–Hilliard phase-field model of the two-phase incompressible flow system with variable density and viscosity. The scheme is linear, decoupled, and unconditionally energy stable. Its key idea is to combine the penalty method of the Navier–Stokes equations with the Strang operator splitting method, and introduce several nonlocal variables and their ordinary differential equations to process coupled nonlinear terms. The scheme is highly efficient and it only needs to solve a series of completely independent linear elliptic equations at each time step, in which the Cahn–Hilliard equation and the pressure Poisson equation only have constant coefficients. We rigorously prove the unconditional energy stability and solvability of the scheme and carry out numerous accuracy/stability examples and various benchmark numerical simulations in 2D and 3D, including the Rayleigh–Taylor instability and rising/coalescence dynamics of bubbles to demonstrate the effectiveness of the scheme, numerically.

Gas-liquid two-phase flows simulation based on weakly compressible scheme with interface-adapted AMR method

A weakly compressible scheme using the interface-adapted adaptive mesh refinement (AMR) method is proposed for simulating low Mach number gas-liquid two-phase flows in large-scale and high-resolution simulations. A pressure evolution equation is derived from the iso-thermal process and low-Mach-number assumptions. We solve the Navier-Stokes equation and the pressure equation explicitly including acoustic waves. In the conventional semi-implicit method to solve incompressible two-phase flows, the computational cost of the linear solver increases as the problem size becomes larger. The conservative phase-field equation was solved using the finite volume method (FVM) to capture the gas-liquid interface while conserving the total mass of each phase. We developed GPU code for a tree-based AMR that greatly reduces the computational cost of assigning a high-resolution mesh around the moving interface region. The data structure of the memory pool has extra space allocated to prevent frequent memory allocations and deallocations. Results for the two-phase flow benchmark problems, 2D rising-bubble and dam-breaking onto dry floor problem, are in good agreement with the results of incompressible solvers. The two-phase flow simulation of the dam-breaking onto a wet floor problem maintained the shape interface regardless of whether the AMR method was used or not. We have computed a two-dimensional soap bubble growing from 16mm to 120mm in diameter with about 0.9 million total cells. The final thickness of the liquid film was about 0.1mm with the finest cell size being 0.0195mm, equivalent to 8,192×8,192 uniform cells. This computation showed that high-resolution two-phase flow simulation with the explicit scheme and AMR method reproduced the behavior of the thin liquid film.

Numerical analysis of fully discrete energy stable weak Galerkin finite element Scheme for a coupled Cahn-Hilliard-Navier-Stokes phase-field model

The Cahn-Hilliard phase-field model of two-phase incompressible flows, namely the Cahn-Hilliard-Navier-Stokes (CH-NS) problem represents the fundamental building blocks of hydrodynamic phase-field models for multiphase fluid flow dynamics. Since finding solution (numerically and theoretically) of the CH-NS system is non-trivial (owing to the coupling between the Navier-Stokes equation and the Cahn-Hilliard equation), in this paper, we propose and analyze a numerical scheme with the following properties: (1) first-order in time, (2) nonlinear, (3) fully coupled and (4) energy stable; for solving CH-NS system, in the framework of weak Galerkin (WG) method. More precisely, we employ the WG method, which uses discontinuous functions to construct the approximation space, and the first-order backward Euler (implicit) method for space and time discretizations, respectively. We first recall the corresponding variational formulation, and then summarize the main WG method ingredients that are required for our discrete analysis. In particular, in order to define the weak discrete bilinear (and trilinear) form, whose continuous version involves classical differential operators, we propose two well-known alternatives for gradient and divergence operators onto a suitable polynomial subspace. Next, we show that the weak global discrete bilinear form satisfies the hypotheses required by the Lax-Milgram lemma. In this way, we derive the associated a priori error estimates for phase field variable, chemical potential, velocity and further the pressure in L2 norm. Finally, several numerical results confirming the theoretical rates of convergence and illustrating the good performance of the method for simulation influence of surface tension, small density variations, and the driven cavity are presented.

Energy stable modeling of two-phase flow in porous media with fluid–fluid friction force using a Maxwell–Stefan–Darcy approach

Two-phase incompressible flow in porous media plays an important role in various fields including subsurface flow and oil reservoir engineering. Due to the interaction between two phases flowing through the pores, the fluid–fluid friction force may have a significant effect on each phase velocity. In this paper, we propose an energy stable (thermodynamically consistent) Maxwell–Stefan–Darcy model for two-phase flow in porous media, which accounts for the fluid–fluid friction. Different from the classical models of two-phase flow in porous media, the proposed model uses the free energy to characterize the capillarity effect. This allows us to employ the Maxwell–Stefan model to describe the relationships between the driving forces and the friction forces. The driving forces include the pressure gradient and chemical potential gradients, while both fluid–solid and fluid–fluid friction forces are taken into consideration. Thermodynamical consistency is the other interesting merit of the proposed model; that is, it satisfies an energy dissipation law and also obeys the famous Onsager's reciprocal principle. A linear semi-implicit numerical method is also developed to simulate the model. Numerical simulation results are provided to show that the fluid–fluid friction force can improve the oil recovery substantially during the oil displacement process.

Fluid-structure interactions in a flexible pipe conveying two-phase flow

We present three-dimensional direct numerical simulation (DNS) results of a horizontal flexible pipe conveying two-phase incompressible flow at moderate Reynolds numbers (400 ∼ 700). We find that the Reynolds number (Re) and the void fraction (α) essentially determine the onset of self-sustained oscillations in the two-fluid/pipe system. We employ a phase-field formulation to solve the Navier-Stokes coupled with the Cahn-Hilliard equation and the structure equation in an arbitrary Lagrangian Eulerian (ALE) framework. A spectral/hp element method is adopted for the fluid solver and a Galerkin method for the structure solver. We construct stability diagrams showing the transitions between different states of the two-fluid/pipe system. Specifically, we focus on the transition from a stable state to a planar oscillatory state and also to an out-of-plane oscillatory state as a function of four non-dimensional parameters, the Reynolds number Re, the void fraction α, the fluid-tension parameter Ip, and the averaged flow velocity U. By comparing simulation results with those for a pipe conveying single-phase flow under similar conditions, we find that the two-fluid/pipe system loses stability at a lower Reynolds number. We also find that, under similar Re,Ip, and U, at values of void fraction α∈0.6∼0.8, corresponding to the slug/churn flow regimes, the pipe experiences the largest vibrations, which is consistent with previous experimental and theoretical studies. We also perform a spectral analysis and identify the different excited modes, which are found to vary depending on the values of the void fraction α.

Numerical simulation of falling film flow hydrodynamics over round horizontal tubes

The performance of a falling film heat exchanger depends significantly on the flow characteristics and the associated heat transfer in the liquid film flowing over its horizontal tube surfaces. In this study, three-dimensional numerical simulations based on the Volume of Fluid method are conducted to investigate the hydrodynamic behavior of a falling film over horizontal round tubes with diameter and tube-to-tube spacings each of 16 mm. Using a mixture of water and ethylene glycol, the simulations comprise a range of flow rates (Re = 15 to 210) that covers all basic intertube flow modes i.e., droplet, jet (in-line and staggered), and sheet mode. InterFoam, a two-phase incompressible flow solver in OpenFOAM, has been used to simulate the flow field. The computational model agrees well with the available experimental data. The numerical results present the variation of the liquid film thickness and the film interface velocity over the tube surface for all the flow modes. In the droplet mode, the impact of the droplet on the tube produces liquid waves that interact with the neighboring waves and are seen to travel on the tube surface downstream of the impingement region and the region between the impingement sites. The movement of these waves causes a significant change in the film thickness over the tube surface, quantitatively, by more than 350% during a short interval of the droplet mode cycle. After the droplet impact, the liquid bridge formed in between the tubes undergoes necking and finally bifurcates under the Raleigh-Plateau instability mechanism at a location of spacing-to-tube-diameter ratio of 0.8 from the top surface of the impinging tube and produces satellite droplets in between the tubes. In the steady in-line jet mode, the liquid from the neighboring jets do not interact with each other, and the bulk liquid flows downstream of the impingement region. A dimple around the jet impingement region is formed where the film thickness changes by around 40%. The base of the impinging jets was seen to possess ripples due to the Raleigh-Plateau instability with wavelengths ranging from 0.3 - 1.0 times the capillary length scale. For the steady staggered jet mode, the liquid jet impinges on the tube surface, which then spreads rapidly and interacts with the adjacent jet liquid to develop crest and stable segments along the tube length direction. The ratio of the film thickness between the crest and the stable segments was found to be 1.7. Finally, for the sheet mode, due to the formation of surface waves in the liquid sheet, interfacial waves are seen to travel along the tube periphery with amplitudes of about ±20% of the nominal film thickness. These flow mode characteristics alter the film thickness and flow velocity distribution over the tube surface and are expected to play an essential role in the surface heat transfer behavior over round horizontal tubes. A set of correlations have been presented to predict the film thickness and interfacial velocity over the tube surface for each of the studied modes of flow with reasonable accuracy.

An efficient stabilized multiple auxiliary variables method for the Cahn–Hilliard–Darcy two-phase flow system

Herein, we develop a totally decoupled, linear, and temporally second-order accurate numerical scheme for the Cahn–Hilliard–Darcy system which models the two-phase incompressible fluid flows in porous medium or in Hele–Shaw cells. Our proposed scheme is based on a simple and efficient stabilized multiple scalar auxiliary variables (S-MSAV) approach. Two time-dependent variables are defined to change the original governing equations to be the equivalent forms and then the Crank–Nicolson type approximation and the explicit Adams–Bashforth approximation are used to temporally discretize the equivalent equations. All nonlinear parts and auxiliary variables are explicitly treated, thus we can decouple the phase-field variable and auxiliary variables in time. The velocity and pressure are decoupled by using a second-order accurate pressure correction method. Therefore, our numerical scheme is very simple and efficient. We analytically prove the discrete energy dissipation law and unique solvability of our scheme with the absence of external force. The benchmark tests are performed to show that our method has desired accuracy and energy stability. Moreover, our method works well in simulating buoyancy-driven pinchoff and viscous fingering.