A consistent rescaled momentum transport method for simulating large density ratio incompressible multiphase flows using level set methods

Many multiphase flows occur at ambient conditions resulting in large density ratios between the fluids and involve significant shear at the phase interface. A prime example for such a scenario involves the atomization of liquids. Although many atomization devices for technical applications under real operating conditions operate in pressurized gaseous environments, such devices are often studied experimentally under ambient pressure conditions to lower experimental cost. As such, significant more experimental data exists for ambient, i.e. large density ratio conditions than for high pressure, i.e. low density ratio conditions. Experimental datasets necessary for validation of simulation tools are thus more readily available and often of higher fidelity for high density ratio than for low density ratio applications.

The case is reversed for detailed numerical simulations of atomizing flows. There, a large density ratio can result in a more stiff system of equations that are more costly to solve. Unfortunately, some classes of numerical techniques describing the dynamics of the phase interface in incompressible flows, mostly level set based methods, are prone to numerical instabilities if the liquid-to-gas density ratio is large, i.e. of the order of 100 or more, and the flow exhibits a significant shear at the phase interface, common to many atomization devices. This numerical instability manifests itself in a sudden spike in local velocity that can grow unbounded. A good review of past work on different numerical methods to overcome this issue can be found in [2, 14].

Interestingly enough, some numerical methods to describe the phase interface motion appear not to be susceptible to the numerical instability, among them the geometric transport volume of fluid (VOF) methods, see [16] and references therein. The important difference of these methods to most level set based approaches lies not only in the fact that the VOF approaches solve the momentum equation in conservative form, but more importantly that they employ discrete operators for the convection terms in the momentum and VOF-scalar equation that are discretely identical. They thus ensure that mass, in the form of the VOF scalar, and momentum are transported in exactly the same discrete manner. For example Rudman [15] introduced a method to solve multiphase flows with high density ratio using a projection method on staggered grids with mass conservation based on a volume tracking method using a grid twice as fine as the velocity-pressure grid. Consistent mass and momentum advection was achieved by solving the conservative form of the momentum advection term with mass flux densities obtained from the volume-of-fluid method. Bussmann et al [1] extended Rudman's algorithm to three-dimensional unstructured flow solvers, using a collocated arrangement of variables on a single mesh.

Standard level set based methods, on the other hand, transport mass and momentum in entirely different ways. Mass is transported by solving the level set equation, whereas momentum is transported using typically a non-conservative formulation of the Navier–Stokes equations. Thus, even a small error in the position of the phase interface can lead to strong generation of artificial momentum in the presence of large density ratios and large shear, see figure 1. A literature survey shows that different numerical methods have been suggested to overcome this issue.

Zoom In Zoom Out Reset image size

Smiljanovski et al [17] proposed an in-cell reconstruction hybrid tracking/capturing method for deflagration discontinuities. With it, the front geometry is explicitly computed using a level set scalar. This is then used to reconstruct consistent values in both fluids using the known jump conditions across the interface and consistent cell face fluxes using only values in the respective fluids. The method was successfully applied to density discontinuities in the context of Richmyer–Meshkov instabilities [4].

Hu et al [6] constructed a method based on a standard Cartesian finite volume method and the level set technique for multi-fluid problems with complex boundaries. Two sets of equations are solved for each fluid, and an interface interaction method [5] is used to exchange momentum between the two fluids. Cell-face apertures are calculated according to the level set distribution along the Cartesian grid cell faces, in order to determine the amount of momentum transferred along the interface.

To avoid the artificial generation of momentum depicted in figure 1 one could of course switch from the non-conservative form of the Navier–Stokes equations to the conservative form and use a discretely conservative numerical scheme. However, this approach is equally bound to fail, since the density necessary to reconstruct velocity from momentum is again prone to phase interface position errors resulting in large over-shoots in velocity. Minimizing position errors of the phase interface can alleviate the problem, however, even if a level set method were mass conserving, there still exists the mechanism of artificial momentum/velocity creation since momentum and mass are not guaranteed to be transported in a discretely consistent manner.

The key in avoiding the numerical instability is thus to ensure a discretely consistent transport of mass and momentum. To this end, Sussman et al [19] suggested a method that uses the coupled level set and volume-of-fluid method [18], where they extend the velocity of the fluid with the higher density to the fluid with the lower density in a band around the fluid interface and store two different velocity fields. Transport of level set and volume-of-fluid functions are done using the extrapolated (higher density fluid) velocity.

Raessi [12] and Raessi and Pitsch [14] introduced a method to construct flux densities [15] from level set scalar information and use these flux densities to transport momentum in a consistent manner. However, their method is presently limited to one- and two-dimensional cases and not straightforward to extend to three dimensions.

In this paper we introduce an alternative approach to ensure discretely consistent transport of mass and momentum for level set based methods that is simple to implement, viable in three dimensions, applicable to unstructured, collocated finite volume formulations of the governing equations and consistent with the balanced force formulation for level set methods [3].

The equations governing the motion of an unsteady, incompressible, immiscible, two-fluid system are the Navier–Stokes equations in conservative form

$$\begin{equation} \frac{\partial \rho \boldsymbol{u}}{\partial t} + \nabla \cdot ( \rho \boldsymbol{u} \boldsymbol{u}) = - \nabla p + \nabla \cdot \left( \mu \left( \nabla \boldsymbol{u} + \nabla^{\mathrm{T}} \boldsymbol{u} \right) \right)+ \rho \boldsymbol{g} + \boldsymbol{T}_{\sigma} \end{equation} \tag{ 1 }$$

or in non-conservative form

$$\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \cdot \nabla \boldsymbol{u} = - \frac{1}{\rho} \nabla p + \frac{1}{\rho} \nabla \cdot \left( \mu \left( \nabla \boldsymbol{u} + \nabla^{\mathrm{T}} \boldsymbol{u} \right) \right) + \boldsymbol{g} + \frac{1}{\rho} \boldsymbol{T}_{\sigma}, \end{equation} \tag{ 2 }$$

where u is the velocity, ρ the density, p the pressure, μ the dynamic viscosity, g the gravitational acceleration, and T_σ the surface tension force which is non-zero only at the location of the phase interface x_f,

$$\begin{equation} \boldsymbol{T}_{\sigma} (\boldsymbol{x}) = \sigma \kappa \delta(\boldsymbol{x}-\boldsymbol{x}_{f}) \boldsymbol{n} , \end{equation} \tag{ 3 }$$

with σ the assumed constant surface tension coefficient, κ the local mean surface curvature, n the local surface normal, and δ the delta-function. In addition, conservation of mass results in the continuity equation

$$\begin{equation} \frac{\partial \rho}{\partial t} +\nabla \cdot \left(\rho \boldsymbol{u}\right) = 0 . \end{equation} \tag{ 4 }$$

The phase interface location x_f between the two fluids is described by a level set scalar G, with

$$\begin{equation} G(\boldsymbol{x}_{f},t) = 0 \end{equation} \tag{ 5 }$$

at the interface, G(x,t) > 0 in fluid 1, and G(x,t) < 0 in fluid 2. Differentiating (5) with respect to time yields the level set equation

$$\begin{equation} \frac{\partial G}{\partial t} + \boldsymbol{u} \cdot \nabla G = 0. \end{equation} \tag{ 6 }$$

For numerical accuracy of geometric properties of the phase interface it is advantageous, although not necessary, to define the level set scalar away from the interface to be a signed distance function

$$\begin{equation} |\nabla G| = 1. \end{equation} \tag{ 7 }$$

Assuming ρ and μ constant within each fluid, density and viscosity at any point x can be calculated from

$$\begin{equation} \rho(\boldsymbol{x}) = H(G) \rho_1 + (1-H(G)) \rho_2, \end{equation} \tag{ 8 }$$

$$\begin{equation} \mu(\boldsymbol{x}) = H(G) \mu_1 + (1-H(G)) \mu_2, \end{equation} \tag{ 9 }$$

where indices 1 and 2 denote values in fluid 1, respectively 2, and H is the Heaviside function. Finally, the interface normal vector n and the interface curvature κ can be expressed in terms of the level set scalar as

$$\begin{equation} \boldsymbol{n} = \frac{\nabla G}{| \nabla G|} \, , \quad \kappa = \nabla \cdot \boldsymbol{n}. \end{equation} \tag{ 10 }$$

The Navier–Stokes equations are solved using a fractional step method [7] on unstructured collocated meshes using a finite volume approach, where control volume material properties like density and viscosity are defined using equations (8) and (9) as

$$\begin{equation} \rho_{\mathrm{cv}} = \psi_{\mathrm{cv}} \rho_1 + (1-\psi_{\mathrm{cv}}) \rho_2, \end{equation} \tag{ 11 }$$

$$\begin{equation} \mu_{\mathrm{cv}} = \psi_{\mathrm{cv}} \mu_1 + (1-\psi_{\mathrm{cv}}) \mu_2 , \end{equation} \tag{ 12 }$$

with the control volume fraction ψ_cv given by

$$\begin{equation} \psi_{\mathrm{cv}} = 1/V_{\mathrm{cv}} \int_{V_{\mathrm{cv}}}\! H(G)\,\mathrm{d}V . \end{equation} \tag{ 13 }$$

Here V_cv is the volume of the control volume.

In the consistent rescaled momentum transport (CRMT) method, we first solve the continuity equation, equation (4). In discrete form this results in

$$\begin{equation} V_{\mathrm{cv}} \frac{{{\rho }}^*_{\mathrm{cv}} - {{\rho }}^n_{\mathrm{cv}}}{\Delta t} + \sum_f u_f^{n} {{\rho }}_f^\prime A_{f} = 0 , \end{equation} \tag{ 14 }$$

where A_f is the cell face area, uⁿ_f is the face normal velocity, and ρⁿ_cv is calculated from equations (11) and (13) using the level set solution at time tⁿ. In equation (14), ρ'_f is defined as

$$\begin{equation} {\rho}_{f}^\prime = \left\{ \begin{array}{@{}ll@{}} \rho ^{n}_{\rm Upwind}, & \epsilon < \psi_{\mathrm{cv}}\! < 1- \epsilon \\ & \mathrm{or} \quad \epsilon < \psi_{\mathrm{nbr}} < 1- \epsilon, \\ \displaystyle\frac{\rho^n_{\mathrm{cv}}+\rho^n_{\mathrm{nbr}}}{2} & {\textrm{elsewhere}}. \end{array} \right. \, \end{equation} \tag{ 15 }$$

Here is a small number, the index nbr denotes the neighbor control volume to cv sharing the same face, and ρⁿ_Upwind is calculated using a simple first-order upwind approach

$$\begin{equation} \rho^{n}_{\rm Upwind} = \left\{ \begin{array}{@{}ll@{}} \rho^{n}_{\mathrm{cv}}, & u_{f} \geqslant 0, \\ \rho^{n}_{\mathrm{nbr}}, & u_{f} < 0 . \end{array} \right. \, \end{equation} \tag{ 16 }$$

The reason for using a simple first-order approach here lies in the fact that this will guarantee boundedness of ρ*_cv, provided the face normal velocities uⁿ_f are discretely divergence free. The choice of ρ'_f away from the phase interface in equation (15) is dictated by the discrete operator for the momentum equation away from the interface described below and does not result in an unconditionally unstable method for equation (14), since away from the phase interface ρⁿ_cv = ρⁿ_nbr.

To conserve momentum discretely, we choose the conservative form of the Navier–Stokes equations. The discrete form of the conservative Navier–Stokes equations on collocated unstructured meshes reads

$$\begin{eqnarray} &&V_{\mathrm{cv}} \frac{({{\rho u})}^{*}_{i,\mathrm{cv}} - ({{\rho u})}^n_{i,\mathrm{cv}}}{\Delta t} + \sum_{f} u_{f}^{n} ({{\rho u})}^{\prime}_{i,f} A_{f}\nonumber\\ &&\quad = \sum_{f} A_{f} \mu_{f} (({\nabla u})_{i,f} + ({{\nabla u}})^{\mathrm{T}}_{i,f}) +V_{\mathrm{cv}} \rho_{\mathrm{cv}}^{n} g_{i}+V_{\mathrm{cv}} \rho_{\mathrm{cv}}^{n} F^{n}_{i,\mathrm{cv}} ,\nonumber\\ \end{eqnarray} \tag{ 17 }$$

where Fⁿ_i,cv is a surface tension force induced acceleration, the subindex i indicates a spatial direction, μ_f = (μ_cv + μ_nbr)/2, and (ρu)ⁿ_i,cv = ρⁿ_cvuⁿ_i,cv.

To ensure discrete consistency with equation (14), (ρu)'_i,f in equation (17) is defined as

$$\begin{eqnarray} ({{\rho u})}^\prime_{i,f} = \left\{ \begin{array}{@{}ll@{}} ({{\rho u})}^{n}_{i,\mathrm{Upwind}}, & \epsilon < \psi_{\mathrm{cv}} < 1- \epsilon \\ & \mathrm{or} \quad \epsilon < \psi_{\mathrm{nbr}} < 1- \epsilon, \\ \\ \displaystyle\frac{({{\rho u})}^n_{i,\mathrm{cv}}+({{\rho u})}^n_{i,\mathrm{nbr}}}{2} & {\textrm{elsewhere,}} \end{array} \right. \, \end{eqnarray} \tag{ 18 }$$

with (ρu)ⁿ_i,Upwind calculated using the first-order upwind scheme

$$\begin{equation} ({{\rho u})}^{n}_{i,\mathrm{Upwind}} = \left\{ \begin{array}{@{}cl} ({{\rho u})}^{n}_{i,\mathrm{cv}}, & u_{f} \geqslant 0, \\ ({{\rho u})}^{n}_{i,\mathrm{nbr}}, & u_{f} < 0. \end{array} \right. \, \end{equation} \tag{ 19 }$$

This ensures that the resulting method is discretely energy conserving away from the phase interface and discretely identical to the continuity equation method.

To calculate Fⁿ_i,cv in equation (17), we follow the balanced force approach for collocated, finite volume methods [3]. At the cell face, the surface tension force is

$$\begin{equation} T_{\sigma_{f}}^{n} = \sigma \kappa_{f}^{n} (\nabla \psi)_{f}^{n} , \end{equation} \tag{ 20 }$$

resulting in

$$\begin{equation} F_{f}^{n} = T_{\sigma_{f}}^{n}/\rho_{f}^{n} , \end{equation} \tag{ 21 }$$

with ρⁿ_f = (ρⁿ_cv + ρⁿ_nbr)/2. To ensure discrete consistency between the surface tension force at the control volume centroid and the pressure gradient evaluated there, Fⁿ_i,cv is calculated from Fⁿ_f using a face-area weighted least-squares method [9] by minimizing

$$\begin{equation} \epsilon_{\mathrm{cv}} = \sum_f \left( F_{i,\mathrm{cv}}^{n} n_{i,f} - F_{f}^{n}\right)^2 A_{f} , \end{equation} \tag{ 22 }$$

where n_i,f is the face normal vector.

Using the solutions to equations (17) and (14), (ρu)^⋆_i,cv and ρ^⋆_cv, we can then calculate the predicted velocity

$$\begin{equation} u_{i,\mathrm{cv}}^{\star}=\frac{{(\rho u)}_{i,\mathrm{cv}}^{\star}}{{\rho}_{\mathrm{cv}}^{\star}}. \end{equation} \tag{ 23 }$$

Next, we project the predicted velocity field u^⋆_i,cv into the subspace of divergence free velocity fields, by first solving the pressure Poisson equation

$$\begin{equation} \nabla \cdot \left(\frac{1}{{\rho}^{n}}\nabla {p}\right)=\frac{1}{\Delta t} \nabla \cdot {\boldsymbol{u}}^{\star} \end{equation} \tag{ 24 }$$

which in discrete form is

$$\begin{equation} \sum_f \frac{1}{\rho_{f}^{n}} \frac{\partial p_{\mathrm{cv}}}{\partial n} A_{f} = \frac{1}{\Delta t} \sum_f u_{f}^* A_{f} , \end{equation} \tag{ 25 }$$

with

$$\begin{eqnarray} &&u_{f}^* = \textstyle\frac{1}{2}\left( u_{i,\mathrm{cv}}^* + u_{i,\mathrm{nbr}}^* \right) n_{i,f} - \textstyle\frac{1}{2} \Delta t \left( F_{i,\mathrm{cv}}^{n} + F_{i,\mathrm{nbr}}^{n} \right) n_{i,f} + \Delta t F_{f}^{n}.\nonumber\\ \end{eqnarray} \tag{ 26 }$$

The projection, i.e. correction step, is then

$$\begin{equation} {\boldsymbol{u}}^{n+1}={\boldsymbol{u}}^{\star}- \frac{\Delta t}{{\rho}^{n}} \nabla{p} \end{equation} \tag{ 27 }$$

which for the face velocities in discrete form is

$$\begin{equation} u_{f}^{n+1} = u_{f}^* - \Delta t P_{f} , \end{equation} \tag{ 28 }$$

with

$$\begin{equation} P_{f} = \frac{1}{\rho_{f}^{n}}(\nabla p_{\mathrm{cv}}^{n})_{f} = \frac{1}{\rho_{f}^{n}}\frac{p_{\mathrm{nbr}}^{n} - p_{\mathrm{cv}}^{n}}{| \boldsymbol{s}_{\mathrm{cv},\mathrm{nbr}}| }. \end{equation} \tag{ 29 }$$

Here, s_cv,nbr is the vector connecting the cv and nbr control volume centroids.

To correct the control volume velocities, first the control volume centroid-based density weighted pressure gradient P_cv is calculated from the face-based density weighted gradient P_f using the same face-area weighted least-squares method employed in calculating F_f, see equation (22),

$$\begin{equation} \epsilon_{\mathrm{cv}} = \sum_f \left( P_{i,\mathrm{cv}} n_{i,f} - P_f\right)^2 A_f . \end{equation} \tag{ 30 }$$

Then, the control volume centroid velocity is corrected as

$$\begin{equation} u^{n+1}_{i,\mathrm{cv}} = u^*_{i,\mathrm{cv}} - \Delta t P_{i,\mathrm{cv}} . \end{equation} \tag{ 31 }$$

Finally, we discard the solution to the continuity equation, ρ*, and reset the density at tⁿ⁺¹ using the level set solution Gⁿ⁺¹ obtained from solving the level set equation (6), as

$$\begin{eqnarray} &&\rho_{\mathrm{cv}}^{n+1} = \psi_{\mathrm{cv}}^{n+1} \rho_1 + (1-\psi_{\mathrm{cv}}^{n+1}) \rho_2. \end{eqnarray} \tag{ 32 }$$

This density is then used to update the momentum as

$$\begin{equation} ({{\rho u})}^{n+1}_{i,\mathrm{cv}} = {{\rho}}^{n+1}_{\mathrm{cv}}{{ u}}^{n+1}_{i,\mathrm{cv}}. \end{equation} \tag{ 33 }$$

In this section, we present results obtained with the new method for a range of test cases involving large fluid to gas density ratios. In all cases, we track the position of the phase interface using the refined level set grid method [3].

4.1. Collapse of a water column

A 2D water column with initial height and width of h = 5.715 cm is placed inside a container of size 40 × 10 cm as shown in figure 2. The density of the water and air are ρ_l = 1000 kg m⁻³, respective ρ_g = 1.226 kg m⁻³, the viscosities are μ_l = 1.137 × 10⁻³ kg ms⁻¹, respective μ_g = 1.78 × 10⁻⁵ kg ms⁻¹, the surface tension coefficient is σ = 0.0728 N m⁻¹, and the gravitational acceleration is g = −9.81 m s⁻².

Zoom In Zoom Out Reset image size

Figure 2 shows the phase interface shape at Δt = 0.1 s, time intervals obtained using the CRMT method compared to the shape obtained using a standard non-conservative method [3]. The non-conservative method shows significantly slower lateral spread of the water column and some unphysical deformations of the phase interface as compared to the results of the CRMT method. The CRMT results are comparable to the improvements reported using the method of Raessi [12] and Raessi and Pitsch [13].

The non-dimensional front position and non-dimensional height of the water column as a function of non-dimensional time are shown in figure 3, where reference length is h and the reference time is $\sqrt {\frac {h}{g}}$. Figure 3 also shows the results of a grid refinement study resolving the container by 512 × 128, 1024 × 256, and 2048 × 512 equisized hexahedral control volumes. The lateral front position converges under grid refinement, however the results appear to converge to a slightly faster spread rate as that observed experimentally [10], whereas the height of the water column is well captured even on the coarsest mesh.

Zoom In Zoom Out Reset image size

Table 1 shows both the order of convergence for the relative error of the front position (Z) with respect to the experimental data, $E_{\mathrm {exp}}=\frac {Z_{\mathrm {exp}}-Z}{Z_{\mathrm {exp}}}$ and with respect to the finest grid solution, $E_{\mathrm {finest}}=\frac {Z_{\mathrm {finest}}-Z}{Z_{\mathrm {finest}}}$ at time t = 3.9.

Table 1. Error of the front position in the collapse of the water column.

Grid	Eexp	Order	Efinest	Order
128×32	0.356	–	0.360	–
256×64	0.190	0.90	0.195	0.88
512×128	0.058	1.71	0.064	1.60
1024×256	0.015	1.95	0.021	1.60
2048×512	0.006	1.32	–	–

4.2. Damped surface wave

The dynamics of a small amplitude damped surface wave between two superposed immiscible fluids are described by the initial value theory of Prosperetti [11]. The initial surface position inside a [0,2π] × [0,2π] box is given by a sinusoidal disturbance of wavelength λ = 2π and amplitude A₀ = 0.01λ,

$$\begin{equation} G(\boldsymbol{x}, t=0) = y-y_0 + A_0 \cos (x-h_G/2) , \end{equation} \tag{ 34 }$$

with y₀ = π. Periodic boundary conditions are used in the x-direction and slip walls are imposed in the y-direction. The initial value solution for two fluids with equal kinematic viscosity ν and λ = 2π can be written as [11]

$$\begin{eqnarray} A_{\mathrm{ex}}(t) &=& \frac{4(1-4\beta)\nu^2}{8(1-4\beta)\nu^2 +\omega_0^2} A_0 \,\mathrm{erfc} \sqrt{\nu t} \nonumber\\ & &+ \sum_{i=1}^4 \frac{z_i}{Z_i}\left(\frac{\omega_0^2 A_0}{z_i^2-\nu }\right) \exp[(z_i^2-\nu ) t] \mathrm{erfc} (z_i \sqrt{t}) ,\nonumber\\ \end{eqnarray} \tag{ 35 }$$

where z_i are the roots of

$$\begin{eqnarray} &&z^4-4 \beta \sqrt{\nu} z^3 + 2 (1-6 \beta) \nu z^2+ 4 (1-3 \beta) \nu^{3/2} z\nonumber\\ &&\quad + (1-4 \beta)\nu^2 +\omega_0^2 = 0 , \end{eqnarray} \tag{ 36 }$$

the dimensionless parameter β is given by β = ρ₁ρ₂/(ρ₁ + ρ₂)², the inviscid oscillation frequency is $\omega _0 = \sqrt {\frac {\sigma }{\rho _1+\rho _2}}$, and $Z_i = \prod _{j=1 \atop j \neq i}^4 (z_j-z_i)$. In [3] results for ρ_l = 1000, ρ_g = 1, σ = 2 and ν_l = ν_g = 0.064 720 863 were reported using the non-conservative formulation.

Figure 4 shows the temporal evolution of the non-dimensional disturbance amplitude A for a mesh consisting of 128 × 128 equisized hexahedra using both the CRMT and non-conservative method [3], including a zoom of the temporal evolution starting at t = 130 to more clearly see the difference in results. The CRMT method shows noticeably improved results compared to the non-conservative methods.

Zoom In Zoom Out Reset image size

Figure 5 shows the results of a grid refinement study in a zoom of the temporal evolution starting at t = 130. Excellent agreement of the CRMT method with the analytical results can be seen. Figure 6 shows the evolution of the corresponding non-dimensional error E(t) = (A(t) − A_ex(t))/A₀ for hexahedral and prism meshes, while table 2 summarizes their root mean squares. At the same grid resolution, the CRMT method shows a significantly lower error as compared to the non-conservative formulation [3].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Rms of amplitude error for damped surface wave.

	CRMT	Ref. [3]
Cartesian mesh
32	1.440×10−2	4.82×10−2
64	4.534×10−3	2.08×10−2
128	2.465×10−3	1.27×10−2
256	2.246×10−3	1.18×10−2
Prism mesh
32	1.263×10−2	5.64×10−2
64	6.743×10−3	1.41×10−2
128	6.188×10−3	1.13×10−2
256	4.703×10−3	1.57×10−2

4.3. Zero gravity column oscillation

To further verify the implementation of the new method, this section presents results for zero gravity oscillating columns. The theoretical oscillation period for columns in the linear regime is given by [8]

$$\begin{equation} \omega^2 = \frac{n ( n^2-1) \sigma}{(\rho_1 + \rho_2) R_0^3} . \end{equation} \tag{ 37 }$$

In all simulations, a column of radius R₀ = 2 is placed in the center of a [ − 10,10] square box with periodic boundary conditions on all sides and σ = 1, ρ₁ = 1, ρ₂ = 0.01, μ₁ = 0.01, and μ₂ = 1 × 10⁻⁴, resulting in a Laplace number of La = 20 000. The column is initially perturbed by a mode n = 2 perturbation with an initial amplitude of A₀ = 0.01R₀. The time step size in all simulations is chosen as Δt = 0.5Δt_cap. Where

$$\begin{equation} t_{\mathrm{cap}} =\sqrt{\frac{(\rho_1 + \rho_2)h^3}{4 \pi \sigma}}. \end{equation} \tag{ 38 }$$

Table 3 shows the period of oscillation error E_T  = |T_calcω/2π − 1| for the oscillating column together with the results reported in [3]. On fine hexahedral meshes, the CRMT method gives noticeably improved results as compared to those of the non-conservative method, whereas the results on prism meshes are comparable.

Table 3. Zero gravity 2D column oscillation. Error in oscillation period as compared to linear theory [8].

Grid size	ET (hexahedra)	Ref. [3]	ET (prism)	Ref. [3]
20/64	7.47×10−2	4.04×10−2	5.66×10−2	5.91×10−2
20/128	7.32×10−3	1.05×10−2	1.61×10−2	1.65×10−2
20/256	3.44×10−4	3.7×10−3	1.35×10−2	1.36×10−2

4.4. Convection of high density droplet

In this test case initially proposed by Bussman et al [1], a 2D liquid droplet of diameter D = 0.4 is placed in the center of a 1 × 1 periodic domain filled with gas. The density ratio is chosen as ρ_l/ρ_g = 10⁶ and the fluids are assumed inviscid and without surface tension. The drop is given an initial homogeneous velocity of u = (1,0) while the gas is initially at rest. We have employed different structured equisized hexahedral and unstructured prism grids. Because of the large density ratio, the drop should stay essentially undeformed while passing through the computational domain multiple times.

Figure 7 shows the drop shape after passing the domain once (t ≈ 1) using the non-conservative method [3] on a 128 × 128 hexahedral mesh. Erroneous transfer of momentum from the liquid to the gas has caused significant interface deformation, resulting in an unphysical shattering of the drop. Figure 8 shows the drop shape after passing the entire domain once (t = 1) obtained using the CRMT method together with the expected solution for varying mesh resolutions and table 4 summarizes the corresponding shape errors. While there are some minor deformations of the drop visible, the drop stays nearly spherical and no erroneous large scale interface deformations are visible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Shape errors for convection of high density droplet.

Grid size	ET (hexahedra)	Order	ET (prism)	Order
1/32	5.33×10−3	–	5.42×10−3	–
1/64	3.73×10−3	0.51	2.93×10−3	0.88
1/128	2.64×10−3	0.49	1.65×10−3	0.82
1/256	1.85×10−3	0.51	–	–

A new CRMT scheme for modeling incompressible, multiphase flows with high density ratios in the context of level set interface capturing methods has been presented. In this method, the conservative form of the Navier–Stokes equations is solved using an unstructured, collocated, incompressible, fractional step flow solver. Instead of replacing the continuity equation in its entirety by a level set equation to keep track of the phase interface separating fluids of different, but constant, density, we solve the continuity equation to obtain a predicted density using operators that are discretely consistent with those used in the solution of the conservative momentum equation. Using the predicted momentum and density, we recover a predicted velocity that is then projected into the subspace of divergence free velocity fields. To avoid undue dissipative errors from solving the continuity equation directly, we then reset the density and momentum according to the density obtained from a level set solution.

It should be pointed out that this new approach is neither mass nor momentum conserving. The key idea is in fact to allow mass and momentum errors, where the former are unavoidable in a pure level set method, but to ensure that they are discretely consistent.

The new method has shown excellent results for a range of test cases involving large density ratios up to 10⁶ and offers a path to simulate atomizing flows with large density ratio fluids.

This work was supported by NSF grant number CBET-0853627. The authors would also like to thank M Raessi, O Desjardins and F Ham for many helpful discussions.