A robust and accurate approach to computing compressible multiphase flow: Stratified flow model and AUSM⁺-up scheme

Abstract

In this paper, we propose a new approach to compute compressible multifluid equations. Firstly, a single-pressure compressible multifluid model based on the stratified flow model is proposed. The stratified flow model, which defines different fluids in separated regions, is shown to be amenable to the finite volume method. We can apply the conservation law to each subregion and obtain a set of balance equations¹. Secondly, the AUSM⁺ scheme, which is originally designed for the compressible gas flow, is extended to solve compressible liquid flows. By introducing additional dissipation terms into the numerical flux, the new scheme, called AUSM⁺-up, can be applied to both liquid and gas flows. Thirdly, the contribution to the numerical flux due to interactions between different phases is taken into account and solved by the exact Riemann solver. We will show that the proposed approach yields an accurate and robust method for computing compressible multiphase flows involving discontinuities, such as shock waves and fluid interfaces. Several one-dimensional test problems are used to demonstrate the capability of our method, including the Ransom’s water faucet problem and the air–water shock tube problem. Finally, several two dimensional problems will show the capability to capture enormous details and complicated wave patterns in flows having large disparities in the fluid density and velocities, such as interactions between water shock wave and air bubble, between air shock wave and water column(s), and underwater explosion.

Introduction

Seeking an accurate method to simulate compressible multifluid and multiphase flows has been an important research topic for many engineering applications. Examples of compressible multifluid flows include the cooling system in the conventional nuclear reactor, the fuel transport system in which the fuel and gas are transported simultaneously, and the generation, deformation and collapse of cavities around the underwater propeller or within the high pressure injector. The flow phenomena involved in these systems are very complicated, entailing reliable mathematical description (modeling) and numerical methods.

In this paper, we will concentrate on the numerical algorithm for compressible multifluid equations in which the fluids are assumed inter-penetrating, non-homogeneous and non-equilibrium; that is, each fluid has its own velocity and temperature fields at the same location, but all fluids share the same pressure. As shown by Buyevich [1] and Ishii [2], one can use two sets of Euler/Navier–Stoke equations to describe the motion of gas and liquid phase fluids respectively, as given in the following:

$$\begin{array}{cc}\hfill & \frac{\mathrm{\partial }({\alpha }_i{\rho }_i)}{\mathrm{\partial }t}+\nabla ·({\alpha }_i{\rho }_i{\overrightarrow{\text{v}}}_i)=S_{\rho }\text{,}\hfill \\ \hfill & \frac{\mathrm{\partial }({\alpha }_i{\rho }_i{\overrightarrow{\text{v}}}_i)}{\mathrm{\partial }t}+\nabla ·({\alpha }_i{\rho }_i{\overrightarrow{\text{v}}}_i{\overrightarrow{\text{v}}}_i)+\nabla ({\alpha }_{i}p)=p\nabla {\alpha }_i+{\overrightarrow{S}}_v\text{,}\hfill \\ \hfill & \frac{\mathrm{\partial }({\alpha }_i{\rho }_i{E}_i)}{\mathrm{\partial }t}+\nabla ·({\alpha }_i{\rho }_i{H}_i{\overrightarrow{\text{v}}}_i)=-p\frac{\mathrm{\partial }{\alpha }_i}{\mathrm{\partial }t}+S_e\text{,}\hfill \end{array}$$

where the subscript “i”=“g” or “l”, representing gas or liquid phase fluid respectively. α_i is the void fraction of fluid “i”, and must satisfy the constraint, ${\alpha }_{\text{g}}+{\alpha }_{\text{l}}=1$. The RHS of Eq. (1) represents the interactions that couple both fluid motions together. The “S” terms on the right hand side represent a group of terms that arise from interfacial physics, viscous effects, phase change, body forces, etc. They are expressed in differential (first or higher derivatives) or non-differential form. The system of equations is quite general and has been widely used to describe multiphase flows in a variety of applications. However, the inviscid limit of the multifluid model has been known to be problematic, giving rise to instability, loss of accuracy, and non-convergence in numerical solution. It is primarily attributable to the fact that the system can become non-hyperbolic and ill-posed. Hence, the root of the problems points to the first derivative terms that differ from that for the single fluid equations (which are known to be hyperbolic), namely, p∇α_i and $-p\frac{\mathrm{\partial }{\alpha }_i}{\mathrm{\partial }t}$. Additionally, these terms are not in conservative form. It is not clear that a discontinuity in the sense of weak solution still remains valid. Some comprehensive reviews of the compressible multifluid models can be found in [3], [4].

It has been recognized that non-hyperbolicity is a major reason for causing numerical instability. Thus, a great deal of efforts have been focused on how to improve the hyperbolicity of the system and make the multifluid equations well-posed; clearly altering first derivative terms with some physical basis is necessary, such as interfacial pressure correction [5], virtual mass [6], or separate pressures [3], [7]. However, we will show that non-hyperbolicity is not the only reason for causing numerical instability. It is as important to properly handle the discretization and define numerical fluxes that includes all relevant interactions terms between the same and different phases in order to obtain a stable and accurate numerical solution.

To highlight the importance of a numerical method for solving a multifluid model, we consider the case of an 1D moving fluid interface with constant velocity and pressure. Since all convection fluxes cancel out with each other, Eq. (1) reduces to a simple equation that the numerical solution must satisfy the so called pressure non-disturbing condition [7] given as follows:

$$\nabla ({\alpha }_{i}p)=p\nabla {\alpha }_i\text{.}$$

The LHS of the above equation is in conservative form and similarly exists in the single fluid equations, hence it is rather clear what to do about its discretization. However, it is not so obvious how to discretize the RHS term, ∇α_i, which is in nonconservative form, because it will have to be compatible with how the pressure flux α_ip is evaluated in order to satisfy the above equation. For example, a central differencing of ∇α_i is unlikely to be compatible with an upwind differencing of ∇(α_ip). On the other hand if the discretization of the LHS requires analytical form of the eigenstructure of Eq. (1), which unfortunately is not available in general, then it will be difficult to come up with a compatible RHS.

While Eq. (1) is usually derived from an averaging procedure [1], [2], Stewart and Wendroff [3] gave another approach for derivation based on the stratified flow model. Adopting the concept of the stratified flow model has several advantages. Firstly, it gives a clear view of the mathematical representation of physics involved in multifluid flow. Secondly, we find that it provides a clue as to the construction of numerical fluxes. It has led us to recognize various types of interactions, not only that occurring within a cell between different phases, but also that at the cell boundaries, as illustrated in Fig. 1. The former is the in-cell interaction that gives rise to the nonconservative term in Eq. (2), marked by A. The latter, marked by B in Fig. 1b, is a natural consequence of the finite volume method and is consistent with the stratified model; it is however not observed previously in the literature. We refer these interactions due to different phases collectively to as inter-phasic terms. We first explore the stratified flow model in our earlier paper [8] and further refine the numerical procedure in [9] by including the Riemann solver and applications to various problems. Independently, the idea of recognizing the presence of different phases in each sub-volume appears in [10]. Their method is further employed for different studies [11], [12] in the multifluid framework.

To our knowledge, all numerical methods for multifluid flows are extended from the ones for single fluid. The extension is not necessarily straightforward; in fact, difficulties arise because we need be concerned with additional issues. These include: (1) disparities in fluid velocities and properties, (2) non-hyperbolicity of the partial differential equations, (3) terms in nonconservative form, (4) surface tension force, etc. The manifestation of difficulties can be in stability, accuracy (e.g., unwanted oscillations and smearing), or uniqueness of solution. Moreover, even if the multifluid model is rendered hyperbolic, its eigensystem is still too complicated to be put in an analytical form, hence making it difficult to use the characteristic-based approximate Riemann solvers such as the Roe’s scheme or the Osher’s scheme. On the other hand, the simplicity of the flux vector splitting scheme [13] or the AUSM-family schemes [14] makes them an attractive alternative for the current multifluid model.

The AUSM⁺ scheme proposed by Liou [15] is known to be accurate and robust for compressible gas flows, especially for its ability in capturing shock and contact discontinuities. It can be easily extended to multispecies equations, and it can also handle flows of very low Mach number with the help of a pre-conditioning matrix and the numerical speed of sound [16]. While being successful for computing gas flows, the AUSM⁺ scheme is found to yield oscillatory solutions for liquid fluid, for which the equation of state, such as the stiffened gas model [17] or the Tait’s model, is stiff. To overcome this problem due to stiffness, new diffusion terms based on the pressure and velocity fields are introduced to the AUSM⁺ scheme [8], [18]. The modified scheme is essentially a variation of the AUSM⁺-up scheme of Liou [14], [19]. These diffusion terms are used to enhance the coupling between the pressure and velocity fields. We will show that they can effectively suppress numerical oscillations behind the pressure waves.

The rest of the paper is organized as follows. Section 2 shows the details of our method. We will discuss the stratified flow model, the discretization, the AUSM⁺-up scheme, exact Riemann solver for gas–liquid interface, the interfacial pressure correction, the updating procedure, and the extension to multidimensional flows. Section 3 presents results of several 1D and 2D problems. Finally, concluding remarks are given in Section 4.