Problems of Two-Phase Immiscible Liquids and Methods for Their Solution

This paper provides a literature review on the problems of two-phase immiscible fluids and the methods for their solution. The analysis of such systems has significant practical importance, as it allows for the optimization of separation processes, improvement of heat exchange, and enhancement of device efficiency.

The objective of this literature review is to systematize and analyze existing approaches, theories, and methods used to solve problems related to the dynamics of two-phase immiscible liquids. The review covers various aspects of modeling and solving such problems, including analytical, numerical, and experimental methods. Special attention is given to approaches that account for capillary and interfacial phenomena, as well as surface tension between liquids. This study focuses on various methods, such as the Volume of Fluid (VoF) method, the phase-field method, and the smoothed particle method, along with their applications in modeling complex processes like droplet breakup and coalescence, phase transitions, and the interaction of two liquids under different physical conditions.

Another key objective is to enhance the computational efficiency of these methods, identify their advantages and limitations, and determine the most suitable approaches for practical applications, such as filtration, liquid extraction, and process control in petrochemical and other industries.

4.1 Types of problems for two-phase immiscible liquids

4.1.1 Modeling of two-phase flow systems

Recent studies [1-6] have focused on modeling the flow of two immiscible liquids, such as oil and water, in pipes, reservoirs, or open systems. This problem involves determining the distribution of pressure, velocity, temperature, and phase behavior [7-9].

In investigations examining how submergence ratio affects pump performance, the length of the two-phase flow section was held fixed, while liquid flow rates were recorded across different submergence levels. Conversely, to assess the impact of two-phase section length, the height of the liquid column was kept unchanged, and measurements were carried out for various lengths of the two-phase region (see Figure 1).

1.png

Figure 1. Experimental apparatus of airlift [3]

Because the operation of the airlift pump is inherently unsteady, the liquid flow rate recorded during experiments tends to oscillate around specific levels. To ensure greater consistency and reliability of the data, a time-averaging (TA) approach was employed for the liquid flow rate measurements. The averaged values were determined using the flowing procedure:

$T A=\frac{1}{n} \sum_{i=1}^n E X P_i$                        (1)

The graphs and illustrations are provided with detailed explanations. For example, Figure 2 illustrates the averaged profile of the liquid surface velocity under different flow regimes.

2.png

Figure 2. The averaged profile of the liquid surface velocity under different flow regimes [3]

Under the assumption that the air expands isothermally within the pipe, the energy supplied to the airlift pump corresponds to the work required for air compression. In contrast, the useful output is represented by the gravitational potential energy gained by the lifted liquid. The pump's efficiency is thus characterized by the ratio of the useful output to the energy input, and can be expressed as follows:

$\eta=\frac{\rho g Q_L(L-H)}{P_a Q_G \ln \frac{P_{i n}}{P_a}}$                             (2)

where, $Q_L$ is the volumetric flow rate of the liquid, $Q_G$ is the volumetric flow rate of the gas, $P_{\text {in }}$ is the inlet air pressure, $\rho$ is the liquid density, and $P_a$ is the atmospheric pressure.

The interaction between the liquid and gas surface velocities was analyzed for submergence depths of 2.8, 2.5, and 2.0 meters. The results indicate that operating the pump at greater submergence levels generally leads to improved performance. However, since the points of maximum liquid flow and peak efficiency do not coincide, the choice of gas flow rate during real-world application should be tailored to the operational objective—whether it prioritizes higher efficiency or greater liquid throughput.

The proposed model has a significant error at high gas flow rates. To improve accuracy under such conditions, this study introduces a more precise and simplified calculation formula for the airlift pump, based on a model that accounts for pipeline losses.

The Bernoulli equation for the gas injection interface and the horizontal liquid level is established as follows:

$P_1=P_a+\rho_L g H+\frac{P_L}{2}\left(V_3^2-V_1^2\right)-W_h-W_j$                         (3)

where, $P_1$ and $P_a$ are the section and atmospheric pressures $\left(\mathrm{N} / \mathrm{m}^2\right), \rho_L$ is the liquid density $\left(\mathrm{kg} / \mathrm{m}^3\right), \mathrm{g}$ is the gravitational acceleration $\left(\mathrm{m} / \mathrm{s}^2\right), H$ is the distance between the gas inlet and the horizontal liquid level, $V_1$ and $V_3$ are the liquid flow velocities at respective sections $(\mathrm{m} / \mathrm{s}), W_h$ and $W_j$ represent pipeline and local losses between sections. The indices 1,2 , $3, \mathrm{G}$, and L denote sections 1-1, 2-2, horizontal section 3-3, gas phase, and liquid phase, respectively.

$W_h=\sum_{i=1}^n \frac{l_i f V_i^2 \rho_L}{2 D_i}$                        (4)

$W_j=\sum_{j=1}^n \frac{\varepsilon_j V_j^2 \rho_L}{2}$                         (5)

where, $\varepsilon_1=0.3, \varepsilon_2=0.5, \varepsilon_3=0.45, \varepsilon_4=1, f$ is the pipeline friction coefficient, $l$ is the pipeline length ( m ), D is the pipeline cross-sectional area $\left(\mathrm{m}^2\right), \varepsilon$ is the local loss coefficient, with subscripts i and j indicating the ith and j -th sections of the pipeline, $l_1=3.254+H-L \mathrm{~m}, l_2=$ $0.15 \mathrm{~m}, l_3=0.2 \mathrm{~m}, l_4=0.35 \mathrm{~m}$.

In the injector, gas is injected with a volumetric flow rate. Assuming that the velocity of the mixture leaving the injector is V₂ and neglecting gas density variations, we obtain the volume continuity equation:

$A V_2=Q_G+A V_1=Q_G+Q_L$                       (6)

$V_2=\frac{Q_G}{A}+V_1=V_1\left(1+\frac{Q_G}{Q_L}\right)$                         (7)

where, A is the pipe cross-sectional area (m²).

At the gas injection point, the density of the gas is significantly lower than that of the liquid. As a result, the gas mass flow rate is negligible in comparison to the liquid mass flow rate. This allows the application of the continuity equation in simplified form:

$\rho_2 A V_2=\rho_L A V_1$                           (8)

$\rho_2=\rho_L \frac{V_1}{V_2}=\frac{\rho_L}{\left(1+\frac{Q_G}{Q_L}\right)}$                          (9)

where, $\rho_2$ is the density of the mixture at the gas inlet (kg/m³).

Applying the momentum equation to the injector as a control volume, we obtain:

$P_2=P_1-\rho_L V_1\left(V_2-V_1\right)=P_1-\rho_L V_1 \frac{Q_G}{A}$                          (10)

Thus, combining Eqs. (3) and (10):

$\begin{aligned} P_2=P_a+\rho_L g H & +\frac{\rho_L}{2}\left(V_3^2-V_1^2\right)-W_h-W_j-\rho_L V_1 \frac{Q_G}{A}\end{aligned}$                              (11)

Furthermore, the momentum equation for the two-phase gas-liquid section above the injection port gives:

$P_2=P_a+K \frac{\rho_L V_1^2}{2}\left(1+\frac{Q_G}{Q_L}\right)+\frac{\rho_L g L}{1+\frac{Q_G}{s Q_L}}$                            (12)

$K=\frac{4 f L}{D}$                           (13)

The slip ratio S is given by:

$s=\frac{V_G}{V_L}=1.2+0.2 \frac{Q_G}{Q_L}+\frac{0.35 \sqrt{g D}}{V_1}$                           (14)

Experimental observations reveal that the liquid flow rate rises with an increase in gas injection up to a certain point, after which it tends to level off. Notably, the highest liquid throughput does not always align with the point of maximum pumping efficiency.

4.1.2 Interface dynamics between two liquids

Studies investigate the behavior of the interface between two immiscible liquids, such as when one liquid (e.g., water) drips or bubbles in another (e.g., oil) [10-16]. The modeling of water droplets in oil, their breakup, coalescence, or shape change is also explored [17].

The VoF approach is commonly employed to capture the interface between immiscible phases, whereas the Phase-Field method provides a framework for simulating the evolution of droplets and bubbles. Furthermore, periodic forcing at the phase boundary has been shown to modify cavity behavior, offering potential for advanced fluid manipulation techniques [18]. For mesh-free simulations involving droplets and bubbles in liquid media, the Smoothed Particle Hydrodynamics (SPH) method is widely utilized.

The mathematical model relies on the Navier–Stokes framework, formulated for incompressible, immiscible two-phase flows. The governing equations include the continuity equation and the momentum conservation equation, expressed as follows:

$\nabla \cdot u=0$                        (15)

$\begin{gathered}\frac{\partial(\rho u)}{\partial t}+\nabla \cdot(\rho u \otimes u)=-\nabla p+\nabla \cdot\left[\mu\left(\nabla u+\nabla u^T\right)\right]+\rho g+f\end{gathered}$                             (16)

where, $u$ is the velocity vector, $p$ is pressure, $g$ is the gravitational force vector, $f$ is the external force vector, and $t$ is time. The fluid parameters $\rho$ and $\mu$ represent the density and dynamic viscosity, respectively. From a physical property's perspective, the flow solver includes relationships for density and dynamic viscosity, given as:

$\rho=F \rho^a+(1-F) \rho^b$                          (17)

$\mu=F \mu^a+(1-F) \mu^b$                         (18)

where, the two fluid phases are denoted as a and b, with a representing the volume fraction of phase a within each computational cell.

Several numerical schemes are available for discretizing the governing equations, including the finite difference (FDM), finite volume (FVM), and finite element (FEM) methods. In this study, the computations were performed using a flow solver operating on a Cartesian grid.

4.1.3 Surface tension model

In cases where interfacial forces play a crucial role, surface tension cannot be neglected. The Continuum Surface Force (CSF) model is incorporated into the Navier-Stokes equations to account for surface tension effects. This force can be expressed as:

$f_{C S F}=\sigma \kappa n \delta$                          (19)

where, $\sigma$ is the surface tension coefficient, $n$ is the interface normal, $\kappa=-\nabla \cdot n$ is the interface curvature and $\delta$ is the Dirac delta function. We use the approximations $\delta=|\nabla F|$ and $n=\frac{\nabla F}{|\nabla F|}$ to model surface tension as a function of the volume fraction F.

Recent studies on the method of moment determination from liquid focus on improving computational efficiency and expanding their relevance for real-world applications. Additionally, they aim to address issues related to accuracy, reliability, and efficiency (Figure 3).

3.png

Figure 3. Algorithm for solving the Navier-Stokes equations

4.1.4 Modeling of capillary and surface phenomena

Capillary effects are fundamental to the behavior of two-phase systems, particularly in processes such as liquid extraction and filtration. Prior studies [15, 19, 20] have examined phenomena like capillary rise in narrow tubes and the interaction of droplets at fluid interfaces.

To incorporate capillary forces into numerical models, phase-field approaches are employed, which include formulations for surface tension and Van der Waals interactions. In the case of a liquid film of thickness HHH confined between two parallel surfaces, the Van der Waals attraction manifests either as a surface-normal stress or as a disjoining pressure distributed within the film.

When transitioning to a Lagrangian description, the governing equations for weakly compressible, viscous flows are expressed in the following form [21, 22]:

$\frac{d p}{d t}=-\rho \nabla \cdot \mathrm{u}$                          (20)

$\rho \frac{d u}{d t}=\nabla \cdot \sigma+F_s+F_{\vartheta}$                         (21)

$\frac{d x}{d t}=u$                           (22)

where, $\rho$ is the density, $t$ is time, $u$ is velocity, $\sigma$ is the total stress tensor, $r$ is position. $F_s, F_{\vartheta}$ represent surface tension forces and Van der Waals forces per unit volume, respectively. The symbol $\frac{d}{d t}$ denotes the material derivative, $\frac{d}{d t}=\frac{\partial}{\partial t}+u \cdot \nabla$.

In experimental studies and theoretical analysis [23, 24], real-time tracking of physical fields and intuitive analysis of morphological evolution remain challenging, and the investigation of droplet collision dynamics is still ongoing.

4.1.5 Modeling of phase separation and state change

These problems [25-27] involve analyzing the phase transition of one of the liquids, such as evaporation, condensation, or freezing, which can affect the behavior of a two-phase system. This includes modeling the separation of oil and water when one of the phases evaporates or condenses due to changes in temperature or pressure.

4.1.6 Modeling of two-phase flows in porous media

The study [28] presents modeling of two-phase flows in porous media, which is used, for example, to evaluate oil displacement by water in oil reservoirs or the filtration of liquids through porous membranes.

Modeling oil displacement by water in an oil reservoir requires considering capillary and gravitational forces, as well as interactions with the porous material. Due to the complexity and nonlinear nature of the Stokes-Darcy equations, obtaining an analytical solution is often impossible.

Darcy's Law describes the flow of fluid through a porous medium and is expressed as:

$q=-\frac{k}{\mu} \nabla p$                          (23)

where, $q$ is the fluid flux density, $k$ is the permeability coefficient, $\mu$ is the dynamic viscosity of the fluid, $\nabla p$ is the pressure gradient.

Darcy's Law is applicable to laminar fluid flow through a porous medium [29, 30]. For high-velocity or turbulent flows, more complex models may be required, such as extended Darcy equations accounting for turbulence or two-phase flow equations.

In cases involving complex geometries, multicomponent fluids, or anisotropic materials, modifications of Darcy's Law must be applied, or more advanced approaches should be used, such as solving the Navier-Stokes equations in porous media.

4.2 Methods for solving problems of two-phase immiscible liquids

4.2.1 VoF method

The VoF method is used to track and compute the boundaries between two fluids. This method assumes that each grid cell stores the fraction of each fluid, allowing for the tracking of phase boundaries and the changing shape of droplets. It is well-suited for modeling fluid flows and interactions between liquids [16, 31]. A computational grid is used to solve the Navier-Stokes equations for each fluid. At each time step, the interface position between the fluids is calculated, and the grid is adapted accordingly.

4.2.2 FEM

In the following studies [32-35], the FEM is used for the numerical solution of equations describing fluid behavior, such as equations for multi-fluid systems. This method allows working with unstructured grids, which is useful for modeling complex geometries.

Features: It allows solving equations for each phase separately, as well as their interaction.

It is effective for solving problems involving the mechanical behavior of fluids (e.g., in droplet and bubble dynamics). Each fluid is modeled through its own system of equations, and then interactions between the fluids are determined through phase boundary conditions.

4.2.3 Phase-field method

The phase-field method [36] is used to model the dynamics of the interface between two phases, where the transition between fluids occurs through a continuous variation of a field. It allows tracking changes in the shape and size of liquid droplets or bubbles, as well as their merging or breakup.

Diffuse crack model in the phase-field method.

$\psi_{p o t}(u, \Gamma)=\psi_e+\psi_p+\psi_{f r a c}-\psi_{e x t}$                         (24)

The first part of the total potential energy is the elastic energy $\psi_e$, which can be decomposed into tensile and compressive components. In this implementation, only the tensile component can crack, and the damaged material degrades using the degradation function $g(\phi)$, which is defined as follows:

$g(\phi)=(1-\phi)^2+k$                        (25)

where, k is a small number.

The linear strain tensor $\varepsilon=\varepsilon(u)$ is formulated as:

$\varepsilon_{i j}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$                          (26)

Under the assumption of isotropic linear elasticity, the linear energy density $\psi_e(\varepsilon)$ is expressed in terms of the constants λ and μ as:

$\psi_e=\frac{1}{2} \lambda \varepsilon_{i i} \varepsilon_{j j}+\mu \varepsilon_{i j} \varepsilon_{i j}$                         (27)

Then, the degradation function, which contains the phase-field value, reduces the tensile strain energy density. Since the strain energy consists of tensile strain (+) and compressive strain (–), the elastic energy is expressed as:

$\begin{aligned} \psi_e(\varepsilon)=\int_{\Omega} \psi_{\varepsilon} d \Omega & =\int_{\Omega}\left[(1-\phi)^2+k\right] \cdot \psi_e(\varepsilon)^{+}+\psi_e(\varepsilon)^{-} d \Omega\end{aligned}$                         (28)

The plastic energy density $\psi_p$ for linear isotropic hardening is formulated as:

$\psi_p(\delta, d)=\frac{1}{2} \widetilde{H}(d) \delta^2$                          (29)

where, $\widetilde{H}(d)$ is the hardening modulus of the degraded material, and $\delta$ is the accumulated plastic strain. It should be noted that both compressive and tensile plastic deformation equally contribute to crack propagation, which differs from the state of elastic deformation. The plastic deformation energy can be expressed as:

$\left.\psi_p(\delta)=\int_{\Omega}\left[(1-\phi)^2+k\right]\left(\frac{1}{2} \tilde{H}(\phi) \delta^2+\widetilde{\sigma_y} \delta\right)\right] d \Omega$                         (30)

The surface energy density function in phase-field theory is written as follows:

$\gamma(\phi, \nabla \phi)=\frac{1}{2 l_0} \cdot \phi^2+\frac{l_0}{2} \cdot\left|\nabla_\phi\right|^2$                             (31)

Using the crack surface density function (8), the considered fracture energy can be regarded as the minimum energy required to create a new crack region and can then be expressed with the critical energy release rate G_c as:

$\begin{aligned} \int_{\Gamma} G_c d \Gamma \approx \int_{\Omega} G_c \cdot & \gamma\left(l_0, \phi, \nabla \phi\right) d \Omega =\int_{\Omega} G_c \cdot\left[\frac{1}{2 l_0} \cdot \phi^2+\frac{l_0}{2} \cdot\left|\nabla_\phi\right|^2\right] d \Omega\end{aligned}$                         (32)

This function allows modeling complex interactions between fluids and phase boundaries.

The phase-field method was chosen for its accurate modeling and computational efficiency, reducing the typically high computational costs associated with detailed modeling without losing accuracy. It is effective for the dynamics of droplets, bubbles, and other interfacial phenomena.

4.2.4 SPH method

In studies [37, 38], this method utilizes particles representing the fluid to track its movement in space. The particles are not confined to a grid and can adapt to changes in the shape and position of the interface. It has an advantage over traditional methods (e.g., VoF) as it does not require a grid, making it more suitable for modeling systems with dynamic interfaces. It can be used for modeling droplets and bubbles.

Each particle represents a fluid element that interacts with neighboring particles. The equations of motion and interactions are solved for each particle.

4.2.5 Methods considering capillary forces and surface tension

The studied methods [39-41] are used to model phenomena related to surface tension between two immiscible liquids, such as emulsification, capillary rise, or bubble formation.

In this sensitivity study, two dimensionless numbers are introduced to determine which phenomenon dominates under various surface tension conditions. The first is the Bond number Bo(–), which represents the ratio of gravitational forces to surface tension forces:

$B o=\frac{\rho \underset{g}\rightarrow l^2  }{\sigma}$                         (33)

where, ρ (kg/m³) is the liquid density, g (m/s²) is the gravitational constant, σ (kg/s²) is the surface tension, l (m) is the characteristic length of the molten bath. The second introduced number is the Nusselt number Nu(–), which is defined as:

$N u=\frac{h l}{k_f}$                         (34)

where, $h\left(W /\left(\mathrm{m}^2 \cdot K\right)\right)$ is the convective heat transfer coefficient of the molten bath, $k_f(W /(\mathrm{m} \cdot K))$ is the thermal conductivity of the liquid, $l(\mathrm{~m})$ is the characteristic length, in this case, the length of the molten bath.

The effective thermal conductivity is determined using a high-fidelity numerical simulation of thermo-hydrodynamic behavior. This refined property can then serve as input for a reduced-order model, allowing it to account for neglected multiphysics effects and improve the accuracy of predicted temperature distributions. As a result, reduced-order models offer a reliable and time-efficient alternative for thermal analysis with significantly lower computational costs.

This study presents mesoscale numerical modeling of thermofluid dynamics using specialized software. The simulated shape and size of the molten pool show strong agreement with experimental observations. Building on this validation, two parametric studies are carried out to investigate the influence of recoil pressure and capillary forces under molten pool conditions.

The outcomes of the two parametric studies help define a process–material parameter space, highlighting which physical mechanisms significantly affect the system and which can be considered negligible. This insight enables informed decisions about when a detailed multiphysics model is necessary and when a simplified, computationally efficient reduced-order model can be reliably used without substantial loss of accuracy.

4.2.6 Description of adaptive finite volume method

The adaptive grid strategy is based on error estimation of interface curvature, where refinement is triggered when local curvature exceeds a predefined threshold. Grid cells near the interface are refined using a quadtree/octree structure, maintaining volume conservation. The capillary force is coupled through the CSF model, with surface tension forces applied at reconstructed interfaces using smoothed volume fraction gradients [41-44].