Implementation of BVD (boundary variation diminishing) algorithm in simulations of compressible multiphase flows
Abstract
We present in this work a new reconstruction scheme, socalled MUSCLTHINCBVD scheme, to solve the fiveequation model for interfacial two phase flows. This scheme employs the traditional shock capturing MUSCL (Monotone Upstreamcentered Schemes for Conservation Law) scheme as well as the interface sharpening THINC (Tangent of Hyperbola for INterface Capturing) scheme as two buildingblocks 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 MUSCLTHINCBVD 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 welldefined 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.
keywords:
Compressible multiphase flows, fiveequation model, interface capturing, THINC reconstruction, BVD algorithm[cor]Corresponding author: Dr. F. Xiao (Email: xiao@es.titech.ac.jp)
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 onefluid 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 onefluid model to compressible interfacial multiphase flow.
The new difficulties we face when applying the onefluid model

(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 welldefined interface falls in;

(II) The numerical dissipation in the socalled highresolution 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 NavierStokes equations along with interfaceindication 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 coexist. A simple singlefluid model was reported in (9); (11); (10) for interfacial multiphase compressible flows using either explicit time marching or semiimplicit pressureprojection 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 highMach flows involving shock waves. Conservative formulations, which have been wellestablished 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 MieGrüneisen equations of state (EOS) (46); (45). A more general fiveequation model (43) was developed for a wide range fluids. These models apply to multiphase compressible flows with either spread interfaces or sharp interfaces. The fiveequation 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 wellbalanced mixing closure without spurious oscillations in thermal variables, we can in principle implement numerical methods for single phase compressible flow (e.g. standard shockcapturing schemes) to solve multiphase ones. TVD (Total Variation Diminishing) schemes, such as the MUSCL (Monotone Upstreamcentered 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 longterm computation. Applying high order schemes like WENO (Weighted Essentially NonOscillatory) 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 MieGrü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 longterm 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 monotonicitypreserving 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 interfacetracking and interfacecapturing. Interfacetracking methods like Arbitrary LagrangianEulerian (ALE) (35); (36), freeLagrange (37), fronttracking (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 freeLagrange 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.
Interfacecapturing 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 postprocessing approach, antidiffusion 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 jumplike 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 interfacesharping 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 singlefluid 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, socalled MUSCLTHINCBVD, implements the boundary variation diminishing (BVD) algorithm (32); (33) with the traditional MUSCL scheme and the interfacesharpening THINC scheme as two buildingblocks 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 MUSCLTHINCBVD 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 welldefined 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 fiveequation model and closure strategies are stated. In Section 3, after a brief review of the finite volume method in wavepropagation form for solving the quasiconservative fiveequation model, the details of the new MUSCLTHINCBVD scheme for spatial reconstruction are presented. In Section 4, numerical results of benchmark tests are presented in comparison with other highorder methods. Some concluding remarks end the paper in Section 5.
2 Computational models
2.1 Governing equations
In this work, the inviscid compressible twocomponent flows are formulated by the fiveequation model developed in (43). By assuming that the material interface is in mechanical equilibrium of mixed pressure and velocity, the fiveequation model consists of two continuity equations for phasic mass, a momentum equation, an energy equation and an advection equation of volume fraction as follows
(1)  
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 twophases, the fiveequation 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 MieGrüneisen equation of state,
(2) 
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
(3)  
Derived in (45), under the isobaric assumption the mixture Grüneisen coefficient and , can be expressed as
(4)  
The mixture pressure is then calculated by
(5) 
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 dimensionwise 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 MUSCLTHINCBVD reconstruction scheme.
3.1 Wave propagation method
We rewrite the one dimensional quasiconservative fiveequation model (1) as
(6) 
where the vectors of physical variables q and flux functions f are
q  (7)  
f 
respectively. The matrix is defined as
(8) 
where denotes the velocity component in direction.
We divide the computational domain into nonoverlapping cell elements, , , with a uniform grid with the spacing . For a standard finite volume method, the volumeintegrated average value in the cell is defined as
(9) 
Denoting all the spatial discretization terms in (6) by , the semidiscrete version of the finite volume formulation can be expressed as a system of ordinary differential equations (ODEs)
(10) 
In the wavepropagation method, the spatial discretization for cell is computed by
(11) 
where and , are the right and leftmoving 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 leftmoving fluctuations can be calculated by
(12) 
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
(13) 
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
(14)  
where and are sound speeds calculated by reconstructed values and respectively. Then the speed of the middle wave is estimated by
(15) 
The leftside intermediate state variables is evaluated by
(16) 
where the vector , and the intermediate pressure may be estimated as
(17) 
Analogously, the rightside intermediate state variables is
(18) 
Then we calculate the jumps as
(19)  
Given the spatial discretization, we employ threestage thirdorder SSP (Strong StabilityPreserving) RungeKutta scheme (49)
(20)  
to solve the time evolution ODEs, where and denote the intermediate values at the substeps.
3.2 MUSCLTHINCBVD 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 MUSCLTHINCBVD 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 volumeintegrated average values , which reads
(21) 
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 steplike discontinuity. The THINC reconstruction function is written as
(22) 
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
(23)  
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
(24)  
where is a small positive (e.g.,).
In summary, the effective reconstruction function of MUSCLTHINCBVD scheme reads
(25) 
where the minimum value of total boundary variation (TBV) , for reconstruction function , is defined as
(26)  
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 MUSCLTHINCBVD scheme with substantially reduced numerical dissipation in comparison with other existing methods. The material interface can be captured sharply and any extra step, like antidiffusion or other artificial interface sharpening techniques used in the existing works (23); (29); (26); (28), is not needed here. More importantly, the MUSCLTHINCBVD 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 antidiffusion or artificial compression methods aforementioned. For example, in (27); (28); (29) antidiffusion postprocessing steps are required to adjust other state variables across interfaces to get around the oscillations.
As discussed in (52); (53); (21); (22), when highorder 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 highorder 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 onedimensional 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
(27)  
By using Eq. (11) and Euler onestep forward time scheme, the cell average can be updated by
(28) 
or a component form,
(29) 
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,
(30) 
From the momentum equation, we have
(31) 
where we assume the operator should satisfy in which represents the reconstructed variable, is constant, then the velocity is retrieved by
(32) 
It is clear that the consistent condition requires
(33) 
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.
Proof.
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.
Proof.
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
(34) 
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. Nonlinear 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
(35) 
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
(36) 
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 highresolution 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
(37)  
Recall (23), we get and , which leads to and .
From the momentum equation, we have
(38) 
Again we can find . Then, we finally get