An accelerated sharp-interface method for multiphase flows simulations

TL;DR

This paper introduces an accelerated sharp-interface method for multiphase flow simulations that allows different fluids to be updated with distinct time steps, significantly improving computational efficiency while maintaining accuracy.

Contribution

The novel method enables separate time stepping for each fluid in multiphase simulations, enhancing efficiency without sacrificing stability or accuracy.

Findings

Achieves comparable accuracy to traditional methods

Provides significant computational speedup

Maintains conservative properties with flux correction

Abstract

In this work, we develop an accelerated sharp-interface method based on (Hu et al., JCP, 2006) and (Luo et al., JCP, 2015) for multiphase flows simulations. Traditional multiphase simulation methods use the minimum time step of all fluids obtained according to CFL conditions to evolve the fluid states, which limits the computational efficiency, as the sound speed c of one fluid may be much larger than the others. To address this issue, based on the original GFM-like sharp interface methods, the present method is developed by solving the governing equations of each individual fluid with the corresponding time step. Without violating the numerical stability requirement, the states of fluid with larger time-scale features will be updated with a larger time step. The interaction step between two fluids is solved for synchronization, which is handled by interpolating the intermediate states of fluid with larger time-scale features. In addition, an interfacial flux correction is implemented to maintain the conservative property. The present method can be combined with a wavelet-based adaptive multi-resolution algorithm (Han et al., JCP, 2014) to achieve additional computational efficiency. A number of numerical tests indicate that the accuracy of the results obtained by the present method is comparable to the original costly method, with a significant speedup.

Figures

Tables

Table 1. Table 1: Periodic motion of a water droplet: Speedup ratios for different resolutions. In this paper, for all cases, the speedup performance of the accelerated method is obtained on a desktop with the Intel Core i7-4790 Processor (8 M Cache, 3.60 GHz, 8 cpus, 8GB of RAM). The intel Threading Building Blocks(TBB) library [ 5 ] is used to map logical tasks to physical threads.

Resolution	Calculation time ($s$)		Speedup ratio $r_s$
	Original	Modified
Uniform grid - 2D
$32\times 32$	140.585	91.3541	1.539
$64\times 64$	618.656	321.757	1.923
$128\times 128$	3736.34	1609.05	2.322
$256\times 256$	26601.7	9259.04	2.873
Uniform grid - 3D
$128\times 128\times 128$	1001450.6	349911.1	2.862
Multi-resolution grid - 2D
$256\times 256$	12309.6	5970.85	2.062

Table 2. Table 2: Static drop with surface tension: Speedup ratios for different resolutions.

Resolution	Calculation time($s$)		Speedup ratio $r_s$
	Original	Modified
Uniform grid
$32\times 32$	230.491	144.468	1.595
$64\times 64$	993.682	461.515	2.153
$128\times 128$	6200.05	2449.27	2.531
$256\times 256$	43680.1	14068.5	3.1048
Multi-resolution grid
$256\times 256$	9441.02	19285.4	2.043

Table 3. Table 3: Recovering a circle shape: Speedup ratios for different resolutions.

Resolution	Calculation time($s$)		Speedup ratio $r_s$
	Original	Modified
Uniform grid - 2D
$32\times 32$	112.889	71.778	1.573
$64\times 64$	604.497	348.205	1.736
$128\times 128$	3141.19	1270.66	2.472
$256\times 256$	21889.7	7168.9	3.053
Uniform grid - 3D
$128\times 128\times 128$	946140.3	328314.1	2.882
Multi-resolution grid - 2D
$256\times 256$	9206.3	5174.4	1.779

Table 4. Table 4: Liquid filament contraction: Speedup ratios for different resolutions

Resolution	Calculation time($s$)		Speedup ratio $r_s$
	Original	Modified
Uniform grid - 2D
$16\times 32$	236.415	175.885	1.344
$32\times 64$	930.915	230.498	1.661
$64\times 128$	4325.97	2116.39	2.044
$128\times 256$	28849.9	11673.9	2.471
Uniform grid - 3D
$64\times 64\times 128$	668940.1	275521.6	2.428
Multi-resolution grid - 2D
$128\times 256$	17085.6	9735.21	1.755

Table 5. Table 5: Equilibrium shape of a water drop resting on a wall: Speedup ratios for different static contact angles

Static contact angle(${\theta }_s$)	${30}^{\circ }$	${45}^{\circ }$	${60}^{\circ }$	${90}^{\circ }$	${120}^{\circ }$	${150}^{\circ }$
Uniform grid
Speedup ratio	2.171	2.078	2.039	2.0	2.162	2.139
Multi-resolution grid
Speedup ratio	1.714	1.724	1.726	1.717	1.718	1.723

Equations

$$\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot \mathbf{F}=0,$$

$${\alpha }_{i,j}^{n+1}{\mathbf{U}}_{i,j}^{n+1}=$$

$$+\frac{\Delta t}{\Delta y}[A_{i,j-1/2}{\hat{\mathbf{F}}}_{i,j-1/2}-A_{i,j+1/2}{\hat{\mathbf{F}}}_{i,j+1/2}].$$

$$\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot \mathbf{F}=\nabla \cdot {\mathbf{F}}_v+{\mathbf{a}}_s,$$

$$p$$

$$B$$

$${\alpha }_{i,j}^{n+1}{\mathbf{U}}_{i,j}^{n+1}=$$

$$+\frac{\Delta t}{\Delta x}[A_{i-1/2,j}{\hat{\mathbf{F}}}_{i-1/2,j}-A_{i+1/2,j}{\hat{\mathbf{F}}}_{i+1/2,j}]$$

$$+\frac{\Delta t}{\Delta y}[A_{i,j-1/2}{\hat{\mathbf{F}}}_{i,j-1/2}-A_{i,j+1/2}{\hat{\mathbf{F}}}_{i,j+1/2}]$$

$$+\frac{\Delta t}{\Delta x}[A_{i-1/2,j}{\widehat{{\mathbf{F}}_v}}_{i-1/2,j}-A_{i+1/2,j}{\widehat{{\mathbf{F}}_v}}_{i+1/2,j}]$$

$$+\frac{\Delta t}{\Delta y}[A_{i,j-1/2}{\widehat{{\mathbf{F}}_v}}_{i,j-1/2}-A_{i,j+1/2}{\widehat{{\mathbf{F}}_v}}_{i,j+1/2}].$$

$$\frac{\partial \Phi }{\partial t}+\mathbf{u}\cdot \nabla \Phi =0,$$

$$(\alpha \mathbf{U}{)}_1^{interpolation}=(1-\frac{t_2}{t_1})(\alpha \mathbf{U}{)}_1^{previous}+\frac{t_2}{t_1}(\alpha \mathbf{U}{)}_1^{future},$$

$${\hat{\mathbf{X}}}^{correction}$$

$${\widehat{{\mathbf{X}}_v}}^{correction}$$

$${\hat{\mathbf{X}}}^{correction}$$

$${\widehat{{\mathbf{X}}_v}}^{correction}$$

$$[p{]}_{\Gamma }=p_{I,1}-p_{I,2}=\sigma \kappa ,$$

$${\widehat{{\mathbf{X}}_s}}^{correction}=({\widehat{{\mathbf{X}}_s}}^{correction}{)}^{\ast }-0.5\times \widehat{{\mathbf{X}}_1}\times {\Delta t}_1+0.5\times \frac{p_{I,1}}{p_{I,2}}\times \widehat{{\mathbf{X}}_2}\times {\Delta t}_2$$

$${\widehat{{\mathbf{X}}_s}}^{correction}=({\widehat{{\mathbf{X}}_s}}^{correction}{)}^{\ast }+0.5\times \frac{p_{I,1}}{p_{I,2}}\times \widehat{{\mathbf{X}}_2}\times {\Delta t}_2$$

$$r_s=\frac{n}{[n]-([n]-1)m},$$

$$\{\begin{array}{rl}if& I\%2^{L_{max}-l}=0\text{leaves of fluid1 at the level l},\\ if& I\%(2^{L_{max}-l}\times n)=0\text{leaves of fluid2 at the level l},\end{array}$$

$$\{\begin{array}{rl}if& (I+1)\%2^{L_{max}-l}=0\text{leaves of fluid1 at the level l},\\ if& (I+1)\%(2^{L_{max}-l}\times n)=0\text{leaves of fluid2 at the level l}.\end{array}$$

$$(\rho ,u,p,\gamma )=\{\begin{array}{rlr}& (0.01,0,1000,2& )ifx<0.5,\\ & (1,0,1,7.15& )ifx>0.5,\end{array}$$

$$(\rho ,u,p,\gamma )=\{\begin{array}{rlr}& (1.0376,6.0151,1000,7.15& )if0<x<0.7,\\ & (0.001,0,1,1.4& )if1>x>0.7,\end{array}$$

$$\{\begin{array}{rlr}\rho & =1.2,p=1,u=\frac{\sqrt{10}}{100},v=0,\gamma =1.4& air,\\ \rho & =1000,p=1,u=\frac{\sqrt{10}}{100},v=0,\gamma =7.15& water,\\ \varphi & =-0.1+\sqrt{(x-0.5{)}^2+(y-0.5{)}^2}& levelset,\end{array}$$

$$\frac{M^n}{M^0}=\frac{{\sum }_{i=1}^{i=N}{\alpha }^n{\rho }_i^n}{{\sum }_{i=1}^{i=N}{\alpha }^0{\rho }_i^0},$$

$$\{\begin{array}{rlr}\rho & =1,p_0=10^3,B=100,\gamma =1.4,\mu =1.6\times 10^{-5}& air,\\ \rho & =1000,p_0=10^3,B=5.0\times 10^5,\gamma =7.15,\mu =10^{-3}& liquiddrop,\\ \varphi & =-0.1+\sqrt{(x-0.5{)}^2+(y-0.5{)}^2}& levelset.\end{array}$$

$$\Delta p_{exact}=p_{liquid}-p_{air}=\frac{\sigma }{R}=360,$$

$$E(\Delta p)=\frac{|\Delta p_{exact}-\Delta p_{numerical}|}{|\Delta p_{exact}|}.$$

$$\{\begin{array}{rlr}\rho & =1,p_0=10^3,B=100,\gamma =1.4,\mu =1.6\times 10^{-5}& gas,\\ \rho & =1000,p_0=10^3,B=5.0\times 10^5,\gamma =7.15,\mu =10^{-3}& liquiddrop,\\ \varphi & =max(|x-0.5|,|y-0.5|)-0.1& levelset.\end{array}$$

$$\{\begin{array}{rlr}\rho & =1,p_0=10^4,B=500,\gamma =1.4& gas,\\ \rho & =100,p_0=10^4,B=2.0\times 10^5,\gamma =7.15& liquiddrop,\\ \varphi & =-1.0+\sqrt{(x-0.5{)}^2/0.2^2+(y-0.5{)}^2/0.12^2}& levelset.\end{array}$$

$${\omega }^2=(o^3-o)\frac{\sigma }{({\rho }_l+{\rho }_g)R^3},T=\frac{2\pi }{\omega }.$$

$$\left\{\begin{aligned} \rho&=1\ kg/m^{3},p_{0}=10^{4}Pa,B=5.0\times 10^{2}Pa,\quad\text{air}\\ \gamma&=1.4,\mu=1.6\times 10^{-5}Pa\ s\\ \rho&=1000\ kg/m^{3},p_{0}=10^{4}Pa,B=3.0\times 10^{6}Pa,\quad\text{liquid\ filament}\\ \gamma&=7.15,\mu=1\times 10^{-3}Pa\ s\\ \phi&=\left\{\begin{array}[]{ll}{-0.19+|x|\ mm}&{z<1.52}\\ {-0.19+\sqrt{x^{2}+(z-1.52)^{2}}\ mm}&{z\geq 1.52}\end{array}\right.\quad\text{ level-set. }\end{aligned}\right.$$

Full text

1 Introduction

Nowadays, the study of multiphase flows has attracted much attention from academia due to its widely practical usage in many branches of industry such as the bubble column reactors [25], the solidification process of metals [41], the civitation in the naval industry [3], etc. The main challenge in the numerical simulation of multiphase flows is the treatment of the material interface which separates interacting different fluids.

In general, the approaches for modelling interface interactions can be divided into two classes, the diffuse-interface methods and the sharp-interface methods. In the diffuse-interface methods, with limited thickness, the interface is often defined as an artificial region where the transition from one phase to another is smooth. Among those diffuse-interface methods, volume of fluid (VOF) methods [32, 30, 10, 33, 29] may be the most popular ones. The VOF methods use the volume of fluid function, which represents the fraction of the liquid volume in a given computational cell, to capture the interface position. The VOF function is updated by solving the advection equation with the finite volume scheme so that the mass conservation is ensured. Nevertheless, there are some problems in computing the geometrical variables such as the normal vector and curvature because the VOF function is not continuous.

The other approach is to track the interface with a non-smearing representation, i.e., localize the interface as well as the interface interactions at an infinitely thin region. In the arbitrary Lagrangian Eulerian(ALE) methods [12, 19], the interface is represented by the computational mesh which evolves and deforms with the flow while it is represented by interface markers in the front-tracking method [39]. Introduced by Osher et al. [27], the level-set methods share algorithmic similarities with the diffuse-interface methods and are simple to be implemented. Different from those sharp-interface mehtods mentiond above, the level-set methods employ a signed distance function known as the level-set function to represent the interface implicitly and are able to handle interface deformations and topological changes in a straightforward way [38, 26]. The main drawback of the level-set methods is the lack of discrete conservation properties. To address this issue, Hu et al. [14] introduce a fully conservative level-set based sharp-interface method for compressible flows. The standard finite volume scheme on Cartesian grids is used for the far interface region while it is slightly modified for the near interface region. Unlike those hybrid methods such as the particle level-set methods [7, 40] and the coupled level-set and volume of fluid methods(CLSVOF) [24, 37], this method maintains the simplicity of the pure level-set method and can be simply implemented in multi-dimension and multi-level time integrations. In the method of Hu et al., the treatments of surface tension and the viscous force at the interface singularity are not considered, which often play important roles in mulitphase flows. It is a challenge for numerical models of multiphase flows to handle the singular surface tension force only acting at the interface and the jump of material properties across the interface such as density and viscosity. By employing a weakly compressible model, Luo et al. [22] propose a conservative sharp-interface method based on the work of Hu et al. for incompressible multiphase flows. In the cell cut by the interface, an effective viscosity is constructed to model the viscous flux across the interface and a constrained Riemann problem [2] is solved to model the interfacial flux caused by surface tension.

When explicit time discretization is adopted, the simulation suffers from the limitation of time step due to numerical stability reasons, i.e., CFL conditions such as $\Delta t\propto\frac{1}{|\bm{u}|+c}$ for advection terms and $\Delta t\propto\frac{\rho\Delta x^{2}}{\mu}$ for diffusive terms. In many multiphase flows problems, the time-scale features of different fluids are very different because of the differences of density, sound speed, viscosity, etc. A bottleneck of traditional multiphase simulation methods is that the fluid with small time-scale features imposes a small time step which is also used to evovle the states of fluid with larger time-scale features. To alleviate the stringent CFL restriction, the semi-implicit methods [18, 16, 42] are proposed which are suitable for flows containing low Mach number regions. The calculation is divided into two parts: advection and non-advection. Since the non-advection terms containing the acoustic components are solved implicitly, the time step is constrained by material velocity only while the stringent CFL time step restriction imposed by the acoustic waves is avoided. However, in addition to leading to extra oscillations, these semi-implicit methods require lots of changes to the original code when applied to those multiphase flows solvers based on explicit time discretization. When the interface between the fluids can be considered as a free-boundary, the computation can be limited to the more viscous and dense fluid. For example, in [4, 17], the computation takes place only in the water phase and the computational time for the air phase is saved. Nevertheless, this method is not suitable for problems in which the air phase gets pressurized or the stresses on the liquid caused by the air phase cannot be ignored.

The objective of this paper is to develop an accelerated sharp interface method based on [14, 22] for compressible and weakly compressible multiphase flows. Without generating unphysical oscillations, this method is relatively convenient to be coded up based on the original program. In contrast to the free-surface conditions, the present method makes no assumptions for the gas phase (or the less viscous and dense phase). The main idea of the present method is to obtain speedup by introducing larger time step for the evolution of fluid with larger time-scale features without violating the numerical stability requirements. The intermediate states of fluid with larger time-scale features are obtained by interpolation for the interaction step between two fluids, which is solved for synchronization. In addition, an interfacial flux correction is implemented to maintain the conservative property. The present method can be combined with a wavelet-based adaptive multi-resolution algorithm [11] to achieve additional computational efficiency. For those multiphase flows problems that fluid with lager time-scale features occupies most of the computational domain, significant speedup can be obtained by using the present method because the fluid states of most computational cells will be updated with a larger time step.

The paper is organized as follows: In section 2, we briefly review the original conservative sharp interface methods for compressible and weakly compressible multiphase flows. Subsequently we introduce that how to implement this accelerated method on a uniform grid as well as a multi-resolution grid in section 3. To validate the accelerated method, a number of numerical tests are carried out in section 4. Finally, the concluding remarks are given in section 5.

2 Conservative sharp interface method

2.1 Compressible sharp interface method

For the inviscid and compressible fluid, the governing equations can be written as

[TABLE]

where $\bm{U}={(\rho,\rho\bm{v},E)}^{T}$ is the vector of mass, momentum and total energy densities with the relation $E=\rho e+\frac{1}{2}\rho{|\bm{v}|}^{2}$, and $\bm{F}$ denotes the corresponding physical flux functions of $\bm{U}$. This set of equations is closed by an equation of state (EOS), which gives the thermodynamic properties of the materials. As shown in Fig. 1, the flow domain $\Omega$ is separated by an interface $\Gamma(t)$ into two sub-domains ${\Omega}^{1}$ and ${\Omega}^{2}$, which represent the regions occupied by fluid1 and fluid2, respectively. As in Hu et al. [14], with explicit first-order forward time difference, the discretization of Eq. (1) for each phase on a two-dimensional Cartesian grid with spacing $\Delta x$ and $\Delta y$ can be written as

[TABLE]

Here, $\Delta t$ is the time step size, and ${\alpha}_{i,j}{\bm{U}}_{i,j}$ are the conservative variables in the cut cell, where ${\alpha}_{i,j}$ is volume fraction and ${\bm{U}}_{i,j}$ denotes cell-averaged densities of mass, momentum and energy of the considered fluid. ${A}_{i-1/2,j}$, ${A}_{i+1/2,j}$, ${A}_{i,j-1/2}$ and ${A}_{i,j+1/2}$ are cell-face apertures as shown in Fig. 1. $\hat{\bm{F}}$ denotes the reconstructed cell-face numerical flux obtained by a high-order shock-capturing scheme such as the WENO [35] scheme. $\hat{\bm{X}}(\Delta{\Gamma}_{i,j})$ , where $\Delta{\Gamma}_{i,j}$ is the interface segment within the cut cell, denotes the inviscid momentum and energy flux across the interface determined by the interface interaction.

2.2 Weakly compressible sharp interface method

For the weakly compressible model [22], the governing equations can be written as

[TABLE]

where the right-hand side terms $\bm{F}_{v}$ and $\bm{a}_{s}$, which don’t exist in Eq. (1), are viscous fluxes and fluxes caused by surface tension, respectively. Note that the right-hand side terms in Eq. (3) will vanish when it represents the continuity equation. As shown in Luo et al. [22], Tait’s equation is used as the artificial equation of state to close the governing equations:

[TABLE]

Here, $\gamma$ is the the artificial specific-heat ratio and $\bm{v}_{f}$ represents the characteristic flow velocity. Moreover, $\rho_{0}$ and $p_{0}$ are the reference density and pressure, respectively. The artificial sound speed is $c=\sqrt{\frac{\gamma B}{\rho}}$ and $M=\frac{|\bm{v}_{f}|}{c}$ denotes the Mach number. To impose incompressibility, the Mach number should be small enough, which handled by choosing the parameters $\gamma$, $B$ and $\rho_{0}$ specifically. As in section 2.1, applying explicit first-order forward time difference to Eq. (3) leads to the following form:

[TABLE]

The meanings of the terms which are the same as those in Eq. (2) are not repeated here. In the above equation, $\hat{\bm{F}_{v}}$ denotes the reconstructed cell-face numerical viscous flux obtained by central difference, while $\hat{\bm{X}_{v}}$ and $\hat{\bm{X}_{s}}$ are the fluxes across the interface due to the viscous force and surface tension, respectively. Note that all the interfacial fluxes $(\hat{\bm{X}},\hat{\bm{X}_{v}},\hat{\bm{X}_{s}})$ in Eq. (5) will vanish when it represents the continuity equation.

In this paper, to maintain numerical stability, we follow Hu et al. [14] and mix the conservative variables of small cut cells with their corresponding neighbors of larger volume fraction. The interface interaction models and more details about the conservative sharp interface methods can be found in [13, 14, 22].

2.3 Level-set method

The level-set funtion $\Phi(\bm{x},t)$ [27, 36] is used to calculate the geometrical variables such as the volume fraction, cell-face apertures, the normal vector, etc. With $|\nabla\Phi=1|$, the level-set function describes the signed distance from the interface to each cell center. The location of the interface is given by the zero level-set $\Phi=0$ implicitly and the level-set field is updated by the linear advection equation [8]

[TABLE]

where $\bm{u}$ denotes the level-set advection velocity obtained by the interface interaction model [13].

In practice, a reinitialization procedure is needed for maintaining the signed distance property of the level-set function, which will be violated during the evolution of the interface. Besides, in order to complete the interpolation stencils for the near interface cells and solve the interface interaction, we need an extending algorithm to extrapolate the fluid states.

3 The accelerated method

Considering a flow domain occupied by fluid1 and fluid2, we can obtain time step of fluid1 ${\Delta t}_{1}$ and time step of fluid2 ${\Delta t}_{2}$ according to CFL conditions. With the assumption that ${\Delta t}_{1}$ is much larger than ${\Delta t}_{2}$, the ratio of time steps is defined as $n={\Delta t}_{1}/{\Delta t}_{2}$. Note that the second-order TVD Runge-Kutta (RK2) scheme [35] is used as the time integration method in this paper.

3.1 Uniform grid

When the accelerated method is implemented on a uniform grid, it does not matter whether $n$ is an integer or not. However, when $n$ is a decimal number, the computational efficiency will be influenced if the decimal part of $n$ is not large enough. To address this issue, $n$ is truncated to an integer by reducing ${\Delta t}_{1}$ when the decimal part of $n$ is less than the preset threshold. In this paper, we choose 0.7 as the threshold for all cases. In one time cycle of the accelerated method, value of ${\Delta t}_{1}$ will be kept while value of ${\Delta t}_{2}$ will be updated along with the calculation. We assume that ${\Delta t}_{2}$ is always equal to the same value for simplicity. With the assumption that ${\Delta t}_{1}=2.8$ and ${\Delta t}_{2}=1$, the basic ideas of the accelerated method for one time cycle are illustrated in Fig. 2. Note that some procedures such as the extending step, the interaction step, the mixing step, which can be found in [14, 22], and the evolution of the level-set field are not drawn in this figure.

First, after solving the interaction step for obtaining interfacial fluxes, we compute the first step of RK2 for fluid1 with ${\Delta t}_{1}=2.8$. Both previous and future conservative variables of fluid1 are stored for interpolating its intermediate states. Second, we compute the first step of RK2 for fluid2 and the level-set field with the time step ${\Delta t}_{2}=1$. Let $t_{1}$ and $t_{2}$ denote the evolution time of fluid1 and fluid2, respectively. In one time cycle, we have $t_{1}={\Delta t}_{1}$ while $t_{2}$ is accumulated every time ${\Delta t}_{2}$ is updated according to CFL condition. The intermediate states of fluid1 corresponding to fluid2 can be obtained by

[TABLE]

Then, with flow variables of fluid1 at $t_{1}=1$ obtained by interpolation, the interaction step is solved so that we can compute the second step of RK2 for fluid2 and the level-set field. The interfacial fluxes transferred to fluid1 and fluid2 are not equal due to the separated evolution of two fluids. To maintain the conservative property, we record the interfacial fluxes transferred to fluid2 during the whole time cycle and make the interfacial flux correction for fluid1 at last. The correction quantities for both inviscid and viscous interfacial fluxes can be calculated according to two different situations as follows:

i. if the interfacial fluxes will be used to update the states of both fluid1 and fluid2

[TABLE]

ii. if the interfacial fluxes will only be used to update the states of fluid2

[TABLE]

Note that the interfacial flux correction quantities should be set to 0 at the beginning of one time cycle and the superscript $*$ denotes the previous values. The calculation of the correction quantities for interfacial fluxes due to surface tension is slightly different. Caused by mechanical equilibrium, surface tension will result in a pressure jump at the interface, which can be expressed by

[TABLE]

where $\sigma$, $\kappa$ are the surface tension coefficient and curvature, respectively. As in Luo et al. [22], the interfacial fluxes due to surface tension $\hat{\bm{X}_{s}}$ are integrated into the inviscid interfacial fluxes $\hat{\bm{X}}$, and then we have $\hat{\bm{X}_{1}}$ for fluid1 and $\hat{\bm{X}_{2}}$ for fluid2 with the relation $\frac{\hat{\bm{X}_{1}}}{\hat{\bm{X}_{2}}}=\frac{p_{I,1}}{p_{I,2}}$. The correction quantities for interfacial fluxes due to surface tension can be obtained by:

i. if the interfacial fluxes will be used to update the states of both fluid1 and fluid2

[TABLE]

ii. if the interfacial fluxes will only be used to update the states of fluid2

[TABLE]

Next, we repeat updating the states of fluid2 and the level-set field until $t_{2}$ is greater than (or equal to) $t_{1}$. At the end of one time cycle, $t_{2}$ should be equal to $t_{1}$ for synchronization. Hence, we truncate ${\Delta t}_{2}$ by ${\Delta t}_{2}=t_{1}-{t_{2}}^{*}$ and compute the first step of RK2 for fluid2 and the level-set field with it. Finally, after computing the second step of RK2 for fluid1, fluid2 and the level-set field together, we make the interfacial flux correction for fluid1.

In the original method, the minimum time step i.e. ${\Delta t}_{2}=1$ is used to evolve the fluid states of both fluid1 and fluid2. In order to advance the flow field to $t=3$, we have to compute RK2 three times for both fluid1 and fluid2. With the accelerated method, we compute RK2 three times for fluid2 and only once for fluid1 to advance the flow field to $t=2.8$. Compared with the RK2 scheme, time costs of new extra procedures such as interpolation, interfacial fluxes recording, interfacial flux correction are much lower and can almost be omitted. In the original method, let $m$ denotes the proportion of the calculation time of RK2 for fluid1 in the total calculation time of marching a time step. The speedup ratio $r_{s}$, which is defined as the ratio of previous total simulation time to modified total simulation time, can be calculated by

[TABLE]

where $[\ ]$ denotes rounding a number to the next larger integer.

With the assumption that ${\Delta t}_{1}>{\Delta t}_{2}$ , the procedures of the accelerated method for a whole time cycle can be summarized as follows:

Extend the fluid states and calculate interfacial fluxes with the interface interaction models proposed in [14, 22]. Then obtain intefacial flux correction quantities via Eq. (8) or Eq. (11) according to whether the surface tension is considered. 2. 2.

Compute the first step of RK2 for fluid1 with ${\Delta t}_{1}$ and store both previous and future conservative variables of fluid1 for interpolating its intermediate states via Eq. (7). 3. 3.

Compute the first step of RK2 for fluid2 and the level-set field with ${\Delta t}_{2}$. With the interpolated intermediate states of fluid1, extend fluid states, solve the interaction step and update interfacial flux correction quantities via Eq. (9) or Eq. (12) according to whether the surface tension is considered. 4. 4.

Compute the second step of RK2 for fluid2 and the level-set field. Then update ${\Delta t}_{2}$ according to CFL conditions while ${\Delta t}_{1}$ is kept, and accumulate $t_{2}$. 5. 5.

Repeat steps 3-4 until $t_{2}$ is greater than (or equal to) $t_{1}$. Truncate ${\Delta t}_{2}$ and implement the step 3. Note that the interfacial flux correction quantities should be updated via Eq. (8) or Eq. (11). Compute the second step of RK2 for fluid1, fluid2 and the level-set field together and finally make the interfacial flux correction for fluid1.

3.2 Multi-resolution grid

The original MR algorithm [11] employs a storage-and-operation-splitting pyramid data structure which will be adaptively updated through the multi-resolution analysis. The pyramid data structure is defined as a set of tables in which the basic element is called a node. There are two types of nodes, empty node which means there is no computational data and occupied node which contains a block of $n\times m$ computational cells. When it has no substructures, a node is called a leaf, which can be viewd as a sub-domain with a uniform grid and the simulations are carried out on. First, by using a projection operator $P_{l+1\rightarrow l}$ [6, 1], cell-averages at the level $l$ are calculated from its child-level $l+1$ exactly. Then a prediction operator $P_{l\rightarrow l+1}$ [6, 1] is employed to give an approximation of cell-averages at the level $l+1$ by using those values obtained by projection at the level $l$. The prediction errors between the approximation and the exact cell-averages at the level $l+1$ is called details. When the detail is smaller than a certain threshold, the flow filed can be represented by the values on the coarser grid, otherwise the corresponding nodes need to be refined. In order to improve the computational efficiency, the local time-stepping scheme [6] is implemented. The main principle is that the finer leaves are updated more frequently than the coarser leaves, which is handled by introducing larger time step for the evolution of the coarser leaves, e.g., the time step of leaves at the level $l$ is twice the time step of leaves at the level $l+1$. The intermediate states of the coarser leaves are obtained by interpolation. At the coarse-fine resolution boundaries, the discrete conservation is violated because the fluxes at the interface of the coarser leaves are computed by using the cell-averaged values of the adjacent nodes, which are obtained by projecting from the child-level. To ensure the discrete conservation, a flux correction step is applied to the coarser leaves located at the coarse-fine resolution boundaries.

When the accelerated method is combined with the MR algorithm, those procedures of multi-resolution analysis such as projection, prediction and refinement can be implemented without modification. Since all multiphase leaves are located at the finest level according to the MR algorithm [11], those procedures corresponding to the accelerated method such as interpolation for fluid advanced with a larger time step, interfacial flux recording and interfacial flux correction can be applied to the multiphase leaves in a straightforward way, as no resolution jumps are involved. The time steps of fluid1 and fluid2 at the level $l$ are represented by ${\Delta t}_{1}^{l}$ and ${\Delta t}_{2}^{l}$, respectively, and the assumption that ${\Delta t}_{1}^{l}$ is much larger than ${\Delta t}_{2}^{l}$ is still held. Note that the time step of each node is set at the beginning and will be kept during the whole time cycle. With the relations that ${\Delta t}_{1}^{l}=2{\Delta t}_{1}^{l+1}$ and ${\Delta t}_{2}^{l}=2{\Delta t}_{2}^{l+1}$, the ratio of time steps $n$ at every level equals to the same value. To combined with the local time-stepping scheme, $n$ will be truncated to an integer by reducing ${\Delta t}_{1}^{l}$ when it is a decimal number. At the level $l$, the states of fluid1 will be updated $2^{l}$ times while the states of fluid2 will be updated $n\times 2^{l}$ times. Note that the minimum level is $l=0$.

For simplicity, we only consider two scale levels (level 0 and level 1) here. With the assumption that ${\Delta t}_{1}^{1}=2$ and ${\Delta t}_{2}^{1}=1$, the basic ideas of the accelerated method combined with the MR algorithm for one time cycle are illustrated in Fig. 3. It can be observed that after being modified, the evolution of each fluid is almost the same as the original process except the extra procedures implemented at the finest level. As discussed in the above paragraph, these extra procedures related to the accelerated method are implemented at the finest level as same as on a uniform grid. Hence we do not discuss the repetitive details here. Let $L_{max}$ denotes the maximum level of the pyramid data structure and $I$ represents the counter of the cycle, which counts from zero. The general process of the accelerated method combined with the MR algorithm for one time cycle is given as follows:

Compute the first step of RK2 for leaves satisfying the conditions as follows:

[TABLE]

where $\%$ denotes the modulus operator. 2. 2.

Obtain intermediate states of both fluid1 and fluid2 at the coarser level and intermediate states of fluid1 at the finest level by interpolation. Record the corresponding fluxes for conservative correction. 3. 3.

Compute the second step of RK2 for leaves satisfying the conditions as follows:

[TABLE] 4. 4.

When two successive levels reach the same discrete time, a flux correction step is applied to those leaves located at the coarse-fine resolution boundaries. In addition, at the finest level, make the interfacial flux correction for fluid1 every time the evolution time of two fluids are the same. Then the counter $I$ is accumulated once. 5. 5.

Repeat steps $1-4$ until $I$ equals $2^{L_{max}}\times n$.

4 Numerical examples

The following numerical examples are provided to illustrate the potential of the accelerated sharp-interface method for compressible and weakly compressible multiphase flows. Different kinds of flow problems involving shock waves, the viscous force and surface tension are calculated to compare the accuracy and speed of the present method with the original method. For all cases, with the CFL number of 0.6, one-phase calculations are carried out with the fifth-order WENO-LF method [15] and the second-order TVD Runge-Kutta scheme [35]. In section 4.1 - section 4.2, an ideal-gas EOS is used for gas while water is modelled with Tait’s equation Eq. (4). In section 4.3 - section 4.7, for the wealy compressible model, Tait’s equation Eq. (4) is used as the artificial EOS for two phases.

4.1 Gas-water interaction

Two 1D problems with the gas-water interaction are calculated to validate the present method for compressible multiphase flows. With respect to the state of water at 1 atmosphere and length scale 1 m, the non-dimensional parameters in Tait’s equation are given by $B=3310$, $\rho_{0}=1$ and $p_{0}=1$. In the two problems, simulations are carried out on a uniform mesh with 200 grid points and the reference solution is sampled on 1000 grid points. Due to the small calculation costs of the two problems, speedup is slight and will not be discussed here.

4.1.1 Gas-water shock tube problem

The initial conditions are given as

[TABLE]

and the final time is $t=0.0008$. In this problem, the transmitted and reflected wave fronts move rapidly while the high-pressure gas expands slowly. At first, the time step of water is about three times the time step of gas according to CFL conditions. Fig. 4 gives the computed pressure, velocity, density profiles of the present method, which show a good agreement with the original results. No oscillations are observed and the position of the gas-water interface is the same as that in the original result. The numerical dissipation increases because the states of water is evolved with a larger time step. As a result, the shock wave is underestimated slightly, which can be observed from the enlarged velocity profile Fig. 4(d).

4.1.2 Shock impacting on an air-water interface

Taken from [20], this case is an air bubble collapse problem in one dimension. The initial conditions are given as

[TABLE]

and the final time is $t=0.003$. In this problem, the underwater shock wave impacts on the air-water interface so that a weak shock wave transmitted into the air and a very strong rarefaction wave reflected back to the water medium. At first, the time step of air is about four times the time step of water according to CFL conditions. In Fig. 5, the computed pressure, velocity, density profiles of the present method show a good agreement with the original results and no oscillations are observed. The air-water interface is captured at the position as same as in the original result. The shock wave is underestimated slightly because of larger numerical dissipation caused by introducing larger time step for air.

4.2 Periodic motion of a water droplet

We consider the translation of a water droplet in a periodic domain to study the convergence and the speedup of the present method. Embedded in the air, a water droplet of radius $R=0.1\ m$ is placed at the center of a $1\ m\times 1\ m$ square computational domain with periodic boundary conditions. The reference length is $l_{f}=1\ m$ and the reference velocity is $v_{f}=\sqrt{p_{f}/\rho_{f}}$, where the reference pressure and density are given by $p_{f}=1\ atm$, $\rho_{f}=1\ kg/m^{3}$, respectively. In this case, the non-dimensional parameters in Tait’s equation are chosen as $B=3310$, $\rho_{0}=1000$ and $p_{0}=1$. The initial conditions are given by

[TABLE]

and the final time is $t=10\sqrt{10}$ which is exactly a cycle.

In Table. 1, we list the speedup ratios for different resolutions. When the domain is discretized by a uniform mesh with increasing grid points from $32\times 32$ to $256\times 256$, it can be observed that the speedup ratio keeps increasing along with gird refinement. According to CFL conditions, the time step of gas is about four times the time step of water during the whole simulation in this case. As resolution growing, the ratio of the number of grid points occupied by the gas to the number of grid points occupied by the water increases so that $m$ defined in Eq. (13) increases. Hence, according to Eq. (13), the speedup ratio $r_{s}$ will increase along with grid refinement. With a uniform grid involving $128\times 128\times 128$ cells, the 3D version of this case is calculated by slightly modifying the initial conditions. The velocity component in the z-axis direction $w$ is set to 0 and the level-set field is initialized by $\Phi=-0.1+\sqrt{(x-0.5)^{2}+(y-0.5)^{2}+(z-0.5)^{2}}$. As shown in Table. 1, with the resolution as same as in two dimensions, additional speedup ratio of about 0.5 is obtained in three dimensions. The reason why the speedup ratio increases is that the ratio of the number of grids points occupied by the gas to the number of grid points occupied by the water in three dimensions is larger than it in two dimensions when the resolutions are the same.

After one cycle, the water droplet returns the original position and the exact area of water droplet is $S_{exact}=0.01\pi$ due to $\nabla\cdot\bm{u}=0$. Fig. 6(b) depicts the convergence of the relative errors $E(area)=|S_{numerical}-S_{exact}|/S_{exact}$ on the uniform mesh with grid refinement, which suggests that nearly second-order accuracy is obtained for both the present and the original method.

This case is also calculated on a multi-resolution grid where the effective resolution at the finest level is $256\times 256$ with the maximum level of resolution at $L_{max}=4$. The multi-resolution representations at different time are depicted in Fig. 6(c)-(d), where the resolutions of three levels can be observed. As shown in Table. 1, the speedup ratio is reduced by about 0.8 compared to the uniform grid situation when the resolutions are the same. There are two main reasons: First, the extra procedures related to the data structure lead to the decrease of $m$ in Eq. 13. Second, as shown in Fig. 6(c)-(d), some single-phase leaves of gas are located at coarser levels where calculations are carried out with a larger time step.

Let $M$ represents the total mass of the droplet and $\rho_{i}$ denotes the cell-averaged mass density of cell $i$. Assuming that the droplet occupys cells ranging from the first to the $N$th, the relative mass is calculated by

[TABLE]

where the superscripts 0 and n represent the initial condition and the nth time step, respectively. Fig. 6(a) depicts the relative mass of the water droplet during the computation. It can be observed that the relative mass is always equal to 1, which confirms that the present method has no conservation error.

4.3 Static drop with surface tension

In this case, we consider a static liquid drop of radius $R=0.1$ embedded in the air to test the Laplace law. With the surface tension coefficient $\sigma=36$, the initial conditions are given by

[TABLE]

The computational domain is a unit square with no-slip boundary conditions, which is discretized by a uniform grid with increasing resolution from $32\times 32$ to $256\times 256$ as well as a multi-resolution grid where the effective resolution at the finest level is $256\times 256$. Speedup ratios for different resolutions are given in Table. 2.

Because of the capillary effects, the pressure distribution of steady state should satisfy the Laplace law Eq. (10). The theoretical pressure jump can be calculated as

[TABLE]

and the relative error is obtained by

[TABLE]

Fig. 7(a) shows the steady state pressure distribution obtained at the finest resolution $256\times 256$. On both the uniform grid and the multi-resolution grid, the results of the accelerated method are almost the same as the original results. From Fig. 7(b), it can be observed that nearly second-order accuracy is obtained for both the present and the original method.

4.4 Recovering a circle shape

Taken from [34], a square drop deforming in a gas-like fluid is computed for dynamic verification of surface tension. We consider a $0.2\times 0.2$ square liquid drop placed at the center of a $1\times 1$ square computational domain with no-slip boundary conditions. With the surface tension coefficient $\sigma=73$, the initial conditions are given by

[TABLE]

Note that the level-set function given above is not a strict signed distance function so it has to be reinitialized before the simulation. In this case, the simulations are carried out on the uniform grid with increasing resolution from $32\times 32$ to $256\times 256$ and the multi-resolution grid where the effective resolution at the finest level is $256\times 256$. Due to the capillary effects, the shape of the liquid drop will evolve to a circle and the pressure of the steady state will satisfy the Laplace law Eq. (10). The flow field reaches the steady state at time $t=2.0$. In Fig. 8(a)-(d), with the same resolution $256\times 256$, the interfaces of the droplet obtained by the present method are compared with the original results at different time, which are almost identical.

Unlike Schmidmayer et al. [34], we employ a weakly compressible model in this case so that mass conservation implies volume conservation. With the initial area of the square drop $S_{initial}=0.04$, the radius of the final circle can be calculated by $R=\sqrt{\frac{0.04}{\pi}}=0.112838$. Then the theoretical pressure jump can be obtained by $[\,p\,]=\frac{\sigma}{R}=646.94566$. Fig. 8(e) depicts the steady state pressure distribution obtained by the present method, which shows a good agreement with the original results. In Fig. 8(f), relative error $E(area)=|S_{final}-S_{initial}|/S_{initial}$ is plotted versus $\Delta x$, which suggests about second-order convergence for both the present and the original method.

To repeat this case in three dimensions, the velocity component in the z-axis direction $w$ is set to 0 and the level-set field is initialized by $phi=max(|x-0.5|,|y-0.5|,|z-0.5|)-0.1$. After reinitializing the level-set field, the simulation is carried out on a uniform grid involving $128\times 128\times 128$ cells. Speedup ratios for different resolutions are listed in Table. 3. With the resolution as same as in two dimensions, additional speedup ratio of about 0.4 is obtained in three dimensions. Fig. 9 depicts the 3D interfaces of the droplet at different time. At the cross section of $z=0$, we can find that the interfaces obtained by the present method are consistent with those in the original results. Note that the interfaces in Fig. 9 are slightly different from those in Fig. 8 due to their lower resolution.

4.5 Drop oscillation

In this case, we consider a $1\times 1$ square computational domain with a ellipsoidal drop located at its center to access the prediction of capillary instabilities [21]. Surrounded by a gas-like fluid, the ellipsoidal drop is initially at rest with zero kinetic energy. Then the existence of surface tension results in the oscillatory behaviour of the liquid drop along with the translation between the potential energy and the kinetic energy. As shown in Fig. 10(a), in on cycle, the shape of the liquid interface changes from an ellipse to a circle, and then it deforms into an ellipse again. With the surface tension coefficient $\sigma=341.642$, the initial conditions are give as

[TABLE]

Note that the level-set function $\phi$ has to be reinitialized before the simulation because the elliptic equation given in Eq. (20) is not a signed distance function. This case is calculated on the uniform grid involving $128\times 128$ cells and the multi-resolution grid with the effective resolution $128\times 128$ at the finest level. The speedup ratio obtained on the uniform grid is $s=2.196$ and decreases to $s=1.864$ on the multi-resolution grid.

The shapes of the droplet at different time are given in Fig. 10(b)-(c), which show a good agreement between the results obtained by the present method and the orignal method. Fig. 10(d) depicts the evolution of the kinetic energy versus simulation time. It can be observed that the peaks of the kinetic energy obtained on the multi-resolution grid are lower than those obtained on the uniform grid due to the extra errors generated by the MR algorithm. However, as long as the simulations are carried out on the same grid, the results of the present method are almost identical to the original results. According to the Rayleigh formula which has been extended to two-phase flows [9], the oscillation frequency can be obtained by

[TABLE]

Here $\rho_{l}$ and $\rho_{g}$ are the density of the liquid drop and the surrounding gas-like fluid, respectively. In addition, $o$ is the oscillation mode and $R$ is the radius of the liquid drop at the equilibrium state. Following [28, 22], we choose $o=2$ and $R=0.15825$ so that the theoretical value of the oscillation period is $T=0.0878$. Evaluated from Fig. 10(d), on both the uniform gird and the multi-resolution grid, the oscillation period obtained by the present method is $T=0.0881$ as same as in the original results. The error is about $0.3\%$ as the same in [22].

4.6 Liquid filament contraction

Proposed by Schulkes et al. [31], the liquid filament contraction is considered to validate the present method for problems with coupled effects of viscosity and surface tension. As shown in Fig. 11, the initial shape of the filament is a cylinder with two hemispherical caps at both of its ends. In addition, the half-length of the filament is $L$ and the radius of the cylinder as well as the hemispherical is $R$. Two important dimensionless parameters for this problem are the aspect ratio of the filament $L_{0}=L/R$ and the Ohnesorge number $Oh=\mu/\sqrt{\rho\sigma R}$, which measures the relative importance of the viscous force to the surface tension force. In this case, we choose $L_{0}=9.0$ and $Oh=8.6\times 10^{-3}$ so that the filament will contract into a single drop. With the the surface tension coefficient $\sigma=72\ mNm^{-1}$, the initial conditions are given as

[TABLE]

Due to the symmetry property of this case, the computational domain is $1/4$ of the physical domain with 1.0 mm in the x-axis direction and 2.0 mm in the z-axis direction. The computations are carried out on the uniform grid with increasing resolutions from $16\times 32$ to $128\times 256$ and the multi-resolution grid where the effective resolution at the finest level is $128\times 256$. To increase the computational accuracy, a two dimensional axisymmetric model is used in this case.

Fig. 12 depicts the evolution of the contracting liquid filament at the finest resolution $128\times 256$. The filament contracts quickly and forms bulbous ends which become larger and larger during the evolution. The filament fails to pinch off its bulbous ends because the capillary forces are not strong enough to overcome the viscous force. After a series of oscillations with diminishing amplitudes, the filament finally become a sphere (not shown here). It can be observed that during the evolution the shapes of the filament obtained by the present method are consistent with the original results on both the uniform grid and the multi-resolution grid.

To calculate the 3D version of this case, we modify the initial conditions by setting the velocity component in the z-axis direction to 0 and initializing the level-set field by

[TABLE]

Involving $64\times 64\times 128$ computational cells, the computational domain is $1/8$ of the physical domain with 1.0 mm in x-axis direction, 1.0 mm in the y-axis direction and 2.0 mm in the z-axis direction. Speedup ratios for different resolutions are given in Table. 4. With the resolution as same as in two dimensions, additional speedup ratio of about 0.4 is obtained in three dimensions. The 3D interfaces of the filament obtained by two methods at different time are shown in Fig. 13. For a better comparison, at the cross section of $y=0$, we depict the interfaces of the filament obtained by the accelerated method and the original method, which are well matched. Note that the interfaces in Fig. 13 are slightly different from those in Fig. 12 because their resolutions are lower.

4.7 Equilibrium shape of a water drop resting on a wall

By employing a curvature boundary condition proposed by Luo et al. [23], we consider an equilibrium water drop on a horizontal wall to validate the present method for problems with a static contact angle. As illustrated in Fig. 14(a), the initial shape of the water drop is a semi-circle with the radius $R_{0}=0.5$, which initially rests on the wall with the contact angle $90^{\circ}$. Due the the difference between the transient contact angle and the static contact angle, the water drop will deform and its final shape is determined by two parameters, the static contact angle and the $E\ddot{o}tv\ddot{o}s$ number $Eo=\rho gR_{0}^{2}/\sigma$. In this case, we ignore the gravity so that $Eo=0$ and the equilibrium shape of the water drop is a circular cap determined by the static angle only. The computational domain is a $2.5\times 2.5$ square with no-slip boundary conditions. With the surface tension coefficient $\sigma=36$, the initial conditions are given by

[TABLE]

With different static contact angles, this case is calculated on the uniform grid involving $128\times 128$ cells and the multi-resolution grid where the effective resolution at the finest level is $128\times 128$. When the water drop reaches the steady state, the wetting length can be calculated by

[TABLE]

according to the volume conservation. Fig. 14(b)-(c) show the equilibrium shapes of the water drop with different static contact angles. It can be observed that the interfaces obtained by the accelerated method are almost identical to those obtained by the original method. The numerical results of the wetting length are compared with the theoretical values in Fig. 14(d). It is found that the simulation results of both the present and the original method agree with the theoretical values quite well. Speedup ratios for different static contact angles are listed in Table. 5.

5 Concluding remarks

In this paper, we have developed an accelerated sharp-interface method for compressible and weakly compressible multiphase flows. This method inherits the simplicity of the original GFM-like methods of Hu et al. [14] and Luo et al. [22]. Different from the semi-implicit methods [18, 16, 42] and the free-surface conditions [4, 17], the present method will not generate extra oscillations and makes no assumptions for the gas phase (or the less viscous and dense phase). The conservative property of the original methods is violated due to the separated evolution of two fluids, which is handled by the interfacial flux correction. Furthermore, the present method offers a fairly direct way to be combined with the wavelet-based adaptive multi-resolution algorithm [11]. A number of numerical examples are calculated on both the uniform grid and the multi-resolution grid to validate the accelerated method. It is shown that the simulations with the present method is in good agreement with the original results. The higher the grid resolution is, the greater the acceleration performance we obtain. When the number of grid points is around $10000$, for all cases, the computational time of our accelerated method is only about $1/2$ of the original one. In future work we plan to extend the accelerated method to handle the solid-fluid coupling problems as the time step of fluid is usually much larger than the time step of solid.