Central Upwind Scheme for a Compressible Two-Phase Flow Model

In this article, a compressible two-phase reduced five-equation flow model is numerically investigated. The model is non-conservative and the governing equations consist of two equations describing the conservation of mass, one for overall momentum and one for total energy. The fifth equation is the energy equation for one of the two phases and it includes source term on the right-hand side which represents the energy exchange between two fluids in the form of mechanical and thermodynamical work. For the numerical approximation of the model a high resolution central upwind scheme is implemented. This is a non-oscillatory upwind biased finite volume scheme which does not require a Riemann solver at each time step. Few numerical case studies of two-phase flows are presented. For validation and comparison, the same model is also solved by using kinetic flux-vector splitting (KFVS) and staggered central schemes. It was found that central upwind scheme produces comparable results to the KFVS scheme.


Introduction
Multiphase flows are commonly observed in nature and science, from stand storms, to volcanic clouds, blood flow in vessels and motion of rain droplets. There are also numerous examples where multiphase flow occurs in industrial applications, for example energy conversion, paper manufacturing, food manufacturing as well as in chemical and process engineering. Due to their wide range applications, suitable models are required for the accurate prediction of the physical behavior of such flows. However, modeling and simulation of flows is a complex and challenging research area of the computational fluid dynamics (CFD).
Multiphase flow problems involve the flow of two or more fluids separated by sharp interfaces. The coupling of interfaces with the flow model is a challenging part in the simulation of such flows, as coupling miss-match may introduce large errors in the numerical simulations. It is important to mention that this work is only concerned with two-phase flow problem.
Several two-phase flow models exist in the literature for describing the behavior of physical mixtures. These models use separate pressures, velocities and densities for each fluid. Moreover, a convection equation for the interface motion is coupled with the conservation laws of flow models. In the literature such models are known as seven-equation models. One of such models for solid-gas two-phase flows was initially introduced by Baer and Nunziato [1] and was further investigated by Abgrall and Saurel [2,3], among others. The seven-equation model is considered as the best and established two-phase flow model. However, it inherit a number of numerical complexities. To resolve such difficulties researchers have proposed reduced models containing three to six equations [4][5][6].
The Kapila's five-equation model [4] deduced from Baer and Nunziato seven-equation model [1] is a well known reduced model and has been successfully implemented to study interfacing compressible fluids, barotropic and non-barotropic cavitating flows. The model contains four conservative equations, two for mass conservation, one for total momentum and one for total energy conservation. The fifth equation is a non-zero convection equation for the volume fraction of one of two phases.
Although, this five-equation model is simple, but, it involves a number of serious difficulties. For example, the model is non-conservative and hence it is difficult to obtain a numerical solution which converges to the physical solution. In the presence of shocks, the volume fraction may become negative. Another issue is related to non-conservative behavior of the mixture sound speed [7].
To overcome the associated difficulties of Kapila's five-equation model, Kreeft and Koren [5] introduced a new formulation of the Kapila's model. The new model [5] is also non-conservative containing two equations for the conservation of mass, one each for mixture momentum and total energy respectively. The fifth equation is the energy equation for one of the two phases and it includes source term on the right hand side which represents the energy exchange between two fluids in the form of mechanical and thermodynamical work. In the current model, the first four equations are conservative and the non-conservativity in the models is due to the energy exchange term in the fifth equation. Consequently, the implementation of finite volume type schemes are relatively convenient to such models.
Very recently, diffuse interface method, finite volume WENO scheme and discontinuous Galerkin method have been used to solve the multiphase flow models [8][9][10]. However, in this work, the central upwind scheme [11] is proposed to solve the same five-equation model [5]. The proposed scheme uses information of local propagation speeds and estimates the solution in terms of cell averages. Further, the scheme has an upwind nature, because it takes care of flow directions by means of one-sided local speeds. Moreover, this scheme can be extended to incompressible flow problems e.g. it can be extended to solve incompressible two-phase shallow flow model [12]. The suggested scheme is applied to both one and two-dimensional flow models. For validation, the results of central upwind scheme are compared with those obtained from the KFVS [13][14][15][16] and the non-oscillatory staggered central scheme [17,18]. The numerical results of the schemes are analyzed qualitatively and quantitatively. It was found that the proposed central upwind scheme has comparable results to the KFVS scheme and are more accurate than the staggered central scheme.

One-dimensional two-fluid flow model
In this section, the one-dimensional two-fluid flow model [5] is presented. The selected model is the reformulation of original five-equation model of Kapila et al. [4]. Here, it is assumed that both fluids are mass conservative and have same velocity and pressure on both sides the interface. Moreover, heat conduction and viscosity are not considered. In this model, first four equations describe the conservative quantities: two for mass, one for over all momentum and one for total energy. The fifth equation is the energy equation and it includes source term on the right hand side which gives the energy exchange between two fluids in the form of mechanical and thermodynamical work. The state vector q of primitive variables has the form q = (ρ, u, p, α) T . Here, the bulk mixture density is denoted by ρ, u = (u, 0, 0) are the bulk velocities along each characteristic direction, p denotes the bulk pressure and α represents the volume fraction of fluid 1. This means that a part α of a small volume dV is filled with fluid 1 and a part (1-α) with fluid 2.
For bulk quantities, such as mixture density ρ and mixture total energy E, we assume that α is a volume fraction of fluid 1 and (1-α) of fluid 2. Using these conventions, we can define and the total energies of each fluid as where e 1 and e 2 denote the internal energies of fluid 1 and fluid 2, respectively. The internal energies e 1 and e 2 are given in terms of their respective densities and pressure through equations of state In one space dimensions, the two-fluid flow model can be written as [5] w t þ f ðwÞ x ¼ s ; where w ¼ ðr; ru; rE; r 1 a; r 1 E 1 aÞ T ; ð4bÞ sðwÞ ¼ ð0; 0; 0; 0; 0; s 5 Þ T : Here, w represents the vector of conservative variables, f is the vectors of fluxes, s is a vector of source terms with only last non-zero term. The term s 5 represents the total rate of energy exchange per unit volume between fluid 1 and fluid 2 and is equal to the sum of rates of mechanical s 5 M and thermodynamic s 5 T works [5], i.e. s 5 ¼ s 5 M þ s 5 T , with The term b ¼ r 1 a r represents the mass fraction of fluid 1, while the relations t 1 ¼ 1 denote the isentropic compressibilities of both fluids. Here, c 1 and c 2 represent the sound speeds of fluid 1 and fluid 2. The bulk isentropic compressibility is defined as Assume that the equations of state in Eq (3) are the stiffened equations of state [19] r where γ i , π i are the material specific quantities. Therefore, the sound speeds in each fluid are given as The expressions for the sound speeds are normally obtained from the second law of thermodynamics. The total energies of fluids 1 and 2 can be given as Using Eqs (4), (10) and (11), we obtain the primitive variables as where Here, w i ; i = 1,. . .,5, represent the components of w, the vector of conserved variables. Moreover, in Eqs (12)-(14) the primitive variables are expressed in terms of conserved variables. While in Eq (13) the positive sign is chosen for (π 2 γ 2 −π 1 γ 1 ) > 0 and negative otherwise. Because of Eq (9) One dimensional Central upwind scheme In this section, the central upwind scheme of Kurganov and Tadmor [11] is derived for the one-dimensional five-equation two-fluid flow model Eq (1). Let N represents the total number of discretization points and ðx iÀ 1 2 Þ i 2 f1; Á Á Á ; N þ 1g denotes the divisions of the given domain [0, x max ]. A uniform width Δx for each cell is considered, while, x i represents the cell centers and x iþ 1 2 refer to the cell boundaries.
Further, we take Moreover, The cell average values of conservative variables w i are defined as Integration of Eq (4) over the interval O i gives where H iþ 1 2 ðtÞ is the numerical flux defined by The first four components of the source vector s i in Eq (21) are zero and the fifth non-zero component is given as The numerical derivatives w x i are approximated through a nonlinear limiter which guarantees the positivity of the reconstruction procedure Eq (24).
Here, MM denotes the min-mod non-linear limiter MMfx 1 ; x 2 ; :::g ¼ Moreover, a iþ 1 2 ðtÞ represents the maximal local which in the generic case could be To achieve the second order accuracy in time, a second order TVD RK-method is applied to the Eq (21). For simplicity if the right hand side of the Eq (21) is taken as L(w), then two stage TVD RK-method to update w is given as under where w n is a solution at previous time step and w n+1 is updated solution at next time step. Moreover, Δt represents the time step.

Two-dimensional two-fluid flow model
In two-dimensional space, the two-fluid flow model can be written as [5] w t þ f ðwÞ x þ gðwÞ y ¼ s ; where w ¼ ðr; ru; rv; rE; r 1 a; r 1 E 1 aÞ T ; ð29bÞ f ðwÞ ¼ ðru; ru 2 þ p; ruv; ruE þ pu; r 1 ua; r 1 E 1 ua þ puaÞ; ð29cÞ Here, w represents the vector of conservative variables, f, g are vectors of fluxes in x and y directions, s is a vector of source terms with only last non-zero term. The term s 6 represents the total rate of energy exchange per unit volume between fluid 1 and fluid 2 and is equal to the sum of rates of mechanical s 6 M and thermodynamic s 6 T work [5], i.e. s 6 ¼ s 6 M þ s 6 T , with , the total energies of fluids 1 and 2 are given as Since the energy equation is directional independent, therefore for one-and two-dimensional problems the procedure of calculating primitive variables are the same. In two-dimensional space, the primitive variables can be retrieved in the following manner. Using Eqs (29), (32) and (33), we obtain > > > > > < > > > > > : where Similarly, as in case of one-dimensional model, w i ; i = 1,. . .,6, represent the components of w, the vector of conserved variables. In Eq (35) the positive sign is chosen for (P 2 γ 2 −P 1 γ 1 ) > 0 and negative otherwise.

Two dimensional Central upwind scheme
In this section, the central upwind scheme [11] is extended to two-dimensional five-equation two-phase flow model Eq (1). To implement the scheme, first we need to discretize the computational domain.
Let N x and N y be the large integers in x and y-directions, respectively. We consider a cartesian grid with rectangular domain [x 0 , x max ] × [y 0 , y max ] and it is covered by the cells C ij Here, the representative coordinates in a cell C ij are denoted by (x i , y j ). Further, we take and The cell average values conservative variable w i, j at any time t are given as The following linear piecewise interpolant is constructed as under wðx; y; tÞ ¼ X Here, χ i, j is the characteristics function corresponding to the cell x iÀ 1 2 ; x iþ 1 2 Â y iÀ 1 2 ; y iþ 1 2 , ðwÞ x i;j and ðwÞ y i;j are the approximations of x and y-derivatives of w at cell centers (x i , y j ). In two-dimensional case to compute the derivative a generalized MM limiter is used The intermediate values in the present case are given as For details and complete derivation of the scheme, see [11].

Numerical test problems
This section presents some numerical test problems (both one and two-dimensional) to check the capability of central upwind and KFVS schemes to compressible two-phase reduced fiveequation flow model. The results are compared with those obtained from central scheme [18].

One-dimensional test problems
In this section six one-dimensional test problems are considered to verify the efficiency and accuracy of the proposed schemes.
Sod's problem. The Sod's problem [6] is the well known test problem in the single phase gas dynamics. In this problem, gases are separated by a very thin membrane placed at x = 0.5 and are initially at rest. The left side gas has high density and pressure as compared to right side gas. After removing the membrane, the gases evolution in time take place. The initial data for the problem are given as Here, γ L = 1.4 and γ R = 1.2, P L = 0 = P R , and CFL = 0.5. This problem is also considered in [14]. It is a hard test problem for a numerical scheme. From the solution we can see a left moving rarefaction wave, a contact discontinuity, and a right moving shock wave. The right moving shock hits the interface at x = 0.5. The shock continues to move towards right and a rarefaction wave is created which is moving towards left. The results are simulated on 400 mesh cells and the final simulation time is taken as t = 0.012. The solutions are presented in the Fig 3. All the schemes give comparable results. However, from the zoomed graph it can be noted that KFVS scheme gives better resolution of peaks and discontinuities.
No-reflection problem. The initial data are given as The ratio of specific heats are γ L = 1.667 and γ R = 1.2. Moreover, P L = 0 = P R and CFL = 0.4. We discretize the computational domain [0, 1] into 500 mesh cells and the final simulation time is t = 0.02. This is a hard test problem for a numerical scheme due to large jumps in pressure at the interface. The choice of pressure and velocity jump over the shock prevents the creation of a reflection wave. Therefore, a shock wave moves to the right. The results are depicted in Fig 4. Wiggles can be seen in the velocity and pressure plots of all schemes, representing small waves that are reflected to the left. However, unlike real velocity and pressure oscillations, these wiggles reduces on refined meshes. Similar type of wiggles are also reported in the results of [5].
Water-air mixture problem. This one-dimensional problem corresponds to the water-air mixture [5,20]. The initial data are given as  5. Although the initial composition of the mixture is constant, it evolves in space and time. It can be observed that the three schemes give comparable results. Moreover, our results are in good agreement with the results in [20]. Water-air mixture problem. Again a one-dimensional water-air mixture problem [5,20] is considered. However, this problem differs from the previous problem by allowing changes in the mixture composition. The initial data are given as The ratio of specific heats are given as γ L = 1.4 and γ R = 1.6. We have chosen 200 mesh cells and the final simulation time is t = 0.1. This problem is a contact discontinuity of water-air density ratio. The numerical results are shown in Fig 7. The same problem was also considered in [5,6]. In this problem, both pressure and velocity are the same. Therefore, the interface is moving to the right with uniform speed and pressure. The numerical results show that KFVS

Two-dimensional test problems
To check the performance of proposed numerical scheme in two-dimensional space, we considered two test problems. In these problems the impact of a shock in air on a bubble of a lighter and a heavier gas is studied. Initially, Haas and Sturtevant [21] investigated these problems. Later, Quirk and Karni [22], Marquina and Mulet [23], Kreeft and Koren [5] and Wackers and Koren [6] also discussed these test cases. A schematic computational setup for these two problems is sketched in Fig 8. A shock tube of length 4.5 and width 0.89 is considered. The top and bottom walls of the tube are solid reflecting walls, while both ends of the tube are open. A cylinder of very thin cellular walls filled with gas is placed inside the tube. A shock wave is generated in the right end of the shock tube and is moving from right to left. After hitting by shock, the walls of the cylinder ruptures and the shock interacts with the gas of the cylinder. Due to fast interaction both gases do not mix in large amount, leading to a two-fluid flow problem. As the shock approaches to the surface of the bubble a reflected shock is generated from the surface of the bubble which moves towards right back in the air. At later time, this interaction become more and more complicated. The shock continues to move towards right in the air after passing through bubble and produces secondary reflected waves in the bubble when it hits the surface of the bubble. The wave patterns generated by interaction are strongly depending on the density of the gas inside the bubble. However, some of the waves can be observed in almost all cases [5,6]. Here, a light helium gas and a heavy R22 gas are considered inside the cylindrical bubble.
Helium bubble. In this problem, we study the interaction of Ms = 1.22 planar shock, moving in air, with a cylindrical helium bubble contaminated with 28% of air. The bubble is   The position of key features occurred during the time evolution are well explained in [5,6,23]. Therefore, we omit discussion on these features. The computational domain is discretized into 800 × 200 mesh cells. The contours for density, pressure and volume fraction are depicted in Figs 9, 10 and 11 at time: 0.25, 0.30, 0.35, 0.40. These results agree closely with the plots given in [5,6,21,22] at times: 32 μs, 52 μs, 62 μs, 82 μs. In Figs 10 and 11 the contours of pressure and volume fraction show a perfect splitting of the pressure waves and the interface. The shocks and interface are sharp during the simulation. As observed in [6], the last interface is slowly bending inwards in Fig 11. The phenomena will continue at later times until the bubble split in two vortices. The comparison between KFVS and central schemes can be clearly observed in Fig 12. R22 bubble. Here, the same Ms = 1.22 planar shock moving in air hits a cylindrical R22 bubble which has higher density and lower ratio of specific heats than air. This results in about two times lower speed of sound. For more details, the reader is referred to [5,6]. The initial data are given as The computational domain is discretized into 800 × 200 mesh cells. Due to the lower speed of sound, the shock in the bubble and the refracted shock lag behind the incoming shock. Moreover, due to the circular shape of the bubble the refracted, reflected and shock waves are curved. The results for density, pressure and volume fraction are displayed in Figs 13, 14 and 15 at times: 0.35, 0.60, 0.70, 0.84, 1.085, 1.26. These results shows good agreement with the results [5,6,21,22] at times: 55 μs, 115 μs, 135 μs, 187 μs, 247 μs, 318 μs. The splitting in pressure and interface is observed in flow pattern of density contours. Moreover, no wiggles are visible in our results and the pressure is continuous over the interface. Hence, the numerical results of our scheme reflect all key features as explained in [5,6,21].

Conclusions
A central upwind finite volume scheme was extended to solve the compressible two-phase reduced five-equation model in one and two-dimensional space. The suggested scheme is based on the estimation of cell averages by using the information of local propagation speed. In twodimensional space the scheme is implemented in a usual dimensionally split manner. The nondifferential part of the source terms are approximated by the cell averaged values, whereas the differential part terms are approximated similar to the convective fluxes. To preserve the positivity of the scheme a min-mod non-linear limiter is used. To achieve the second order accuracy in time a TVD Runge-Kutta method is utilized. For validation, the results of the proposed numerical scheme are compared qualitatively and quantitatively with those of KFVS and staggered central schemes. Good agreements were observed in the results of all three schemes. It was found that in some test problems central upwind scheme produced more accurate results, while KFVS scheme performed well in other problems. Perhaps, this is due to the reason that both the schemes are upwind biased. The staggered central scheme was found diffusive in all test problems.