Abstract
In this article, two-phase compressible six equation flow model is numerically investigated. The six-equation model consists of velocity, pressure equations, and also relaxation terms. An extra seventh equation is included describing the total energy of the mixture to ensure the correct treatment of the sharp discontinuities. The model is hyperbolic and poses numerous difficulties for numerical schemes. An efficient and well-balanced scheme can handle the numerical difficulties related to this model. The second order space-time CE/SE scheme is extended to solve the model. This scheme offers an effective numerical method for several continuum mechanics problems. The suggested scheme suppresses the numerical oscillations and dissipation effects. Several numerical test cases have been carried out to reveal the efficiency and performance of the proposed approach. The results are compared with the exact solution and also with Runge-Kutta Discontinuous Galerkin (RKDG) and central (NT) schemes.
1. Introduction
A lot of research work has been drafted on multiphase flows caused by their various applications in divergent scientific and engineering fields [1].
In many cases, multiphase flow problems arise in technical and advanced heat transfer systems. They have industrial implementation in energy conversion sector, paper creation, food making, and medical applications. Furthermore, an effective and secure evaluation of nuclear reactors accident situation needs an accurate forecasting of two-phase flow structure. Moreover, two-fluid flow has some daily operations in geysers and boilers. Ishii classified two-phase flow into four different groups on the grounds of phase mixture in the flow [2]. The groups contain gas-liquid flow (i.e., gas droplet flows, bubble flows, and separated flows), gas-solid flow (i.e., fluidized flows and flow of gas particle), liquid-solid flow (i.e., flow of slurry and transport of sediment), and immiscible liquid-liquid flow. As two-phase flows have substantial implementations, appropriate mathematical templates are mandatory to be constructed and associated with hypothetical outcomes. The averaging methodology is usually the foundation of multiphase flows as there is intricate connection between the fluids in a multiphase flow approach [3]. Researchers and engineers have proposed a number of multiphase mathematical flow models in the past to study the physical behavior of such flows. One such type of those models is seven-equation models; they use different velocity, density, and pressure for two phases. The seventh equation is obtained by coupling of conservation laws and the convection equation for the interface motion. Baer and Nunziato [4] were the pioneers who formulated a solid gas two-phase seven-equation model; and later, this model was modified and further studied by Abgrall and Saurel [5, 6] and other investigators. Although, such models are treated as the best two-phase flow models. However, there exist a number of physical and numerical difficulties in such models. To overcome those problems present in such flow models, the investigators have suggested models with a smaller number of equations ranging from three to six [7, 8].
In this work, we have primarily considered the six-equation two-phase flow model studied by [9]. As mentioned earlier that because of complex nature of the seven-equation model, Kapila et al. [7] extracted a reduced five-equation model from seven-equation model [4]. The extracted model [7] has been used successfully to study the two-phase flows. Despite the fact that this model is simple but the presence of nonconservative term and shocks make it difficult to obtain the convergent physical solution.
To resolve difficulties, present in the Kapila’s model, Saurel and coauthors [9] proposed a new model by incorporating the nonequilibrium pressure terms present in the seven-equation model. This new model is the six-equation model with same velocity, separate pressure, and relaxation terms for each phase. To treat the shock correctly in the single phase limit, an additional equation of mixture energy has been incorporated in the model. Importantly, this new model is hyperbolic in nature with three wave propagation speeds. Moreover, the positivity of the volume fraction is an important feature present in the model. Saurel and coauthors [9] applied Godunov-type method and HLLC-type solver to solve the six-equation model. Zia et al. [10] applied KFVS and central-upwind schemes to study the model under discussion.
In this article, a space-time CE/SE scheme [11] is extended to solve the reduced six-equation model equations. The space-time CE/SE scheme is different from the existing methods in literature. In CE/SE scheme, the space and time variables are discretized simultaneously, whereas in other methods, the space variable is discretized first, and the resulting system of ODE’s is solved by any other method. Because of the unique idea of discretization, the CE/SE scheme effectively handles the numerical complexities of the two-layer equations. This scheme has been successfully implemented in different areas, e.g., see unsteady flows [12, 13]. The solution of two-phase six-equation model is also available in the literature by upwind scheme and kinetic flux vector splitting (KFVS) scheme [10]. The results obtained by CE/SE scheme are compared with RKDG scheme [14] and central (NT) scheme [15].
This article is outlined as follows. The 1-D two-phase compressible six-equation flow model is presented in Section 2. A review of RKDG scheme is given in Section 3. Section 4 is devoted to numerical test problems. Conclusions are explained in Section 5.
2. Mathematical Model
Here, we present reduced six-equation compressible flow model which has been deduced from the seven-equation model by using asymptotic limit condition of zero velocity relaxation time in [7]. It is observed that the considered model has the ability to resolve those difficulties which were involved in five-equation model presented by [7, 16]. It is important to mention that we are considering the model without heat and mass transfer mechanisms. The one-dimensional model is given as [10].
Here, , , and are the velocity, the mixture density, and the mixture pressure, respectively; and volume fraction for phase is given by . The mixture energy and mixture internal energy are defined in the following equations
For each phase, mass fraction is . Furthermore, the density of the mixture of both phases is defined by . The equation-of-state (EOS) formula is given by
Here, , , and are constants (material) subject to the fluid. Utilizing equation (8), the mixture pressure can be found as under
Here, the interfacial pressure is denoted by and is given by the formula
The acoustic impedance of phase is given by the term . For the seventh equation in the model, the mass and momentum equations are combined with the internal energy equation, and this additional equation becomes where . To handle the nonconservative terms on the RHS of the model and for convergence of the schemes, nonconventional jump conditions are needed which are based on Reimann solver, see [9]. It is worth specifying that the numerical schemes proposed to solve the current model do not need these conditions. Moreover, the frozen sound speed is given as under where is sound speed of phase and is given as
Now, equation (1) can be reformulated as
On utilizing equations (1) to (6) along with equation (15), the model takes the following form, where
The expressions and for involve derivatives; therefore, these expressions are treated in the similar fashion as the flux vectors. For closure of the model equations (1) to (6) and (12), extra equations are required.
Thus, in the present situation, these equations are obtained from the stiffened gas EOS (9). In the next section, we will present a brief derivation of RKDG and CE/SE schemes.
3. Proposed Schemes
3.1. RKDG Scheme
A short overview of RKDG scheme is presented in this section. The RKDG scheme is applied on system of hyperbolic conservation laws of the form where with initial conditions
We subdivide the domain into cells , where We denote the center of cell by and defined as 0. . The size of cell is defined as , for all . Multiply equation (18) by a smooth function and integrate the resulting equation over the cell . By integrating, we obtain
Consider for each time , the solution approximation of belongs to discontinuous Galerkin finite element space
The is the polynomial space of degree at most in cell . We choose the local orthogonal basis , named as scaled Legendre polynomials , over cell as defined below
Then, the approximate solution is defined as where by using the orthogonal property
Now, replace a smooth function by the test function and flux function by numerical flux in equations (20) and (21) to obtain a numerical scheme as follows with
The flux is any monotone numerical flux. Here, we use the LFF as a monotone numerical flux which is defined below where and are the limited values of approximate solution at the cell interface . By using above definitions, equations (26) and (27) simplify to
Next, to obtain nonoscillatory discontinuity propagation and high order accuracy by RKDG method, the limiter procedure is divided into two parts. First, identify the cells, named as “troubled cells,” that may need the limiting procedure. Second, reconstruct the polynomial solutions in these troubled cells by using WENO reconstruction.
Here, the troubled cells are identified by total variation bounded (TVB) limiter, which is described as where
Now, and are modified by the modified minmod (MM) function, as follows where MM is defined as where is a positive constant and the minmod function is defined as
Using explicit, nonlinearly stable high order Runge-Kutta time discretization. [Shu and Osher, JCP, 1988].
The semidiscrete scheme (26) is written as is discretized in time by a nonlinearly stable Runge-Kutta time discretization, e.g., the third order version.
3.2. The CE/SE Scheme
In this section, we give a brief derivation of scheme [1, 11].
Let the coordinates of a 2-D Euclidean space () be and . The equivalent integral form can be obtained by employing the Gauss-divergence theorem, on equation (15). in which (a) , i.e., for all and be the elements of the vector in the and directions, respectively, and (b) def in which represents the area whereas being the typical unit outward vector of a surface component on . Over the space-time region CE (conservation element), equation (38) is implemented. The actual numerical combination uses solution components (SEs) in a discreet manner.
In space, the arrangement of work focuses is represented by where and for each . There is a SE (solution element) corresponding to each , as depicted in Figure 1 (dashed curve). It includes both a horizontal and a vertical line and the immediate neighborhood.
(a)
(b)
(c)
The correct size of immediate vicinity does not make a difference. For any and are estimated by , , and separately as takes the following form
Moreover,
Using chain rule, we acquire where
Here, , where and , respectively, represent the column and row indices. Note that , and are constants in SE . The numerical analogs of these values are the following and at , respectiyely. Since , therefore,
Using equations (39), (40), (41), and (45), we get
The summation gives where
Clearly from the source term, . Moreover,
For more details, the reader is referred to [1, 11].
4. Numerical Test Problems
In the present section, the scheme is used to study the single velocity six-equation model along with a pressure relaxation procedure to accommodate the nonconservative terms in the model. It has been mentioned earlier that two-phase flow models inherit serious computational difficulties. Some one-dimensional test problems carried out in this section and their results are compared with each other, to authenticate the results of two proposed schemes. To further authenticate the efficiency of suggested schemes, we have used the RKDG [14] and central (NT) [15] schemes.
Problem 1 No reflection test problems. This test case was also studied in [8] by implementing KFVS scheme. The Riemann left and right data about interface is given as
The problem is simulated at 500 mesh cells, and is the final simulation time. Here, and are values of specific heat ratios. Furthermore, , and . Because of significant pressure shifts at the interface, it is one of the tough test problems for a numerical scheme. Choosing the velocity and jumping pressure over the shock makes it difficult to create reflection wave. Figure 2 displays the results obtained from two schemes. Clearly, wiggles can be observed in pressure and velocity plots from all the schemes. These wiggles disappear on refined mesh. Similar type of results is reported in Kreeft [17]. Furthermore, we have also computed the error of and RKDG schemes. The errors are given in Table 1. We can clearly observe that has less error as compared to RKDG scheme.
The -errors for CE/SE, KFVS, and central (NT) schemes at different grid points are listed in Table 1. We can clearly observe from Table 1 that CE/SE scheme has less error as compared to KFVS and central (NT) schemes.
Problem 2 Sode’s test problem. The problem considered is like the Sod problem of single-phase gas dynamics, which is simulated by [10] as well. A thin membrane fixed at separates the two gases which are at rest initially. The gases on both sides of the interface mixes with each other as soon as membrane is removed. The initial state of problem is as under
The specific heat ratio is on the of interface, and is on the right of interface. The results at are presented in Figure 3. From the simulated results, a left-going rarefaction wave can be seen. Also, a right-moving two-fluid interface and right-going shock wave can be observed. In Figure 3, the results by proposed schemes are shown at 400 cells. For this comparison, the reference solution is simulated for scheme at 2000 meshes. Both the schemes are in close agreement capturing the shocks effectively.
Problem 3 Two fluid mixture problem. This two-fluid mixture problem has been considered in [10] and was simulated by five-equation model. The values of the other physical variables are given as
(a)
(b)
(c)
(d)
In this case study, and . In the figure, one can see a right moving shock wave, a left going rarefaction wave, and a contact discontinuity. The right going shock hits the interface at . The shock continues to move towards right, and a rarefaction wave is created which is moving towards left. For simulation of the results, we take 400 cells, and is the final time. The results are presented in Figure 4. The results reveal that two schemes give comparable results. However, the results of scheme are better than the other schemes.
Problem 4 Comparison test problem. The present test of water-air mixture is presented here with and RKDG schemes, respectively. The initial values of the problem are [10] where , and . The results at are depicted in Figure 5. We can observe that all the schemes are in close agreement with each other. However, among all performance is better.
5. Conclusions
In this article, reduced six-equation flow model was numerically investigated by the space-time and RKDG schemes. The numerical complexities related to the model include accurate discretization of source term. The proposed methods have the ability to resolve sharp discontinuous profiles. The numerical test problems show the performance of the suggested method to handle a wide range of flow conditions. The results of suggested schemes are compared with those obtained from the central (NT) scheme and exact solutions. A good agreement was observed among the results of and RKDG schemes. However, it was found that our suggested scheme effectively resolves sharp discontinuities better than central (NT). The scheme generally proposes a good strategy for solving both nonconservative and conservative system of equations.
Data Availability
No data is required to perform this research.
Conflicts of Interest
The authors declare that they have no conflicts of interest.