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

Preprint

Article

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

Altmetrics

Downloads

133

Views

Comments

A peer-reviewed article of this preprint also exists.

Simon EL OUAFA^*,

Stephane VINCENT,Vincent LE CHENADEC,Syphax FEREKA,

Benoît TROUETTE,

Amine CHADIL

Simon EL OUAFA^*,

Stephane VINCENT,Vincent LE CHENADEC,Syphax FEREKA,

Benoît TROUETTE,

Amine CHADIL

This version is not peer-reviewed

Submitted:

01 February 2024

Posted:

08 February 2024

You are already at the latest version

Alerts

Abstract

In this paper, we introduce a compressible formulation for dealing with 2D/3D compressible interfacial flows. The formulation is combined with an original monolithic solver, according to which a strong velocity-pressure coupling is achieved. The results show that the scheme can ensure accuracy and stability for various fluid flow scenarios: from the standard incompressible conditions to compressible flows, spanning both single-phase and two-phase flows.

Keywords:

Subject: Physical Sciences - Other

1. Introduction

The modeling and simulation of two-phase compressible flows is a highly dynamic field of study due to their crucial involvement in energy systems, such as aerospace engineering, oil and gas industry, nuclear and biomedical engineering, geophysical studies, and chemical processes.

The use of numerical simulation proves necessary for analysing and understanding compressible two-phase flows involving separated phases. Firstly, these flows are intrinsically complex and involve challenging interactions that occur between different phases of fluids or gases in particular in the vicinity of interfaces. Additionally, safety considerations also drive the use of the numerical simulation. In applications like nuclear reactors or chemical processes, conducting experiments may not be safe or feasible. Further, the numerical simulation facilitates parametric studies, enabling engineers to efficiently explore a large range of operating conditions and design parameters. The reasons mentioned above have prompted the CFD community to develop a two-phase compressible flow model throughout the last decades. Our efforts have been directed towards encompassing a wide array of Mach formulations documented in prior research. Noteworthy references include the work of [1,2,3,4,5,6], and other relevant sources.

In this paper, a compressible formulation is developed to simulate such flow situations. In this formulation, we solve the conservation of mass (in two forms), momentum, and total energy in each of the two phases, as well as an equation for the volume fraction. To close the system of equations, an equation of state is used to take into account the variations in density as a function of pressure and possibly temperature. The original formulation was introduced by [7] and recently extended to all-Mach flows by [8] to account for heat diffusion between two different compressible phases. In the present work, we extend this formulation by maintaining a strong coupling between velocity and pressure variables through the use of a monolithic solver [9,10,11] and by preserving the consistency between mass transport and momentum via a momentum conserving scheme [12].

The remainder of the manuscript is structured as follows: Section 2 introduces the governing equations, including conservation of mass, momentum, and total energy in each of the two phases. Section 3 describes the employed numerical schemes as well as the monolithic solver used for the solution of the saddle point system on the velocity-pressure couple. In Section 4, we introduce different cases to check how accurate is the incorporation of our compressible formulation. These test cases will be divided into two parts: one to deal with an adiabatic case for an inviscid flow, such as Sod’s shock tube problem whereas a second part will be devoted to isothermal cases for viscous flows, like the liquid injection in a closed cavity, compression of an air bubble by water, and the drop impact on viscous liquid film. Ultimately, we will sum up our work and offer a perspective on future research.

2. Governing Equations

The derivation of the compressible formulation follows the work of [7]. The governing equations, suitable for modelling compressible two-phase flows using a 1-fluid model, are presented. In their conservative forms, the mass and momentum equations, read:

\frac{\partial ρ}{\partial t} + \nabla \cdot (ρ v) = 0

(1)

and

\frac{\partial ρ v}{\partial t} + \nabla \cdot (ρ v \otimes v) = - \nabla p + \nabla \cdot \bar{\bar{τ}} + ρ g + F_{s}

(2)

where

v

is the fluid velocity, p the pressure field, t the time, and

ρ

and

μ

the properties of the fluid. In addition,

\bar{\bar{τ}} = μ (\nabla v + {(\nabla v)}^{T}) - \frac{2}{3} μ \nabla \cdot v

is the viscous stress tensor,

g

is the gravity acceleration, whereas

F_{s}

is the capillary term acting on the interface, modeled in this study by the continuum surface tension force (CSF) [13]. The tracking of the spatio-temporal evolution of the interface requires solving an advection equation for the phase indicator color function C:

\frac{\partial C}{\partial t} + v \cdot \nabla C = 0

(3)

By definition,

C = 1

in one phase and 0 in the other phase. In the framework of a finite volume approximation of the solution, the color function in the control volumes cut by the interface is such as

0 < C < 1

. As the velocity and the pressure fields are coupled, a relation between

v

and p is needed to obtain the evolution equation for the pressure. Thereby, the conservation of mass equation is not used in its original form (Equation (1)) but is transformed into a pressure equation that is combined with the velocity variable. Within the framework of a compressible flow between two successive instants t and

t + d t

, the differential of the density with respect to the pressure p and the temperature T leads to:

\frac{d ρ}{d t} = {(\frac{\partial ρ}{\partial p})}_{T} \frac{d p}{d t} + {(\frac{\partial ρ}{\partial T})}_{p} \frac{d T}{d t}

(4)

Introducing the coefficients of isothermal compressibility and isobaric thermal expansivity (or sometimes called the expansion or dilatation coefficient), respectively

χ_{T} = ρ^{- 1} {(\partial ρ / \partial p)}_{T}

and

β = ρ^{- 1} {(\partial ρ / \partial T)}_{p}

, the mass conservation (Equation (1)), combined with Equation (4), can be rewritten as

χ_{T} \frac{d p}{d t} + β \frac{d T}{d t} + \nabla \cdot v = 0

(5)

For incompressible (

χ_{T} = 0

) and isothermal (

β = 0

) flows, Equation (5) gives

\nabla \cdot v = 0

When considering compressible flow, the energy equation is also addressed. Within the scope of this paper, we focus on scenarios where thermal diffusion and mass transfer at the interface do not play a significant role and are then neglected. Consequently, under these specific circumstances, the total energy equation is formulated as follow:

\frac{\partial ρ e}{\partial t} + \nabla \cdot (ρ e v) = - \nabla \cdot (v p) + \nabla \cdot (\bar{\bar{τ}} \cdot v) + ρ g \cdot v + F_{s} \cdot v

(6)

where

e = u + e_{k}

denotes the total energy, which is the sum of the internal energy u and kinetic energy per unit mass

e_{k} = {∥ v ∥}^{2} / 2

To complete the set of equations, an equation of state (EoS) that establishes a relationship among the thermodynamic variables pressure p, density

ρ

and temperature T is required. For any phase, gas or liquid, the Noble-Abel Stiffened-Gas (NASG) [14] can give general formulation of an EoS. In this work, liquid phase is always assumed incompressible while the ideal gas model is adopted for the gas phase. The variation of the density as a function of pressure and temperature is classically expressed by

ρ = p / (r T)

wherein r is the specific gas constant. In the other hand, from the NASG EoS of an ideal gas, it comes

ρ = p / (u (γ - 1))

, with

γ = c_{p} / c_{v}

the isotropic gas coefficient. The combination of the two previous expressions provides a relation for the temperature as a function of the internal energy:

T = r u (γ - 1)

needed to close the system of Equations (5), (2) and (6), where the interface dynamics is provided by Equation (3).

3. Numerical Scheme

In this section, the global algorithm used to solve the coupled mass, momentum, and energy equations detailed in the previous section is presented.

The algorithm is designed with the underlying idea of a full implicit formulation of system equations. For example, in momentum conservation Equation (2), all variables (

ρ

v

, ...) and terms would like to be resolved simultaneously. As it is difficult to achieve because of non-linearity, some quantities, such as the physical properties, the inertial term or the geometrical properties of the interface are estimated by values that are expected to be good approximations of the ones that would have been obtained by an implicit solving. To that purpose, a consistent reformulation and discretization of the inertial term, based on a momentum preserving approach [12], is used. The spacial discretization relies of a classical conservative finite volume approach and is not detailed in the paper. In practice, solving the full system introduce a sequential resolution of equations that themselves combine explicit variables (inertial contribution) with implicit variables resolved by inversion of a linear systems.

After time discretization, all variables at time

t^{n} = t_{0} + n Δ t

, where n is the iteration number and

Δ t

the constant (non restrictive) time step, are supposed known from a previous solution, directly (

v^{n}

p^{n}

C^{n}

and

e^{n}

) or from a reconstruction (

ρ^{n}

μ^{n}

, ...). At this step, density, and more generally all the physical properties are synchronised with the phase indicator function

C^{n}

. The density is for example deduced from a mixing rule,

ρ^{n} = (1 - C^{n}) ρ_{1} + C^{n} ρ_{2}

, where

ρ_{i}

is density of phase i. According to the phase state, an EoS is used to specify the behaviour of

ρ_{i}

The initial step involves inertial term computation of Equations (2) and (6). As the l.h.s. of mentioned equations has the same mathematical structure $\partial ϕ / \partial t + \nabla \cdot (ϕ v)$ , with $ϕ = ρ v$ , $ρ e$ , a general approach is used to compute temporary variables (denoted $ϕ^{★}$ ) in the operator splitting framework [12]. As the density is also a required variable, Equation (1) is also used in the numerical scheme in order to provide an approximation of the density.

In practice, $ϕ^{n} = ρ^{n}$ , $ρ^{n} v^{n}$ , $ρ^{n} e^{n}$ , is first initialised before the time integration using a third-order accurate strong stability preserving Runge-Kutta (SSP-RK3) time integrator [15]:

$ϕ^{(1)} = ϕ^{n} - Δ t \nabla \cdot (ϕ^{n} v^{n})$

(7a)

$ϕ^{(2)} = \frac{3}{4} ϕ^{n} + \frac{1}{4} ϕ^{(1)} - \frac{1}{4} Δ t \nabla \cdot (ϕ^{(1)} v_{a d v}^{(1)})$

(7b)

$ϕ^{★} = \frac{1}{3} ϕ^{n} + \frac{2}{3} ϕ^{(2)} - \frac{2}{3} Δ t \nabla \cdot (ϕ^{(2)} v_{a d v}^{(2)})$

(7c)

where $v_{a d v}^{(1)} = 2 v^{n} - v^{n - 1}$ and $v_{a d v}^{(2)} = (3 v^{n} - v^{n - 1}) / 2$ are the velocities correctly extrapolated to maintain scheme accuracy [16]. To ensure that $ϕ$ remains bounded during advection, high-order schemes such as CUI, WENO, or even CUBISTA are used (see [16] and references herein for more details).

At the end of this step, discrete consistency between temporary density $ρ^{★}$ , momentum ${(ρ v)}^{★}$ , and energy ${(ρ e)}^{★}$ variables is ensured. Furthermore, as Equation (1) is a pure advection equation, $ρ^{★}$ provides directly a predicted density: $ρ^{n + 1} = ρ^{★}$ .
The next step consists in the color function C advection. A conservative VOF approach, proposed by [17], is used where Equation (3), with time discretization, is formulated as

$\frac{C^{n + 1} - C^{n}}{Δ t} + \nabla \cdot (C^{n} v_{a d v}^{(2)}) = C^{n} \nabla \cdot v_{a d v}^{(2)}$

(8)
With the knowledge of $C^{n + 1}$ , the termophysical properties of the 1-fluid model are updated. For $ξ = μ$ , $χ_{T}$ , $β$ , $γ$ and r by using an arithmetic mixing law:

$ξ^{n + 1} = (1 - C^{n + 1}) ξ_{1} + C^{n + 1} ξ_{2}$

(9)

where $ξ_{i}$ denotes the property corresponding to phase i. Note that the density is not updated as its value is already known form step 1.
Evaluate the geometrical properties, normal $n$ and curvature $κ$ , of the interface using. Here a smoothed color function ${\tilde{C}}^{n + 1}$ is computed from a diffusion step applied on $C^{n + 1}$ and the definition

$n^{n + 1} = \frac{\nabla {\tilde{C}}^{n + 1}}{∥ \nabla {\tilde{C}}^{n + 1} ∥} and κ^{n + 1} = \nabla \cdot n^{n + 1}$

(10)

can be used in the CSF expression of $F_{s} = σ κ n δ_{I}$ where $σ$ is the surface tension and $δ_{I}$ the interface localisation.
From the EoS, here of an ideal gas, the temperature T can be expressed as a function solved variables $v$ and e:

$T = \frac{γ - 1}{r} (e - \frac{{∥ v ∥}^{2}}{2})$

(11)

and can be injected in the mass conservation (Equation (5)) that couples solved quantities $v$ , p and e. However, the velocity norm prevents an implicit coupling due to the non linearity. The total derivative $d T / d t$ then will be explicited in the following for non isothermal flows.
Using the new density $ρ^{n + 1}$ and temporary momentum ${(ρ v)}^{★}$ from step 1, physical and interface properties from steps 3 and 4 and temperature definition from step 5, the mass conservation, augmented by the compressibility and dilatation effects, and the momentum equations reads, with a first order time discretization,

$χ_{T}^{n + 1} (\frac{p^{n + 1} - p^{n}}{Δ t} + v^{n + 1} \cdot \nabla p^{n}) + β^{n + 1} (\frac{T^{n} - T^{n - 1}}{Δ t} + v^{n + 1} \cdot \nabla T^{n}) + \nabla \cdot v^{n + 1} = 0$

(12a)

$\frac{\begin{matrix} ρ^{n + 1} \end{matrix} v^{n + 1} - {(ρ v)}^{★}}{Δ t} = - \nabla p^{n + 1} + \nabla \cdot {\bar{\bar{τ}}}^{n + 1} + ρ^{n + 1} g + σ κ^{n + 1} n^{n + 1} δ_{I}$

(12b)

where only velocity $v^{n + 1}$ and pressure $p^{n + 1}$ fields are implicity coupled and constitute the unknown vector of the underlying linear system. This latter is solved with a BiCGStab(2) solver [18] where an efficient block triangular preconditioner, improving the convergence of the iterative solver as explained in [10,11], is used.

A last discussion concerns the linearization of term $v^{n + 1} \nabla p^{n + 1}$ into $v^{n + 1} \nabla p^{n}$ . This choice relies essentially on the available stencils from the incompressible version of the used solver [11]. Another approach, where $v^{n + 1} \nabla ϕ^{n + 1}$ , $ϕ$ being the pressure or the temperature field, is one again approximated with the application of step 1 to a conservative evolution equation for $ϕ$ , is also considered.
Finally, we update the total energy ${(ρ e)}^{n + 1}$ using the intermediate total energy ${(ρ e)}^{★}$ computed at step 1 and all other variables now known at time $t^{n + 1}$ :

$\frac{{(ρ e)}^{n + 1} - {(ρ e)}^{★}}{Δ t} = - \nabla \cdot (v^{n + 1} p^{n + 1}) + \nabla \cdot ({\bar{\bar{τ}}}^{n + 1} \cdot v^{n + 1}) + ρ^{n + 1} g \cdot v^{n + 1} + σ κ^{n + 1} n^{n + 1} \cdot v^{n + 1} δ_{I}$

(13)

The total energy at the end of the time iteration is deduced from $e^{n + 1} = {(ρ e)}^{n + 1} / ρ^{n + 1}$ .

4. Results

4.1. Sod’s Shock Tube Problem

This configuration was first introduced by [19]. This is a famous test case used to check the ability of the compressible schemes to avoid numerical instability at the shock level and improve the capture of a complex system of shock waves inside the tube. The governing equations are the Euler equations (Equations (12a) and (12b) with

μ = 0

) supplemented by the equation for the conservation of total energy (Equation (6)), with air characterized by an ideal gas model.

Here, a tube of length

L = 1

m is considered. A membrane splits the tube into two regions: a high-pressure gas (

p_{l e f t}, ρ_{l e f t}

) on the left side and a low-pressure gas (

p_{r i g h t}, ρ_{r i g h t}

) on the other. The membrane is removed at time

t = 0

. A complex system of shock waves develops inside the tube, i.e., an expansion wave, a contact discontinuity, and a normal shock wave.

The initial condition is given on the 1D computational domain,

0 \leq x \leq 1

m, as follows:

(ρ, v, p) = \{\begin{matrix} (ρ_{l e f t}, v_{l e f t}, p_{l e f t}) = (1, 0, 1) & if x \leq 0.5 \\ (ρ_{r i g h t}, v_{r i g h t}, p_{r i g h t}) = (0.125, 0, 0.1) & if x > 0.5 \end{matrix}

(14)

Concerning the numerical parameters, the simulations are carried out on a 250, 500 and 1000 cartesian grid with time step conditions such that the CFL, defined by

max | u | Δ t / Δ x = 0.3

. Whereby the matrices resulting from discretization will be solved with the BiCGSTAB(2) solver with a residual of

ε = 10^{- 9}

The Figure 1 demonstrates well-resolved shocks with correct shock locations and fewer smeared contact discontinuities. The corresponding solutions show the evolution of pressure and density at different times (

0 s

0.025 s

0.05 s

0.075 s

, and

0.1 s

) as a function of the curvilinear abscissa of the shock tube. From this figure, it is clear that an expansion wave propagates to the left, a shock wave propagates to the right, and there is a discontinuity between these two waves.

Since the analytical solution was available, the theoretical curves are overlayed on the numerical ones in Figure 2. This illustrates a clear convergence of the numerical simulations to the exact solution when using a total energy equation, whether for velocity, pressure, internal energy, or density.

Finally, the choice of the solved energy equation is discussed here. As mentioned in Section 3, as the total energy conservation is used, the temperature is deduced from the fluid EoS properties. Other choice are possible, such as the enthalpy conservation, or in a more traditional way, the energy equation formulated with the temperature variable (heat or thermal energy equation). The two first choices involve dealing with conserved quantities, i.e. the total energy e or the entalpy h, while the temperature variable in the heat equation is not a conserved quantity. In the latter case, at the discrete level, the energy conservation is then excepted and no longer ensured by the formulation. The consequences are significant, in particular for shocks. The Figure 3 presents again pressure, density, velocity and internal energy fields, by resolving total energy (Equation (6)) or the heat equation (see [7]). On one hand, and as shown before with the convergence study (Figure 3), conservation of total energy ensures that all quantities are predicted accurately, even for coarse grids. On the other hand, on the finest grid level, using the thermal energy equation implies that all quantities across the shocks are badly predicted.

4.2. Isothermal Case for A Viscous Flow without Capilary Effects

In order to validate the compressible two-phase model presented in the previous section, two academic test cases are studied [3] in the context of isothermal and viscous flows. The purpose is to check the accuracy of the compressible model in a simple isothermal two-phase case by comparing the numerical density and pressure fields with the analytical solutions. The relative error is then estimated for various meshes and convergence order can be extracted. In both cases, the thermophysical characteristics of the two fluids are given in Table 1.

The two configurations are quite similar and are presented in Figure 4.

In a square cavity of side length

L = 0.1

m, air is compressed by a water injection at a constant mass flow rate. Air initial density is

ρ_{0}

. As analytical solution of theses problems are based on a quasi static evolution hypothesis, the velocity

V_{0}

is chosen to be very low. In configuration (a), air initially spans a length

L_{0}

and the water injection velocity is

V_{0} = 0.1

m/s. The mass of air is constant and its volume will vary over time with the volume of water injected according to the following law:

ρ (t) = ρ_{0} {(1 - \frac{V_{0} t}{L_{0}})}^{- 1}

(15)

In the second case (b), an air bubble of initial radius

R_{0} = 30

mm, is compressed by injecting water from all the sides of the square cavity, with

V_{0} = 2.5 \times 10^{- 3}

m/s. The theoretical equation for density evolution over time is:

ρ (t) = ρ_{0} {(1 - \frac{4 V_{0} L t}{π R_{0}^{2}})}^{- 1}

(16)

The analytical pressure field p is obtained from the analytical density using the ideal gas model:

p (t) = ρ (t) r T_{0}

, with

T_{0} = 300 K

the reference temperature.

In the following, the cases are simulated on four cartesian grids meshes composed of

32^{2}

64^{2}, 128^{2}

and

256^{2}

control volumes, with a constant time step

Δ t = 10^{- 4}

s and a residual of

ε = 10^{- 9}

for the iterative solver. Averaged values of solved density and pressure are performed for the gas phase and comprared to the analytical solutions (15) and (16).

The Figure 5 presents the density evolution over time. The final time is chosen such as the final density value (or the gas volume) is about twice (half) the initial one. Before discussions about the convergence study of the simulations, it should be noted that all quantities of interest (density and pressure in the air) well converge, for all times, through the references solution The convergence with the number of control volume is by upper values in case 1 and by lower values in case 2.

Errors on density and pressure fields, estimated with a L1 relative norm, is presented in Table 2 for both cases.

Regarding the relative errors, a Richardson extrapolation is also used on the three finest meshes to find the asymptotic solutions at final times. A first order 1 convergence towards the asymptotic values is again obtained.

Note the convergence orders of pressure and density differ slightly since the errors use quantities at the end of the time step. Density is obtained through step 1 in the global algorithm presented in Section 3 while the pressure field comes from the resolution of the momentum and mass conservation equations (Equations (12a) and (12b)). In the case the fields are synchronised before the error computation, the order is strictly the same.

4.3. Isothermal Case for A Viscous Flow with Capilary Effects: Drop Impact on Viscous Liquid Film

As a continuation of the test cases already presented, the drop impact on the viscous liquid film case is set up here. The main aim is to check the ability of our compressible formulation to simulate two-phase flows with large density and viscosity ratios and strong interface deformation. Here, the Mach number is around

0.007

and the flow configuration is almost incompressible. Following the simulations in [20], we consider a drop of water of diameter

D = 2.8 \times 10^{- 3} m

, perfectly circular, with coordinates

(0, 0, 5.15 D)

, placed into a rectangular domain of dimensions

L x = 10 D

L y = 10 D

, and

L z = 20 D

. The water surface is initialized at y = 4.5D. The boundary conditions of the domain are non-slip on all faces and homogeneous Neumann on the top. The time step is adaptive, and the CFL is

0.3

. The simulations are carried out with the residual

ϵ = 10^{- 5}

of BiCGSTAB(2)solver, both with the compressible formulation and with the incompressible scheme.

Figure 6. Experimental images [20] of the crater formation inside a deep pool and the rising crown above the liquid surface. The views are respectively at 4.3 ms, 22.6 ms and 28.7 ms.

Figure 7. Simulation of drop impact into a deep pool using different schemes with a resolution of 32 cells per drop diameter. Presentation of the VOF interface (

C = 0.5

isosurface) at three different times: 4.3 ms, 22.6 ms and 28.7 ms. The simulation with a resolution of 32 cells per drop diameter was carried out on 512 processors of the TGCC IRENE Rome HPC cluster.

Figure 7. Simulation of drop impact into a deep pool using different schemes with a resolution of 32 cells per drop diameter. Presentation of the VOF interface (

C = 0.5

isosurface) at three different times: 4.3 ms, 22.6 ms and 28.7 ms. The simulation with a resolution of 32 cells per drop diameter was carried out on 512 processors of the TGCC IRENE Rome HPC cluster.

The results obtained from the compressible solver are compared with those from the incompressible solver and the experimental results shown in Figure 8. The figure shows time snapshots captured at 5, 25, 45, and 65 milliseconds after droplet impact. Initially, as inertia dominates capillary forces, the drop penetrates the water film and induces the formation of a crater (see Figure 8a) and a rising corona, which later induces a secondary ejection of droplets. After reaching maximum depth, the crater begins to retract over time due to capillary forces. Finally, compressible and incompressible solvers give identical results that are in agreement with the experiment. In addition, in Figure 8b, we can see an equivalence between the compressible model and the incompressible model on the evolution of the crater, in which a good agreement is observed between numerical results and experimental measurements. We are dealing with the case where the compressible model degenerates into an incompressible model. We then conclude that our compressible scheme is capable of handling two-phase flows with large interface deformations in incompressible limit of the flow motion.

5. Conclusions

In this work, we introduce a compressible formulation integrated with an original monolithic solver to address two and three dimensional compressible interfacial flow, showcasing a concrete progression in computational fluid dynamics. Our developed scheme stands as a notable advancement, providing a dependable and flexible toolkit for accurately and stably simulating diverse fluid flow scenarios, spanning from incompressible to compressible and encompassing both single-phase and two-phase situations with or without capilary effects. Future works will be devoted to the extension of this solver to the validation of the presented model in the case of droplet impact on viscous liquid film at high velocity or also a spray of drops and the interaction between shock/rarefaction waves and fluid/fluid interfaces.

Funding

This research received no external funding.

Acknowledgments

We are grateful for access to the computational facilities of theFrench CINES (National computing center for higher education) andTGCC granted by GENCI, France under project numbers A0112b06115.We thank the technical and administrative teams of these supercom-puter centers and agencies for their kind and efficient help. All authorsread and approved the final version of manuscript to be published.

Conflicts of Interest

The authors declare no conflict of interest.

References

Yoon, S.Y.; Yabe, T. The unified simulation for incompressible and compressible flow by the predictor-corrector scheme based on the CIP method. Computer Physics Communications 1999, 119, 149–158. [Google Scholar] [CrossRef]
Kwatra, N.; Su, J.; Grétarsson, J.T.; Fedkiw, R. A method for avoiding the acoustic time step restriction in compressible flow. Journal of Computational Physics 2009, 228, 4146–4161. [Google Scholar] [CrossRef]
Caltagirone, J.P.; Vincent, S.; Caruyer, C. A multiphase compressible model for the simulation of multiphase flows. Computers & fluids 2011, 50, 24–34. [Google Scholar] [CrossRef]
Shyue, K.M.; Xiao, F. An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC approach. Journal of Computational Physics 2014, 268, 326–354. [Google Scholar] [CrossRef]
Nourgaliev, R.; Luo, H.; Weston, B.; Anderson, A.; Schofield, S.; Dunn, T.; Delplanque, J.P. Fully-implicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change. Journal of Computational Physics 2016, 305, 964–996. [Google Scholar] [CrossRef]
Urbano, A.; Bibal, M.; Tanguy, S. A semi implicit compressible solver for two-phase flows of real fluids. Journal of Computational Physics 2022, 456, 111034. [Google Scholar] [CrossRef]
Fuster, D.; Popinet, S. An all-Mach method for the simulation of bubble dynamics problems in the presence of surface tension. Journal of Computational Physics 2018, 374, 752–768. [Google Scholar] [CrossRef]
Saade, Y.; Lohse, D.; Fuster, D. A multigrid solver for the coupled pressure-temperature equations in an all-Mach solver with VoF. Journal of Computational Physics 2023, 476, 111865. [Google Scholar] [CrossRef]
El Ouafa, M.; Vincent, S.; Le Chenadec, V. Navier-stokes solvers for incompressible single- and two-phase flows. Communications in Computational Physics 2021, 29, 1213–1245. [Google Scholar] [CrossRef]
El Ouafa, M.; Vincent, S.; Le Chenadec, V. Monolithic solvers for incompressible two-phase flows at large density and viscosity ratios. Fluids 2021, 6, 23. [Google Scholar] [CrossRef]
El Ouafa, S.; Vincent, S.; Le Chenadec, V.; Trouette, B. Fully-coupled parallel solver for the simulation of two-phase incompressible flows. Computers & Fluids 2023, 265, 105995. [Google Scholar] [CrossRef]
El Ouafa, M. Développement d’un solveur tout-couplé parallèle 3D pour la simulation des écoulements diphasiques incompressibles à forts rapports de viscosités et de masses volumiques. Theses, Université Gustave Eiffel, 2022.
Brackbill, J.U.; Kothe, D.B.; Zemach, C. A continuum method for modeling surface tension. Journal of computational physics 1992, 100, 335–354. [Google Scholar] [CrossRef]
Le Métayer, O.; Saurel, R. The Noble-Abel stiffened-gas equation of state. Physics of Fluids 2016, 28. [Google Scholar] [CrossRef]
Strong Stability-Preserving High-Order Time Discretization Methods. SIAM Review 2001, 43, 89–112. [CrossRef]
Nangia, N.; Griffith, B.E.; Patankar, N.A.; Bhalla, A.P.S. A robust incompressible Navier-Stokes solver for high density ratio multiphase flows. Journal of Computational Physics 2019, 390, 548–594. [Google Scholar] [CrossRef]
Weymouth, G.D.; Yue, D.K.P. Conservative volume-of-fluid method for free-surface simulations on cartesian-grids. Journal of Computational Physics 2010, 229, 2853–2865. [Google Scholar] [CrossRef]
Dongarra, J.-J.; Duff, L.-S.; Sorensen, D.C.; Vorst, H.A.V. Numerical Linear Algebra for High Performance Computers; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 1998. [Google Scholar]
Sod, G.A. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. Journal of computational physics 1978, 27, 1–31. [Google Scholar] [CrossRef]
Bisighini, A.; Cossali, G.E.; Tropea, C.; Roisman, I.V. Crater evolution after the impact of a drop onto a semi-infinite liquid target. Physical Review E 2010, 82, 036319. [Google Scholar] [CrossRef]

Figure 1. Sod’s shock tube problem: The different color lines represent the solutions on a 1000 grid at different times, as indicated by the legend. From the left to the right plot, we display density and pressure along the horizontal axis.

Figure 2. Sod’s shock tube problem: Comparison of the numerical solution on different grids (using total energy) against the exact solution. The solutions are compared at a physical time of 0.1 seconds and in a line along the horizontal axis.

Figure 3. Sod’s shock tube problem: Comparison of the numerical solution on a 1000 grid (using total and thermal energy) against the exact solution. The solutions are compared at a physical time of 0.1 seconds and in a line along the horizontal axis.

Figure 4. (a) Water injection in a closed cavity initially full of air. (b) Compression of an air bubble by water injection.

Figure 5. Density over time for both cases and different grids. Final time is chosen such as such that the final density is twice the initial one.

Figure 8. Evolution of the dimensionless crater depth

Z_{c r}^{*}

as function of dimensionless time

τ

for the studied schemes. The numerical results are compared to the experimental measurement of [20].

Figure 8. Evolution of the dimensionless crater depth

Z_{c r}^{*}

as function of dimensionless time

τ

for the studied schemes. The numerical results are compared to the experimental measurement of [20].

Table 1. Physical properties of air (initial) and water used in the computations. Here, water is considered to have a constant density whereas air is treated as a compressible fluid with an ideal gas law that couples pressure and density. The surface tension between to two phases is neglected.

	Air	Water
Density $ρ$ (kg $\cdot$ $m^{- 3}$ )	$\begin{matrix} 1.1768292 \end{matrix}$	1000
Viscosity $μ$ (Pa $\cdot$ s)	$\begin{matrix} 1.85 e - 5 \end{matrix}$	$\begin{matrix} 1 e - 3 \end{matrix}$
Compressibility $χ_{T}$ ( ${Pa}^{- 1}$ )	$\begin{matrix} 9.8692322 e - 6 \end{matrix}$	$\begin{matrix} 0.44 e - 9 \end{matrix}$
Specific gas constant r (J $\cdot$ K/kg)	287	−

Table 2. Convergence study for both case. Relative error in

L_{1}

norm is computed from analytical relations Equation (15) and (16).

Table 2. Convergence study for both case. Relative error in

L_{1}

norm is computed from analytical relations Equation (15) and (16).


Case 1

Mesh	$ρ$	Error $L_{1}$	order	p	Error $L_{1}$	order

$32 \times 32$	2.571	2.172 × 10⁻²		2.217 × 10⁵	2.114 × 10⁻²
$64 \times 64$	2.544	1.093 × 10⁻²	0.98	2.194 × 10⁵	1.034 × 10⁻²	1.03
$128 \times 128$	2.531	5.580 × 10⁻³	0.97	2.182 × 10⁵	4.980 × 10⁻³	1.05
$256 \times 256$	2.523812	2.840 × 10⁻³	0.97	2.176 × 10⁵	2.240 × 10⁻³	1.15
Extrapolation	2.517	3.730 × 10⁻⁴	0.97	2.170 × 10⁻⁵	1.200 × 10⁻³	0.97

Case 2

Mesh	$ρ$	Error $L_{1}$	order	p	Error $L_{1}$	order

$32 \times 32$	2.383	4.784 × 10⁻²		2.123 × 10⁵	1.652 × 10⁻²
$64 \times 64$	2.441	2.482 × 10⁻²	0.94	2.149 × 10⁵	4.310 × 10⁻³	0.93
$128 \times 128$	2.474	1.137 × 10⁻²	1.12	2.182 × 10⁵	4.980 × 10⁻³	1.00
$256 \times 256$	2.490	5.210 × 10⁻³	1.12	2.188 × 10⁵	2.030 × 10⁻³	1.08
Extrapolation	2.503	1.630 × 10⁻⁵	1.12	2.189 × 10⁵	1.640 × 10⁻³	1.04

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

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

Abstract

1. Introduction

2. Governing Equations

3. Numerical Scheme

4. Results

4.1. Sod’s Shock Tube Problem

4.2. Isothermal Case for A Viscous Flow without Capilary Effects

4.3. Isothermal Case for A Viscous Flow with Capilary Effects: Drop Impact on Viscous Liquid Film

5. Conclusions

Funding

Acknowledgments

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe