An Accurate Phase Interface Locating Algorithm for Pore-Scale Two-Phase Interfacial Flows

Mei, Qiliang; Wu, Tianjiang; Li, Chao; Zhou, Xiaoying

doi:https://doi.org/10.1155/2023/2703312

Geofluids

On this page

Abstract Introduction Conclusion Data Availability Conflicts of Interest Authors’ Contributions Acknowledgments References Copyright Related Articles

Special Issue

Pore-Scale Fluid Flow and Controlling Mechanisms in Unconventional Oil/Gas Reservoirs

View this Special Issue

Research Article | Open Access

Volume 2023 | Article ID 2703312 | https://doi.org/10.1155/2023/2703312

An Accurate Phase Interface Locating Algorithm for Pore-Scale Two-Phase Interfacial Flows

Qiliang Mei,^1,2Tianjiang Wu ,³Chao Li,^1,2and Xiaoying Zhou^1,2

Academic Editor: Dengke Liu

Received22 May 2022

Accepted29 Jul 2022

Published10 Apr 2023

Abstract

Direct numerical simulation of pore-scale flow in porous media at pore scale is a fast developing technique to investigate pore-scale flow behaviors. However, with the decrease of the spatial scale, the interfacial tension of the two-phase interface will have an important impact on the two-phase flow processes. As a result, spurious currents near the interface with this method have an important impact on the numerical accuracy and computational efficiency, seriously affecting the numerical prediction of pore-scale fluid flow. Ghost cell method can greatly reduce the spurious currents near the interface, but it needs to locate the phase interface accurately. In this work, a new method using constant velocity spiral to approximate the two-phase interface is proposed. The method not only considers the curvature of the phase interface but also considers the influence of the curvature change on the phase interface, which greatly improves the ability of the phase interface position. The new method can improve the accuracy of interfacial tension treatment and then improve the prediction ability of volume of fluid in high interfacial tension driven flow. Numerical simulations of capillary rising, droplet spreading on a plane, and bubble rising show that the method is accurate and has a strong engineering application prospect.

1. Introduction

The two-phase flow process with a phase interface exists widely in nature (such as rivers and lakes) and industrial processes (oil recovery [1]).To deeply study the dynamic behaviors and mechanism of the fluid behaviors at pore scale is of great significance to understand macroscopic behaviors and optimize development scheme. The research methods for the flow with the phase interface mainly include experimental methods [2–4] and numerical simulation methods [5–7]. The experimental methods can provide the macroscopic characteristics and phenomena of phase flow but do little for the research on the driving mechanism of these macroscopic characteristics and phenomena. With the development of related numerical algorithms and improvement of computing ability of the computers, numerical simulation method has gradually become an important tool to study the multiphase flow with the phase interface.

In the numerical simulation of pore-scale two-phase or multiphase flow, Navier-Stokes equation is generally used [7]. For instance, Raeini et al. used Navier-Stokes equation to simulate pore-scale two-phase flow [8]. Guo et al. used Navier-Stokes method to investigate the residual oil distribution and enhanced oil recovery method [9]. Ning et al.’s group used Navier-Stokes Equation in the pore-scale simulation of oil-water flow [10], pore-scale particle transport process [11, 12], or particle flooding process [1] and have gained interesting results of pore-scale multiphase behaviors (for instance, capillary barriers [10, 13]). In the oil-water two-phase flow, there is always interfacial tension which depends on the curvature of interface and the interfacial energy per unit area. Owing to a sharp increase of curvature in the microscale condition, interfacial tension has to be the vital one among kinetic behaviors. Therefore, only accurate computation of interfacial tension can guarantee accurate prediction of oil-water flow behaviors. From numerical point of view, interfacial tension only exists on interface. It has singularity which makes it harder when coupling with Naiver-Stokes equation. In order to compute it more accurately, tracking of the position matters a lot. For now, there are serval ways for tracking interface, front tracking method (FTM) [14], level set method [15], and volume of fluid (VOF) [16]. FTM uses two-dimensional curved surface mesh to represent the phase interface, which makes the position confirmed. However, when there is a change in the surface topology, the mesh needs to be changed accordingly. In extreme cases, we have to reconstruct the mesh, making the treatment intractable. Level set method can track topology accurately and is easier extended to three-dimensional cases, but it needs to reinitiate level set function due to restricted precision and may cause mass increase or loss. VOF uses distribution of volume fraction of the primary phase to represent the phase and interface distribution and pure convection scalar equation to track its change. As a result, the mass is conserved. Moreover, it is easier to deal with the changes of topology and to extend to three-dimensional cases. This method has been extensively used in pore-scale multiphase flow [10, 17]. With this method, in order to trace movement of interface accurately and reduce diffusion, high accurate algorithms and subgrid interface structure is required.

For conventional VOF, nonphysical flow near the interface, also known as “spurious currents,” usually appears under high surface tension [18]. Generally speaking, this problem is caused by two reasons. First, it is difficult to accurately compute the position and curvature of interface. Second, surface tension-induced pressure and surface tension are unbalanced during Navier-Stokes discretization process. The first problem has always been a hot spot. Solutions such as finite difference method [19], ELVIRA method [20], and least square method [21] are employed to improve the accuracy of interface normal calculation, whereas CLSVOF (coupled level set with volume of fluid) method [22, 23] and linear interface reconstruction method [24–26] are used to obtain position, and the method in reference [18] can increase the precision to compute curvature. As for the second problem, reference [27] gives a strategy to make pressure and interfacial tension balanced by introducing interfacial tension-induced pressure into pressure equation rather than dealing with interfacial tension directly in the momentum equation. Moreover, they used ghost cell method to sharply represent the interface. If the curvature is computed accurately, interfacial tension calculation can achieve the accuracy of the computer. In the ghost cell method for interfacial tension treatment, the intersection point of the connection line of interfacial cells with the interface needs to be computed. The interfacial cells here is denoted as two neighbor cells across the phase interface. The accuracy of computing the intersection point matters a lot for interfacial tension computation with high accuracy, which relates with the dynamic behavior prediction in interfacial tension driven flows. Presently, interface distance interpolation is usually employed to compute this point, which does not consider the influence of interface curvature. Sun and Su used an arc to approximate the interface near in order to enhance the accuracy of solving for considering the influence of interface curvature on the position of intersection point [28]. However, the actual interface is not an arc, and the curvature of the interface may vary with the position.

This article presents a new method for computing the intersection point through approximating the subgrid interface using a variable curvature curve. The influence of interface curvature on the intersection point is taken into account. As a result, the intersection point and the interface curvature at this point are obtained with high accuracy.

2. Mathematical Model

The Navier-Stokes equation is used to describe the two-phase motion.

In addition, the two-phase fluid meets the continuity restriction.

where and are the average density and dynamic viscosity coefficient of two phases, calculated from

where and are the properties of phase 1 and phase 2 and is the volume fraction of phase 1. Phase 1 is the primary phase, and phase 2 is the secondary phase; is dynamic pressure; is gravity acceleration; is the average velocity of two phases. The following equation is used to trace .

Equations (1) and (2) are used to describe the motion of two phases.

τ is deformation rate tensor, given by

When the two-phase interface exists, the pressure will have a jump; i.e., the following equation is met: where is the pressure at the interface in the owner side and is the pressure at the interface on the neighbor side. Owner and neighbour are two cells across the interface and sharing the same face. is the interfacial tension coefficient, and is the curvature of the interface, which is calculated from the following equation:

3. Solution Method and Procedure

In this paper, the collocated finite volume method was used to discretize the equations, and the PISO algorithm was used to decouple the pressure and velocity [29]. The pressure difference caused by the interfacial tension is separated from the actual pressure. The half-discrete form of Equation (1) is as follows: where is the value on the diagonal of the coefficient matrix after discretization of the Equation (1) and is expressed as where is the implicit contribution coefficient of neighbor nodes to the current node during discretization of Equation (1) and is the explicit discrete contribution except pressure; p_c is the pressure gradient caused by the presence of interfacial tension, and is the dynamic pressure.

Both sides of Equation (6) are divided by and substituted in Equation (2) to give

Equation (13) is the pressure equation derived by PISO algorithm. The new pressure can be obtained by solving the equation and substituted into Equation (6) to update the velocity. The velocity flux at the cell face is calculated from

It should be pointed out that the pressure gradient, , caused by interfacial tension in Equations (9) and (10) should be further discussed.

The arbitrary polyhedral finite volume method is used to discretize the equation, and the gamma scheme presented in work [30] was used to discretize the convection term, and the Crank-Nicolson scheme was for the time term [31], and the central scheme is used to discretize the diffusion term. The solution procedure of the equations is as follows.

Step 1. Solving velocity equation (Equation (1)) for velocity prediction. Pressure at the previous time step is used.

Step 2. Solving pressure Equation (10) to obtain the pressure.

Step 3. Using Equation (11) to update velocity flux at the cell face.

Step 4. Repeating Steps 2 and 3 to complete the pressure calculation cycle.

Step 5. Calculating the velocity at the cell center using Equation (8).

Step 6. Calculating the volume fraction using Equation (4).

Step 7. Updating the density and viscosity using Equation (3).

Step 8. Going to Step 1 for next time cycle.

4. Interfacial Tension Processing

4.1. Ghost Cell Method

In this paper, the ghost cell method was used to deal with the interfacial tension [27]. Figure 1 illustrates the two-phase interface across two nearby mesh cells. The dotted line represents the position of the interface; and are the centers of the two adjacent cells, and is the position vector of the interface between the two elements. is the intersection point between the line and the phase interface, and is the intersection point between the line and the cell face. It can be seen from Equation (5) that the interfacial tension leads to higher pressure at the primary phase side than the secondary phase side. Therefore, point or in the primary phase results in a difference in the direction of the pressure difference at the interface. If , is in the primary phase. If , is in the primary phase.

When is in the primary phase, the pressure distribution along is shown in Figure 2(a); when is in the primary phase, the pressure distribution along is shown in Figure 2(b).

(a)

(b)

To solve Equations (10) and (11), it is necessary to calculate the pressure gradient (including the pressure gradient, , and dynamic pressure gradient, , caused by interfacial tension) at the interface . When is in the primary phase, the pressure gradient at the interface can be expressed as

When is in the primary phase, it can be expressed as

where is the curvature of the interface; it is calculated from the following equation:

Equations (12) and (13) can be combined to give

where is a sign function. If , ; if , . The first and second terms on the right of Equation (15) are the value of dynamic pressure, , on the cell face and the value, , of the pressure gradient, p_c, on the cell face caused by the interfacial tension.

The center of the element used in Equation (8) is reconstructed by the following equation: where includes all the faces of a cell and is the area vector of the face of the cell concerned and points to the outside of the cell.

4.2. Spiral-Based Interface Location Algorithm

In Equation (14), is the intersection point between and the phase interface. The accurate location of this intersection point is very important to enhance the numerical accuracy of interfacial tension evaluation and further reduce the nonphysical spurious currents near the interface. Assuming that the interface passing through the point is a plane, the position of the point is usually obtained by linear interpolation of the distance from to the interface. Sun and Su used an arc to approximate the interface near in order to enhance the accuracy of solving for considering the influence of interface curvature on the position of intersection point [28]. However, the actual interface is not an arc, and the curvature of the interface may vary with the position. Thus, based on the trajectory equation of constant velocity spiral, the influence of the interface curvature variation on the interface position was considered, and the curvature at the intersection point between the interface and the central line of the mesh cell is considered in this paper.

When the curvature at points and is small, the interface can be approximated to a plane. The intersection point is calculated by linear interpolation. where

where and are the distances from the centers of elements to the interface. and are the locations of the centers of cell and .

In the interface tracking process, the fluid interface between and is usually not a plane, and the interfaces at different positions have different curvature. Therefore, the variable curvature interface is approximated by the constant velocity spiral as follows. (1)The center, , of the constant velocity spiral is found first

Figure 3 shows the approximation of the constant velocity spiral of the curved surface. The grey-curved surface indicates the fluid phase interface; and are the nearest points of and on the phase interface, respectively; and is the intersection point of the phase interface and the straight line . The center of the constant velocity spiral is calculated from the following formula:

where and are the coordinates of the centers of elements and , and are the curvature radius of elements and , and and are the normal directions of the faces calculated by positions and . (2)With the center , the constant velocity spiral is constructed, and the intersection point of constant velocity spiral and is calculated

The coordinates of projection points, and , can be expressed as

If , , and are on a straight line, the curvature radius on the interface can be calculated by the following formula:

If , , and are not on a straight line, the following local coordinate system is constructed: where , , and are the directions of the three axes in the local coordinate system, where the constant velocity spiral is constructed in plane . The projections of and on plane are as follows:

The constant velocity spiral is constructed by and . If the intersection point of the interface and is approximated by the intersection point of the constant velocity spiral and , the following relationship is satisfied: where is the angle between and : where is the angle between CI and . Point is on and can be expressed as

During the calculation, is the angle () between and . The counterclockwise is positive, while the clockwise is negative. The inner product of and is used to determine positive or negative .

Equations (29) and (26) are substituted into Equation (24). The bisection method [32] is used to obtain , i.e. the position of the intersection point of the interface and . (3)Calculation of the curvature at the intersection point

Once is obtained, the curvature radius at the intersection point is calculated from the following formula:

Thus,

5. Test Case

5.1. Capillary Rising

Driven by the interfacial tension, the water rises or falls along the capillary wall and finally stays at a certain height when wetting the capillary. The retention height is related to the wettability, interfacial tension coefficient, and capillary radius as follows:

where is the contact angle between water and capillary wall, is the interfacial tension coefficient, is the density of water, is the acceleration of gravity, and is the capillary radius.

In this paper, the accuracy of the algorithm for simulation of the flow process driven by the interfacial tension was verified by the numerical simulation of the capillary rising process. Table 1 lists the physical parameters used in this case. The capillary radius , and the water height in the capillary is 0.008 m at the initial state. Under the joint action of gravity and interfacial tension, the water-air interface will move and finally stay at a fixed position. In order to study the influence of wettability on the capillary rising height, , , , , and are calculated in this case, and the liquid rising height at different wettability is obtained.

Figure 4 shows the relationship between the final liquid height and wettability. Moreover, the calculations given by Formula (31) are also shown in Figure 4 as a reference. It can be seen from Figure 4 that the final equilibrium height gradually decreases with the decrease of contact angle. The numerical simulation results given in this paper are in good agreement with the analysis results.

5.2. Droplet Spreading

When a droplet drops on a smooth horizontal wall, the droplet spreads on the horizontal surface and forms a stable shape under its wetting action. The spreading process of the droplet is mainly affected by gravity, interfacial tension, and viscous force, and the final stable shape is directly related to the radius and contact angle of the initial droplet.

The initial state of the semicircular droplet is shown in Figure 5(a). The droplet has the radius . Under the joint action of gravity, viscous force, and interfacial tension, the final shape of the droplet is shown in Figure 5(b). is the height of the droplet, is the width after the droplet spreads, and is the contact angle. Considering no gravity effect and assuming that the droplet is arc-shaped, the spreading width, , and height, , are calculated as follows according to the mass conservation of the droplet.

(a)

(b)

In order to compare the numerical simulation results with the analysis results, the gravity effect is ignored in the numerical simulation. The physical parameters of the fluid used are shown in Table 2. Figure 6 shows the shapes of the stable droplet at different contact angles: (a) , (b) , (c) , (d) , and (e) . The final shape of the droplet at different contact angles can be qualitatively seen from Figure 6. The spreading width and height of droplets at different contact angles are measured and compared with the results obtained by Equations (33) and (34), as shown in Figure 7. It can be seen from the figure that with the increase of contact angle, the height of the droplet increases, but the width of the droplet decreases. The numerical simulation results are in good agreement with the analysis results.

5.3. Bubble Rising

Bubble rising is an important multiphase flow phenomenon with a free interface. The bubble shape, which is controlled by viscous force, interfacial tension, and buoyancy, changes with time. The accuracy of interfacial tension processing plays an important role in prediction of bubble shape.

The computational domain and initial position of the bubble are shown in Figure 8. The computational domain is 2-D, 0.25 m wide and 0.3 m tall. The initial bubble has its center 0.05 m from the bottom and is round with a radius of 0.025 m. There are 250 cells in width and 300 cells in height. The nonslip boundary condition is adopted for the walls. The physical parameters used in the simulation are listed in Table 3 and consistent with reference [33].

In order to verify the feasibility of the proposed algorithm, the experimental results in reference [34] are shown in Figure 9. It can be seen from Figure 9 that the initial bubble is round (as shown in Figure 9(a)) and rises under the action of buoyancy (as shown in Figure 9(b)); the center of bubble rises very fast, and both sides of the bubble rises very slowly, so the rising bubble gradually forms a “yoke” (as shown in Figures 9(c) and 9(d)); as the bubble rises, the bubble further deforms; i.e. the upper part of the bubble forms a crescent, and the gas on both sides gradually separates from the main bubble, as shown in Figure 9(e); as the bubble further rises, smaller bubbles formed on both sides separate from the main bubble, as shown in Figure 9(f). The bubble shape predicted by numerical simulation at different time is in good agreement with the experimental data, meaning that the algorithm is accurate and feasible.

In order to further compare the numerical simulation results with the experimental results quantitatively, the position of the bubble and thickness of main bubble were measured at different time and compared with the experimental results in this paper. Figure 10 shows the change of dimensionless displacement, ( is the displacement of the bubble, and is the radius of the initial bubble), of the bubbler with dimensionless time, ( is physical time, and is acceleration of gravity). It can be seen from Figure 10 that the vertical position of the bubble is approximately linear with time. The simulation results in this paper are closer to the results of experiment 2. Due to the limitation of measurement methods, the experimental results fluctuate slightly at the initial time. Figure 11 shows the change of dimensionless bubble thickness, ( is the bubble thickness), with dimensionless time, . It can be seen from Figure 11 that the simulation results of bubble thickness in this paper are in good agreement with the experimental results.

6. Conclusion

Under the framework of volume of fluid for simulation of multiphase flow with a phase interface, a high-precision locating algorithm for two-phase interface based on constant velocity spiral was proposed to solve the problem that the ghost cell method for dealing with the interfacial tension needs accurate location of the phase interface. During location of the phase interface by the algorithm, both the influence of the interface curvature on the location of the interface and the influence of the change of the curvature are taken into account to greatly enhance the location accuracy of the phase interface. The prediction accuracy of the volume of fluid for multiphase flow at high interfacial tension is enhanced by accurately locating the phase interface. Finally, the capillary rising, droplet spreading on a plane, and bubble rising are simulated by using the developed numerical model. The numerical simulation results are consistent with the analysis results or experimental results. This method has a high accuracy and strong engineering application prospect.

Data Availability

The experimental data used to support the findings of this study are included within the manuscript and the supplementary materials.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Authors’ Contributions

Q. Mei and T.J. Wu proposed the main framework of the paper. Q. Mei mainly wrote the paper. Q. Mei, T. Wu, C. Li, and X. Zhou analyze the data. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

This research was funded by the Key R&D Program of Shaanxi Province (No. 2021GXLH-Z-071).

References

J. Su, G. Chai, L. Wang et al., “Direct numerical simulation of pore scale particle-water-oil transport in porous media,” Journal of Petroleum Science and Engineering, vol. 180, no. 9, pp. 159–175, 2019.
View at: Publisher Site | Google Scholar
D. Ren, H. Zhang, Z. Wang, B. Ge, D. Liu, and R. Zhang, “Experimental study on microscale simulation of oil accumulation in sandstone reservoir,” Frontiers in Physics, vol. 10, no. 10, p. 841989, 2022.
View at: Publisher Site | Google Scholar
Z. Li, Z. Gu, R. Li et al., “Investigation on droplet dynamic snap-off process in a short, abrupt constriction,” Chemical Engineering Science, vol. 235, p. 116496, 2021.
View at: Publisher Site | Google Scholar
Z. Li, Z. Gu, R. Li et al., “A geometrical criterion for the dynamic snap-off event of a non-wetting droplet in a rectangular pore–throat microchannel,” Physics of Fluids, vol. 34, no. 4, article 042014, 2022.
View at: Publisher Site | Google Scholar
S. Junwei and X. X. Y. Gu Zhaolin, “Discrete element simulation of particle flow in arbitrarily complex geometries,” Chemical Engineering Science, vol. 66, no. 23, pp. 6069–6088, 2011.
View at: Google Scholar
D. Ren, L. Ma, D. Liu, J. Tao, X. Liu, and R. Zhang, “Control mechanism and parameter simulation of oil-water properties on spontaneous imbibition efficiency of tight sandstone reservoir,” Frontiers in Physics, vol. 10, no. 10, p. 829763, 2022.
View at: Publisher Site | Google Scholar
S. Junwei, W. Le, G. Zhaolin, Y. Zhang, and C. Chungang, “Advances in pore-scale simulation of oil reservoirs,” Energies, vol. 11, no. 5, p. 1132, 2018.
View at: Google Scholar
A. Q. Raeini, M. J. Blunt, and B. Bijeljic, “Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method,” Journal of Computational Physics, vol. 231, no. 17, pp. 5653–5668, 2012.
View at: Publisher Site | Google Scholar
Y. Guo, L. Zhang, G. Zhu et al., “A pore-scale investigation of residual oil distributions and enhanced oil recovery methods,” Energies, vol. 12, no. 19, p. 3732, 2019.
View at: Publisher Site | Google Scholar
T. Ning, M. Xi, B. Hu, L. Wang, C. Huang, and J. Su, “Effect of viscosity action and capillarity on pore-scale oil–water flowing behaviors in a low-permeability sandstone waterflood,” Energies, vol. 14, no. 24, p. 8200, 2021.
View at: Publisher Site | Google Scholar
S. Junwei, C. Guoliang, L. Wang et al., “Pore-scale direct numerical simulation of particle transport in porous media,” Chemical Engineering Science, vol. 199, pp. 613–627, 2019.
View at: Publisher Site | Google Scholar
S. Junwei, C. Guoliang, L. Wang et al., “Direct numerical simulation of particle pore-scale transport through three-dimensional porous media with arbitrarily polyhedral mesh,” Powder Technology, vol. 367, pp. 576–596, 2020.
View at: Publisher Site | Google Scholar
B. Hu, G. Zhaolin, Z. Chenxing, L. Wang, H. Chuanqing, and S. Junwei, “Investigation of the effect of capillary barrier on water–oil movement in water flooding,” Applied Sciences, vol. 12, no. 12, p. 6285, 2022.
View at: Publisher Site | Google Scholar
S. O. Unverdi and G. Tryggvason, “A front-tracking method for viscous, incompressible, multi-fluid flows,” Journal of Computational Physics, vol. 100, no. 1, pp. 25–37, 1992.
View at: Publisher Site | Google Scholar
M. Sussman, P. Smereka, and S. Osher, “A level set approach for computing solutions to incompressible two-phase flow,” Journal of Computational Physics, vol. 114, no. 1, pp. 146–159, 1994.
View at: Publisher Site | Google Scholar
C. W. Hirt and B. D. Nichols, “Volume of fluid (VOF) method for the dynamics of free boundaries,” Journal of Computational Physics, vol. 39, no. 1, pp. 201–225, 1981.
View at: Publisher Site | Google Scholar
X. Liu, L. Wang, J. Wang, and J. Su, “Pore-scale simulation of particle flooding for enhancing oil recovery,” Energies, vol. 14, no. 8, p. 2305, 2021.
View at: Publisher Site | Google Scholar
M. Coquerelle and S. Glockner, “A fourth-order accurate curvature computation in a level set framework for two-phase flows subjected to surface tension forces,” Journal of Computational Physics, vol. 305, pp. 838–876, 2016.
View at: Publisher Site | Google Scholar
D. L. Youngs, “An interface tracking method for a 3D Eulerian hydrodynamics code. Atomic weapons research establishment (AWRE) technical,” Report, vol. 35, 1984.
View at: Google Scholar
J. E. Pilliod and E. G. Puckett, “Second-order accurate volume-of-fluid algorithms for tracking material interfaces,” Journal of Computational Physics, vol. 199, no. 2, pp. 465–502, 2004.
View at: Publisher Site | Google Scholar
R. Scardovelli and S. Zaleski, “Interface reconstruction with least-square fit and split Eulerian–Lagrangian advection,” International Journal for Numerical Methods in Fluids, vol. 41, no. 3, pp. 251–274, 2003.
View at: Publisher Site | Google Scholar
M. Sussman and E. G. Puckett, “A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows,” Journal of Computational Physics, vol. 162, no. 2, pp. 301–337, 2000.
View at: Publisher Site | Google Scholar
N. Balcázar, O. Lehmkuhl, L. Jofre, J. Rigola, and A. Oliva, “A coupled volume-of-fluid/level-set method for simulation of two-phase flows on unstructured meshes,” Computers & Fluids, vol. 124, pp. 12–29, 2016.
View at: Publisher Site | Google Scholar
R. Scardovelli and S. Zaleski, “Direct numerical simulation of free-surface and interfacial flow,” Annual review of fluid mechanics, vol. 31, no. 1, pp. 567–603, 1999.
View at: Publisher Site | Google Scholar
S. Diot and M. M. François, “An interface reconstruction method based on an analytical formula for 3D arbitrary convex cells,” Journal of Computational Physics, vol. 305, pp. 63–74, 2016.
View at: Publisher Site | Google Scholar
J. López, J. Hernández, P. Gómez, and F. Faura, “A new volume conservation enforcement method for PLIC reconstruction in general convex grids,” Journal of Computational Physics, vol. 316, pp. 338–359, 2016.
View at: Publisher Site | Google Scholar
M. M. Francois, S. J. Cummins, E. D. Dendy, D. B. Kothe, J. M. Sicilian, and M. W. Williams, “A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework,” Journal of Computational Physics, vol. 213, no. 1, pp. 141–173, 2006.
View at: Publisher Site | Google Scholar
H. X. Sun and J. W. Su, “A sub-grid interface locating algorithm for free surface two-phase flows,” Journal of xi’an Jiaotong University, vol. 51, no. 1, pp. 79–87, 2017.
View at: Google Scholar
R. Issa, “Solution of the implicitly discretised fluid flow equations by operator splitting,” Journal of Computational Physics, vol. 62, no. 1, pp. 41–65, 1985.
View at: Google Scholar
H. Jasak, H. G. Weller, and A. D. Gosman, “High resolution NVD differencing scheme for arbitrarily unstructured meshes,” International Journal for Numerical Methods in Fluids, vol. 31, no. 2, pp. 431–449, 1999.
View at: Publisher Site | Google Scholar
J. Crank and N. Phyllis, “A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type,” Mathematical Proceedings of the Cambridge Philosophical Society, vol. 43, no. 1, pp. 50–67, 1947.
View at: Publisher Site | Google Scholar
R. L. Burden and J. D. Faires, Numerical Analysis, PWS Publishers, 3rd ed edition, 1985.
P. N. Sun, Y. B. Li, and F. R. Ming, “Numerical simulation on the motion characteristics of freely rising bubbles using smoothed particle hydrodynamics method,” Acta Physica Sinica, vol. 64, no. 17, pp. 1747011–1747015, 2015.
View at: Google Scholar
F. Davidson, “The initial motion of a gas bubble formed in an inviscid liquid Part 1. The two-dimensional bubble,” Journal of Fluid Mechanics, vol. 12, no. 3, pp. 408–461, 1962.
View at: Google Scholar

Copyright

Copyright © 2023 Qiliang Mei et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

202

Downloads

251

Citations