Implementation of BVD (boundary variation diminishing) algorithm in simulations of compressible multiphase flows

Implementation of BVD (boundary variation diminishing) algorithm in simulations of compressible multiphase flows


We present in this work a new reconstruction scheme, so-called MUSCL-THINC-BVD scheme, to solve the five-equation model for interfacial two phase flows. This scheme employs the traditional shock capturing MUSCL (Monotone Upstream-centered Schemes for Conservation Law) scheme as well as the interface sharpening THINC (Tangent of Hyperbola for INterface Capturing) scheme as two building-blocks of spatial reconstruction using the BVD (boundary variation diminishing) principle that minimizes the variations (jumps) of the reconstructed variables at cell boundaries, and thus effectively reduces the numerical dissipations in numerical solutions. The MUSCL-THINC-BVD scheme is implemented to all state variables and volume fraction, which realizes the consistency among volume fraction and other physical variables. Benchmark tests are carried out to verify the capability of the present method in capturing the material interface as a well-defined sharp jump in volume fraction, as well as significant improvement in solution quality. The proposed scheme is a simple and effective method of practical significance for simulating compressible interfacial multiphase flows.

Compressible multiphase flows, five-equation model, interface capturing, THINC reconstruction, BVD algorithm

[cor]Corresponding author: Dr. F. Xiao (Email:

1 Introduction

Compressible multiphase flow is one of active and challenging research areas of great importance in both theoretical studies and industrial applications. For example, shock/interface interactions are thought to be crucial to the instability and evolution of material interfaces that separate different fluids as can be observed in a wide spectral of phenomena(2). The material interfaces greatly complicate the physics and make problems formidably difficult for analytical and experimental approaches in general. In many cases, numerical simulation turns out to be the most effective approach to provide quantitative information to elucidate the fundamental mechanisms behind the complex phenomena of multiphase flows.

In comparison to computation of single phase flow, development of numerical methods for multiphase flow faces more challenging tasks. The major complexity comes from the moving interfaces between different fluids that usually associate with strong discontinuities, singular forces and phase changes. Given the numerical methods developed for multiphase incompressible flow with interfaces having been reached a relatively mature stage, the numerical solvers for compressible interfacial multiphase flow are apparently insufficient. For incompressible multiphase flows with moving interfaces where the density and other physical properties, e.g. viscosity and thermal conductivity, are constant in each fluid, the one-fluid model (1) can be implemented in a straightforward manner with an assumption that the physical fields change monotonically across the interface region. So, provided an indication function which identifies the moving interface, one can uniquely determine the physical property fields for the whole computational domain. Some indication functions, such as volume of fluid (VOF) function (3); (4); (5) and level set function (6); (7); (8), have been proposed and proved to be able to well define the moving interface with compact thickness and geometrical faithfulness if solved by advanced numerical algorithms. However, substantial barrier exists when implementing the one-fluid model to compressible interfacial multiphase flow.

The new difficulties we face when applying the one-fluid model 2 to compressible interfacial multiphase flow lie in two aspects:

(I) Density and energy in compressible flow have to be computed separately in addition to the indication function, hence special formulations are required to reach a balanced state among all variables for the interface cell where a well-defined interface falls in;

(II) The numerical dissipation in the so-called high-resolution schemes designed for solving single phase compressible flow involving shock waves tends to smear out discontinuities in numerical solutions including the material interfaces, which is fatal to simulations of interfacial multiphase flows even if the schemes can produce acceptable results in single phase cases.

For issue (I) mentioned above, mixing or averaging models that consist of Euler or Navier-Stokes equations along with interface-indication function equations for each of fluid components have been derived and widely used as efficient approximations to the state of the interface cell where two or more species co-exist. A simple single-fluid model was reported in (9); (11); (10) for interfacial multiphase compressible flows using either explicit time marching or semi-implicit pressure-projection solution procedure. The latter results in a unified formulation for solving both compressible and incompressible multiphase flows. As the primitive variables are used in these models, the conservation properties are not guaranteed, and thus might not be suitable for high-Mach flows involving shock waves. Conservative formulations, which have been well-established for single phase compressible flows with shock waves, however may lead to spurious oscillations in pressure or other thermal fields (12); (13). It was found that special treatments are required in transporting the material interface and mixing/averaging the state variables to find the mixed state of fluids in the interfacial cell that satisfies pressure balance across material interface for multiple polytropic and stiff gases (14); (15); (16); (17); (18), and for van der Waals and Mie-Grüneisen equations of state (EOS) (46); (45). A more general five-equation model (43) was developed for a wide range fluids. These models apply to multiphase compressible flows with either spread interfaces or sharp interfaces. The five-equation model will be used in the present work as the PDE (partial differential equation) set to solve.

Provided the SEF models with some desired properties, such as hyperbolicity, conservation and well-balanced mixing closure without spurious oscillations in thermal variables, we can in principle implement numerical methods for single phase compressible flow (e.g. standard shock-capturing schemes) to solve multiphase ones. TVD (Total Variation Diminishing) schemes, such as the MUSCL (Monotone Upstream-centered Schemes for Conservation Law) scheme (19), can resolve discontinuities without numerical oscillations, which is of paramount importance to ensure the physical fields to be bounded and monotonic in the transition region. However, TVD schemes suffer from excessive numerical dissipation, which brings the problem (II) listed above to us. The intrinsic numerical dissipation smears out the flow structures including the discontinuities in mass fraction or volume fraction that are used to represent the material interfaces. Consequently, material interfaces are continuously blurred and spread out, which is not acceptable in many applications, especially for the simulations that need long-term computation. Applying high order schemes like WENO (Weighted Essentially Non-Oscillatory) scheme (20) to solve compressible multiphase flow are found in the literatures (21); (22). However, implementing high order schemes might generate numerical oscillations for compressible multiphase flow with more complex EOS, such as the Mie-Grüneisen equation of state. In (21), to reduce numerical oscillation introduced by high order schemes, the state variables have to be cast into characteristic fields. Although with this effort, stability cannot be guaranteed in the long-term computation even using a forbiddingly small of time step. In a recent work (23), to further reduce the numerical oscillations and to deal with complicated EOS, an approximate intermediate state at each cell edges is obtained in a more careful way to conduct characteristic decomposition. Furthermore, high order monotonicity-preserving scheme (24) was used to ensure volume fraction remain bounded. In general, the implementation of high order shock capturing schemes in compressible multiphase flows will increase complexity of algorithm and may invoke computational instability.

In order to keep material interfaces being a compact thickness during computation, special treatments are required to sharpen or steepen the interfaces. The existing methods for this purpose can be categorized into interface-tracking and interface-capturing. Interface-tracking methods like Arbitrary Lagrangian-Eulerian (ALE) (35); (36), free-Lagrange (37), front-tracking (38) and level set/ghost fluid (39); (40); (41); (42) have the distinct advantage of treating interface as a sharp discontinuity. For example, ALE and free-Lagrange methods treat interfaces as boundaries of distorted computational grids, which however makes them complicated and computationally expensive when there are large interface deformations and topological changes. Another drawback is that these methods are typically not numerically conservative at material interfaces, which may lead to a wrong prediction about the position of interfaces and shock waves.

Interface-capturing methods resolve the interfaces on fixed Eulerian grids, using special numerical techniques to sharpen the interface from spreading out. For example, in (26); (28); (25) the advection equation of the interface function is treated by artificial compression method. As a post-processing approach, anti-diffusion techniques have been introduced by (59) and (30). Another approach is to reconstruct the volume fraction under the finite volume framework by THINC (Tangent of Hyperbola for INterface Capturing) function (51). By virtue of the desirable characteristics of the hyperbolic tangent function in mimicking the jump-like profile of the volume fraction field, the sharp interface can be accurately captured in a simple way (29); (23). However, unlike incompressible multiphase flow, a common occurrence when applying various interface-sharping methods explicitly on the SEF models is that velocity and pressure oscillations may occur across the interface (28); (29); (30); (26); (27) due to the inconsistency between the physical variables and the sharpened or compressed volume fraction field. As stated in (26); (28), in contrast to incompressible flows where density of fluid is fixed, artificial interface sharpening scheme cannot be applied alone to volume fraction function in compressible multiphase cases. In compressible multiphase flows, neither fluid density nor volume fraction alone is sufficient to determine the interface location and the fluid density. As some remedies, density correction equations are formulated in (26); (28). In (29), a homogeneous reconstruction has been proposed where the reconstructed volume fraction is used to extrapolate the remaining conservative variables across the interface to ensure the thermal and mechanical consistency across the isolated material interfaces. In (27), a consistent compression method has been discussed to maintain equilibrium across interfaces.

To alleviate the defect of shock capturing scheme when solving single-fluid model for compressible multiphase flows, a novel spatial reconstruction is presented in this work to resolve contact discontinuities including material interfaces with substantially reduced numerical dissipation, which then maintains the sharpness of the transition layer of material interfaces throughout even long term computations. The scheme, so-called MUSCL-THINC-BVD, implements the boundary variation diminishing (BVD) algorithm (32); (33) with the traditional MUSCL scheme and the interface-sharpening THINC scheme as two building-blocks for reconstruction. The BVD algorithm choose a reconstruction function between MUSCL and THINC, so as to minimize the variations (jumps) of the reconstructed variables at cell boundaries, which in turn effectively removes the numerical dissipations in numerical solutions. More importantly, we apply MUSCL-THINC-BVD scheme to all state variables and volume fraction, so sound consistency is achieved among volume fraction and other physical variables. Resultantly, the manipulations to the physical variables according to the volume fraction in other existing methods are not needed in the present method. The numerical model is formulated under a standard finite volume framework with a Riemann solver in the wave propagation form (34). The numerical results of benchmark tests verify the capability of the present method in capturing the material interface as a well-defined sharp jump in volume fraction, as well as significant improvement in solution quality.

The format of this paper is outlined as follows. In Section 2, the governing equations of the five-equation model and closure strategies are stated. In Section 3, after a brief review of the finite volume method in wave-propagation form for solving the quasi-conservative five-equation model, the details of the new MUSCL-THINC-BVD scheme for spatial reconstruction are presented. In Section 4, numerical results of benchmark tests are presented in comparison with other high-order methods. Some concluding remarks end the paper in Section 5.

2 Computational models

2.1 Governing equations

In this work, the inviscid compressible two-component flows are formulated by the five-equation model developed in (43). By assuming that the material interface is in mechanical equilibrium of mixed pressure and velocity, the five-equation model consists of two continuity equations for phasic mass, a momentum equation, an energy equation and an advection equation of volume fraction as follows


where and denote in turn the th phasic density and volume fraction for , the vector of particle velocity, the mixture pressure and the total energy. When considering more than two-phases, the five-equation model can be extended by supplementing additional continuity equations and volume fraction advection equations for each new phase.

2.2 Closures strategy

To close the system, the fluid of each phase is assumed to satisfy the Mie-Grüneisen equation of state,


where is the Grüneisen coefficient, and , are the properly chosen states of the pressure and internal energy along some reference curves (e.g., along an isentrope or other empirically fitting curves) in order to match the experimental data of the material examined (44). Usually, parameters , and can be taken as functions only of the density. This equation of state can be employed to approximate a wide variety of materials including some gaesous or solid explosives and solid metals under high pressure.

The further closure of the system is completed by defining the mixed volume fraction, density and internal energy as


Derived in (45), under the isobaric assumption the mixture Grüneisen coefficient and , can be expressed as


The mixture pressure is then calculated by


It should be noted that the mixing rule of Eq.(4) and Eq.(5) ensure that the mixed pressure is free of spurious osclillations, which is particularly important to prevent the spurious pressure oscillations across the material interfaces (14); (18); (16); (45); (46).

3 Numerical methods

For the sake of simplicity, we introduce the numerical method in one dimension. Our numerical method can be extended to the multidimensions on structured grids directly in dimension-wise reconstruction fashion. We will first review the finite volume method in the wave propagation form (34) used in this work and then give details about the new MUSCL-THINC-BVD reconstruction scheme.

3.1 Wave propagation method

We rewrite the one dimensional quasi-conservative five-equation model (1) as


where the vectors of physical variables q and flux functions f are

q (7)

respectively. The matrix is defined as


where denotes the velocity component in direction.

We divide the computational domain into non-overlapping cell elements, , , with a uniform grid with the spacing . For a standard finite volume method, the volume-integrated average value in the cell is defined as


Denoting all the spatial discretization terms in (6) by , the semi-discrete version of the finite volume formulation can be expressed as a system of ordinary differential equations (ODEs)


In the wave-propagation method, the spatial discretization for cell is computed by


where and , are the right- and left-moving fluctuations, respectively, which enter into the grid cell, and is the total fluctuation within . We need to solve Riemann problems to determine these fluctuations. The right- and left-moving fluctuations can be calculated by


where moving speeds and the jumps () of three propagating discontinuities can be solved by Riemann solver (48) with the reconstructed values and computed from the reconstruction functions and to the left and right sides of cell edge , respectively. Similarly, the total fluctuation can be determined by


We will describe with details about the reconstructions to get these values, and , at cell boundaries in next subsection as the core part of this paper.

In practice, given the reconstructed values and , the minimum and maximum moving speeds and can be estimated by HLLC Riemann solver (48) as


where and are sound speeds calculated by reconstructed values and respectively. Then the speed of the middle wave is estimated by


The left-side intermediate state variables is evaluated by


where the vector , and the intermediate pressure may be estimated as


Analogously, the right-side intermediate state variables is


Then we calculate the jumps as


Given the spatial discretization, we employ three-stage third-order SSP (Strong Stability-Preserving) Runge-Kutta scheme (49)


to solve the time evolution ODEs, where and denote the intermediate values at the sub-steps.

3.2 MUSCL-THINC-BVD reconstruction

In the previous subsection, we left the boundary values, and , to be determined, which are presented in this subsection. We denote any single variable for reconstruction by , which can be primitive variable, conservative variable or characteristic variable.

The values and at cell boundaries are computed from the piecewise reconstruction functions in cell . In the present work, we designed the MUSCL-THINC-BVD reconstruction scheme to capture both smooth and nonsmooth solutions. The BVD algorithm makes use of the MUSCL scheme (19) and the THINC scheme (51) as the candidates for spatial reconstruction.

In the MUSCL scheme, a piecewise linear function is constructed from the volume-integrated average values , which reads


where and is the slope defined at the cell center . To prevent numerical oscillation, a slope limiter (19); (50) is used to get numerical solutions satisfying the TVD property. We denote the reconstructed value at cell boundaries from MUSCL reconstruction as and . The MUSCL scheme, in spite of popular use in various numerical models, has excessive numerical dissipation and tends to smear out flow structures, which might be a fatal drawback in simulating interfacial multiphase flows.

Being another reconstruction candidate, the THINC (31); (51) uses the hyperbolic tangent function, which is a differentiable and monotone function that fits well a step-like discontinuity. The THINC reconstruction function is written as


where , and . The jump thickness is controlled by parameter . In our numerical tests shown later a constant value of is used. The unknown , which represents the location of the jump center, is computed from . Then the reconstructed values at cell boundaries by THINC function can be expressed by


where and , where and is a mapping factor to project the physical fields onto .

The final effective reconstruction function is determined by the BVD algorithm (32); (33), which choose the reconstruction function between and so that the variations of the reconstructed values at cell boundaries are minimized. BVD algorithm prefers the THINC reconstruction within a cell where a discontinuity exists. It is sensible that the THINC reconstruction should only be employed when a discontinuity is detected. In practice, a cell where a discontinuity may exist can be identified by the following conditions


where is a small positive (e.g.,).

In summary, the effective reconstruction function of MUSCL-THINC-BVD scheme reads


where the minimum value of total boundary variation (TBV) , for reconstruction function , is defined as

Figure 1: Illustration of one possible situation corresponding to when calculating .

Thus, THINC reconstruction function will be employed in the targeted cell if the minimum TBV value of THINC is smaller than that of MUSCL. In Fig. 1, we illustrate one possible situation corresponding to when evaluating . As stated in (32), the BVD algorithm will realize the polynomial interpolation for smooth solution while for discontinuous solution a step like function will be preferred.

As shown in numerical tests in this paper, discontinuities including material interface can be resolved by the MUSCL-THINC-BVD scheme with substantially reduced numerical dissipation in comparison with other existing methods. The material interface can be captured sharply and any extra step, like anti-diffusion or other artificial interface sharpening techniques used in the existing works (23); (29); (26); (28), is not needed here. More importantly, the MUSCL-THINC-BVD scheme is applied not only the volume fraction but also to other all physical variables, which automatically leads to the consistent reconstructions among the physical fields. As observed in our numerical results, no suprious numerical oscillation is generated in vicinity of material interfaces. It is usualy not trival to other anti-diffusion or artificial compression methods aforementioned. For example, in (27); (28); (29) anti-diffusion post-processing steps are required to adjust other state variables across interfaces to get around the oscillations.

As discussed in (52); (53); (21); (22), when high-order reconstructions, such as MUSCL or WENO shemes, are applied, special attention must be paid to decide which physical variables should be reconstructed. It is concluded that one should implement high-order reconstructions to primitive variables or characteristic variables to prevent numerical oscillations in velocity and pressure across material interfaces. However, reconstructing conservative variables or flux functions may cause numerical oscillations. We show in the rest part of this section that THINC reconstruction ensures the consistency among the reconstructed variables across material interface even if the reconstruction is conducted for the conservative variables.

We consider one-dimensional interface only problem where initial condition consists of constant velocity , uniform pressure and constant phasic densities and . Across a material interface, other variables, such as mixture densities , mass fraction and volume fraction have jumps. Without loss of generality, a positive velocity is considered here. Then the fluctuations for cell can be calculated as


By using Eq. (11) and Euler one-step forward time scheme, the cell average can be updated by


or a component form,


We denote the reconstruction operator to compute by .

As shown below, when we implement the reconstructions to the conservative variables, spurious oscillations may be generated in velocity . To facilitate discussions, we address the consistency in velocity by

Definition 1. The reconstruction is -consistent if the numerical results from (29) satisfy for isolated material interface where velocity and pressure are uniform.

In the present model, in order to calculate velocity , we first compute density from the two conservation equations for phasic densities,


From the momentum equation, we have


where we assume the operator should satisfy in which represents the reconstructed variable, is constant, then the velocity is retrieved by


It is clear that the -consistent condition requires


to maintain . I.e. the reconstruction operator should be consistent with and at cell faces to ensure .

Concerning -consistent property, we have

Proposition 1. All schemes satisfy -consistent condition if reconstruction are conducted in terms of primitive variables.

Proof. It is straightforward that since the reconstruction is conducted with , which has also been discussed in (14); (22).

Proposition 2. Piecewise constant reconstruction scheme for conservative variables satisfy -consistent condition.


For piecewise constant reconstruction, it is obvious that above condition (33) can be satisfied, as and .

Proposition 3. All linear reconstruction schemes for conservative variables satisfy -consistent condition.


For general linear schemes, the reconstructed values at cell boundaries can be expressed as linear combinations of the cell average values, , on the stencils, . That is


where the coefficients are the same for , and . With the conclusion in proposition 2 given above for piecewise constant reconstruction, we know that , which leads to from Eq. (32).

Proposition 4. Non-linear schemes for conservative variables may not satisfy -consistent condition.

Proof. Taking 5th order WENO scheme as an example, the reconstructed values at cell faces are a combination of three third order linear schemes through nonlinear weights. The reconstructed value can be expressed as


where is the third order approximation on stencil , on and on . is the nonlinear weight corresponding to the . To satisfy -consistent condition, it is necessary that


Because the nonlinear weights are separately determined according the smoothness of the reconstructed variables, there is no guarantee that Eq. (36) always holds. Thus, WENO scheme is not -consistent when applied to conservative variables, which has also been reported in (52); (53). This observation applies to all high-resolution schemes using nonlinear weights to suppress numerical oscillations.

Proposition 5. THINC reconstruction satisfies -consistent condition.

Proof. We consider a material interface in cell with volume fraction , , which divides the cell into two regions sharing uniform physical properties with the neighboring cells respectively. Then, we have . Without loss of generality, we assume and , the reconstructed values of phasic densities at cell face are


Recall (23), we get and , which leads to and .

From the momentum equation, we have


Again we can find . Then, we finally get