A Characteristic Analysis of One-Dimensional Two-Fluid Model with Interfacial Area Transport Equation

Abstract

An interfacial area transport equation (IATE), proposed to dynamically describe the interfacial structure evolution of two-phase flows, could help improve the predictive capability of the two-fluid model. The present study aims to investigate the well-posedness issue of a one-dimensional two-fluid model with the IATE (named “two-fluid-IATE model” hereafter) using a characteristic analysis. The momentum flux parameters, which take into account the coupling of the volumetric fraction of phase and velocity distributions over the cross-section of a flow passage, are employed. A necessary condition for the system to achieve hyperbolicity under an adiabatic flow condition is identified. A case study is performed for an adiabatic liquid-liquid slug flow, which shows that the hyperbolicity of the two-fluid-IATE model is guaranteed if appropriate correlations of the momentum flux parameters are applied in the two-fluid-IATE model.

1. Introduction

In modeling of two-phase flows, it is realized that the one-dimensional two-fluid model is widely applied in many computer codes [1, 2]. Also, interfacial area transport equations (IATEs), which could capture the interfacial structure evolution through modeling the bubbles/droplets interactions (coalescence and/or disintegration) as well as phase change (boiling, evaporation, or condensation) [3–9], are developed in the framework of the two-fluid model. These IATEs belong to two categories: one-group and two-group models. The one-group IATE models small fluid particles, such as spherical bubbles, and therefore is only applicable to bubbly flows, while the two-group IATE takes into account both small and large fluid particles for applications in flows beyond bubbly flow. For simplicity, the two-fluid model with IATEs is named “two-fluid-IATE model” hereafter.

A successful model proposed for solving physical problems must first satisfy four principles, that is, equipresence, well-posedness, frame invariance, and fulfillment of the second law of thermodynamics. Mathematically, a well-posed problem means that the solution should exist, be unique, and depend continuously on variables. Insights into the field of the well-posedness of the two-fluid-IATE model are of great significance. In contrast to the single-phase Eulerian model, which has real characteristic roots inherently guaranteed by the second law of thermodynamics, the standard form of the partial differential equations in the two-fluid model may not be strictly hyperbolic as an initial value problem [10, 11]. Moreover, the infinite instability growth rate may exist when two-phase flows are modeled with the two-fluid model, which violates the observation that the energy originating from perturbation continuously decays due to the local fluctuations (energy cascade and dissipation), resulting in the naturally steady equilibrium of the flow in different flow regimes. These lead to the possibility of the partial differential equations of the two-fluid model being ill-posed.

In the literature, two fundamental technical approaches were utilized to remedy the ill-posedness of the two-fluid model [12, 13]. The first approach is more realistic, in which additional terms representing physical mechanisms, such as the pressure difference between the phases [14], virtual mass force [15], surface tension force [16, 17], and Reynolds stress [18], are incorporated into the conservation equations. However, difficulties were encountered in both constructing reliable constitutive relations for these terms and resolving the short wavelengths at which these physical mechanisms became dominant. The other approach is to regularize the model mathematically or numerically. For example, the ill-posedness of the two-fluid model can be remedied by either applying a preconditioning mass matrix [19] or adding an artificial numerical damping [12]. The disadvantage of the second approach lies in the possibility of introducing nonphysical dynamic behaviors. Furthermore, Song and Ishii [20] reported that the stable region of a well-posed one-dimensional two-fluid model could be evaluated by introducing momentum flux parameters, which helped preserve the information of the local flow structures, in which the void fraction profile is coupled with the velocity profile over the cross-section of a flow passage.

The current study aims at investigating the hyperbolicity of the one-dimensional two-fluid-IATE model by employing the momentum flux parameters. A characteristic analysis of the partial differential equations in the one-dimensional two-fluid-IATE model for an adiabatic flow was conducted to determine the necessary condition of the stability. The criterion was then developed in terms of the momentum flux parameters and local flow variables, such as the void fraction and phase velocities. Finally, a case study was performed for an adiabatic liquid-liquid two-component slug flow, in which the momentum flux parameters are correlated with the total void fraction of the dispersed phase based on experimental data [9]. The characteristic roots of this two-fluid-IATE model are examined graphically under different flow conditions to reveal the effects of the momentum flux parameters on regularizing the two-fluid-IATE model.

2. Governing Equations of the Two-Fluid-IATE Model

The current paper focuses on discussing the well-posedness of the standard one-dimensional two-fluid-IATE model, which can be obtained directly by area averaging the three-dimensional two-fluid-IATE model. To quantify the coupling effect of the local void fraction and velocity distributions, the momentum flux parameter 𝐶_𝑣,𝑘 is defined for field-𝑘 [20]:𝐶_𝑣,𝑘=𝛼_𝑘𝑣²_𝑘𝛼_𝑘𝑣_𝑘², (2.1) where 𝑣_𝑘 is the velocity of field-𝑘. The signs ⟨⟩ and ⟨⟨⟩⟩ are introduced to represent area averaging and phase fraction-weighted area averaging, respectively, which are defined for a parameter 𝜓 as ⟨𝜓⟩=(∫_𝐴𝜓d𝐴)/𝐴 and ⟨⟨𝜓⟩⟩=⟨𝛼_𝑘𝜓_𝑘⟩/⟨𝛼_𝑘⟩, respectively. Here, 𝐴 is the flow cross-sectional area, and 𝛼_𝑘 is the void fraction of field-𝑘. In the conventional one-dimensional two-fluid model, the momentum flux parameter is assumed to be unity, resulting in an ill-posed one-dimensional two-fluid model.

Dropping the signs ⟨⟩ and ⟨⟨⟩⟩ for simplicity, the area-averaged continuity and momentum equations of the two-fluid model and two-group IATEs for an adiabatic incompressible two-phase flow can be written as 𝜕𝛼₁𝜕𝑡+𝑣₁𝜕𝛼₁𝜕𝑧+𝛼₁𝜕𝑣₁Δ𝑚̇𝜕𝑧=−¹²𝜌₁,𝜕𝛼₂𝜕𝑡+𝑣₂𝜕𝛼₂𝜕𝑧+𝛼₂𝜕𝑣₂=Δ𝑚̇𝜕𝑧𝜌¹²₂,−𝜕𝛼₁𝜕𝑡−𝑣₃𝜕𝛼₁−𝜕𝑧𝜕𝛼₂𝜕𝑡−𝑣₃𝜕𝛼₂+𝜕𝑧1−𝛼₁−𝛼₂𝜕𝑣₃𝜌𝜕𝑧=0,₁𝑣₁𝜕𝛼₁𝜕𝑡+𝜌₁𝐶_𝑣,1𝑣²₁𝜕𝛼₁𝜕𝑧+𝛼₁𝜌₁𝜕𝑣₁𝜕𝑡+2𝛼₁𝜌₁𝐶_𝑣,1𝑣₁𝜕𝑣₁𝜕𝑧+𝛼₁𝜕𝑝𝜕𝑧=𝑀^∗_𝑖,1,𝜌₂𝑣₂𝜕𝛼₂𝜕𝑡+𝜌₂𝐶_𝑣,2𝑣²₂𝜕𝛼₂𝜕𝑧+𝛼₂𝜌₂𝜕𝑣₂𝜕𝑡+2𝛼₂𝜌₂𝐶_𝑣,2𝑣₂𝜕𝑣₂𝜕𝑧+𝛼₂𝜕𝑝𝜕𝑧=𝑀^∗_𝑖,2,−𝜌₃𝑣₃𝜕𝛼₁𝜕𝑡−𝜌₃𝐶_𝑣,3𝑣²₃𝜕𝛼₁𝜕𝑧−𝜌₃𝑣₃𝜕𝛼₂𝜕𝑡−𝜌₃𝐶_𝑣,3𝑣²₃𝜕𝛼₂+𝜕𝑧1−𝛼₁−𝛼₂𝜌₃𝜕𝑣₃𝜕𝑡+21−𝛼₁−𝛼₂𝜌₃𝐶_𝑣,3𝑣₃𝜕𝑣₃+𝜕𝑧1−𝛼₁−𝛼₂𝜕𝑝𝜕𝑧=𝑀^∗_𝑖,3,−23𝐷−𝜒^𝑐𝐷_𝑠𝑚,1²𝑎_𝑖,1𝛼₁𝜕𝛼₁−2𝜕𝑡3𝐷−𝜒^𝑐𝐷_𝑠𝑚,1²𝑎_𝑖,1𝛼₁𝑣₁𝜕𝛼₁+1𝜕𝑧3𝐷+𝜒^𝑐𝐷_𝑠𝑚,1²𝑎_𝑖,1𝜕𝑣₁+𝜕𝑧𝜕𝑎_𝑖,1𝜕𝑡+𝑣₁𝜕𝑎_𝑖,1𝜕𝑧=Φ₁,𝐷−𝜒^𝑐𝐷_𝑠𝑚,1²𝑎_𝑖,1𝛼₁𝜕𝛼₁𝐷𝜕𝑡−𝜒^𝑐𝐷_𝑠𝑚,1²𝑎_𝑖,1𝛼₁𝑣₁𝜕𝛼₁𝐷𝜕𝑧−𝜒^𝑐𝐷_𝑠𝑚,1²𝑎_𝑖,1𝜕𝑣₁−𝜕𝑧2𝑎_𝑖,23𝛼₂𝜕𝛼₂−𝜕𝑡2𝑎_𝑖,23𝛼₂𝑣₂𝜕𝛼₂+1𝜕𝑧3𝑎_𝑖,2𝜕𝑣₂+𝜕𝑧𝜕𝑎_𝑖,2𝜕𝑡+𝑣₂𝜕𝑎_𝑖,2𝜕𝑧=Φ₂. (2.2) Here, the subscripts 1, 2, and 3 denote field-1 (group-1 dispersed fluid particles), field-2 (group-2 dispersed fluid particles), and field-3 (continuous phase), respectively. Here, 𝑝,𝜌_𝑘,Δ𝑚̇₁₂,𝑀^∗_𝑖,𝑘,𝜒,𝐷_𝑐,𝐷_𝑠𝑚, and Φ_𝑘 are the pressure, density of field-𝑘, net intergroup mass transfer rate from group-1 to group-2 fluid particles, generalized drag force of field- 𝑘, coefficient accounting for the inter-group void transport at the group boundary, critical bubble diameter at the group boundary, fluid particle Sauter mean diameter, and interfacial area concentration source/sink rate for field-𝑘, respectively. In order to highlight the roles that the momentum flux parameters play in the process of regularization, we assume that the density of each phase does not change appreciably, both spatially and temporally.

In a gas-liquid two-phase flow, group-1 and group-2 fluid particles are small (spherical and distorted) and large (cap, slug, and churn-turbulent) bubbles, respectively [6, 7]. Therefore, it is reasonable to assume that 𝜌₁is identical to 𝜌₂, since field-1 and field-2 fluid particles are the same fluid with different particle sizes. In addition, there is no mass transfer considered between two different kinds of fluids due to the adiabatic assumption. In other words, no mass transfer takes place either between field-1 and field-3 or between field-2 and field-3. It is also noted that the phase-to-interface pressure difference is neglected, which means a single pressure is shared by all the phases.

In the literature, the contributions of 𝑀^∗_𝑖,𝑘 to the regularization of two-fluid model have attracted attention of many researchers [15, 16]. It was suggested that the inclusion of 𝑀^∗_𝑖,𝑘, including the virtual mass force and drag force, could help render the model into a hyperbolic form. However, as emphasized by Song and Ishii [20], to construct an appropriate correlation of 𝑀^∗_𝑖,𝑘 that could be applicable in any circumstances remained a challenge. Noting the argument that 𝑀^∗_𝑖,𝑘 usually tends to stabilize the two-fluid model, the analysis of the model hyperbolicity would be conservative if one adopts an algebraic form for 𝑀^∗_𝑖,𝑘 since the contribution of 𝑀^∗_𝑖,𝑘 to stabilizing the model will then not be credited in the characteristic analysis.

3. Characteristic Analysis

A characteristic analysis for the one-dimensional two-fluid-IATE model was performed to study the hyperbolicity of the model. The analysis provides information on the limiting short wavelength behavior. If a vector 𝐱=(𝛼₁,𝛼₂,𝑣₁,𝑣₂,𝑣₃,𝑝,𝑎_𝑖,1,𝑎_𝑖,2) is defined, then the governing equations of the system can be rewritten in the form of matrix as[𝐴]𝜕𝐱+[𝐵]𝜕𝑡𝜕𝐱=[𝐶]𝜕𝑧, (3.1) where [𝐴] and [𝐵] are the coefficient matrices. The detailed matrices of [𝐴], [𝐵], and [𝐶] are given in the appendix. The characteristic analysis of (3.1), which is a set of quasilinear partial differential equations, is analogous to the analysis of linear partial differential equations. The necessary condition for the model to be a well-posed initial value problem is that the set of partial differential equations is hyperbolic, that is, the characteristic equation |[𝐴]𝜆−[𝐵]|=𝐹(𝜆)=0 has real roots for 𝜆. The complex conjugate values for the eigenvalue 𝜆 lead to the elliptic equations, which are ill-posed mathematically. In this analysis for the system of (2.2), 𝐹(𝜆) is rearranged and simplified to 𝑔(𝜆)=𝑔₁(𝜆)−𝑔₂(𝜆), (3.2) where the corresponding polynomials are given as 𝑔₁(𝜆)=𝜌₃𝑔₁₁(𝜆)𝑔₁₂(𝜆),𝑔₁₁(𝜆)=𝜆²−2𝐶_𝑣,3𝑣₃𝜆+𝐶_𝑣,3𝑣²₃,𝑔₁₂(𝜆)=𝛼₁+𝛼₂𝜆²−2𝛼₁𝐶_𝑣,2𝑣₂+𝛼₂𝐶_𝑣,1𝑣₁𝜆+𝛼₁𝐶_𝑣,2𝑣²₂+𝛼₂𝐶_𝑣,1𝑣²₁,𝑔₂(𝜆)=−1−𝛼₁−𝛼₂𝜌₁𝑔₂₁(𝜆)𝑔₂₂(𝜆),𝑔₂₁(𝜆)=𝜆²−2𝐶_𝑣,1𝑣₁𝜆+𝐶_𝑣,1𝑣²₁,𝑔₂₂(𝜆)=𝜆²−2𝐶_𝑣,2𝑣₂𝜆+𝐶_𝑣,2𝑣²₂. (3.3) Note that 𝐹(𝜆) is replaced by 𝑔(𝜆) since four out of the eight roots for 𝐹(𝜆)=0 are guaranteed to be real roots.

Considering the fact that 𝜌₁ and 𝜌₃ are positive and the constrains of (1−𝛼₁−𝛼₂) being between zero and unity to maintain its physical meaning, we recognize that the coefficient of the highest order of 𝑔(𝜆), that is, the 4th order, should be always positive. All roots of 𝑔(𝜆) are real if the following two requirements are satisfied:(1)𝑔^′(𝜆)=0 has three real roots 𝜆₁,𝜆₂, and𝜆₃ (let 𝜆₁>𝜆₂>𝜆₃);(2)𝑔(𝜆₂)>0, 𝑔(𝜆₁)<0, and 𝑔(𝜆₃)<0.

It is realized that, to solve the eigenvalues, a 3rd-order polynomial equation is required to be solved, from which the expression of its roots are usually very complicated and lengthy. Therefore, this analytical approach is not taken here.

Alternatively, based on available experimental data of a certain two-phase flow configuration, the momentum flux parameters could be constructed as a function of the void fraction (and perhaps other flow variables), reducing the numbers of undetermined variables in 𝑔(𝜆). The plots of 𝑔(𝜆) could therefore be given, from which it becomes straightforward to study graphically the hyperbolicity property of 𝑔(𝜆). If there are four intersections between 𝑔(𝜆) and the real axis, then the characteristic equation 𝑔(𝜆)=0 will have four real roots for 𝜆. More details are given in the next section.

4. Results and Discussions: A Case Study

To demonstrate the feasibility of the studies discussed above, 𝑔(𝜆) is analyzed for a specific two-phase flow as a case study. Since the compressibility of fluids is not considered here for simplicity, it is appropriate to take a liquid-liquid slug flow as an example. One relevant experiment in the literature was performed by Vasavada [9]. In this experiment, two immiscible fluids, namely, water and Therminol 59, flowed upward adiabatically in a vertical pipe with an inner diameter of 25 mm at the room temperature. Therminol 59, an organic engineering fluid, was the dispersed phase with a density of 972 kg/m³ at 20°C. Based on the sizes of droplets, the organic droplets were classified into field-1 and field-2, with small droplets being field-1. Water was the continuous phase as field-3 with a density of 998 kg/m³ at 20°C. Local flow variables in 17 different flow conditions were measured at three measurement elevations along the flow directions with 𝐿/𝐷_ℎ=13, 67, and 121, respectively [9]. Table 1 lists the cross-sectional area-averaged data for each flow condition at 𝐿/𝐷_ℎ=13.

Table 1 

Area-averaged data for Runs 13–29 at 𝐿/𝐷_ℎ=13 [9]*.

The momentum flux parameter 𝐶_𝑣,𝑘 (𝑘=1,2, and 3) are calculated according to the experimental data [9] when the total void fraction of the organic drops, 𝛼(𝛼=𝛼₁+𝛼₂), is in the range between 0.1 and 0.6. The results are given in Figure 1, where 𝐶_𝑣,𝑘 are plotted against α. From the plots, the relations between 𝐶_𝑣,𝑘 and α could be determined. Figure 1(a) illustrates that 𝐶_𝑣,1 is a relatively weak function of α. Most of 𝐶_𝑣,1 data points fall into the range around 1.5 except two points having a magnitude between 2 and 2.5. After a further investigation of the measured data, it is found that the experimental data from the measurements in the flow conditions associated with these two data points had zero values for group-2 void fraction at several radial locations. This indicates that the flow might be in the bubbly-to-slug transition region instead of slug flows. Considering the difficulties in obtaining accurate measurements in flow transition regions, it might be reasonable to discard those two points and assume that 𝐶_𝑣,1≃1.5. In Figure 1(b), although it shows a scattered distribution when α is small, 𝐶_𝑣,2 follows the power law profile of α to a good approximation as α becomes larger, which corresponds to the slug flow conditions of interest. Other data-fitting methods of constructing correlation of 𝐶_𝑣,2, such as a linear function or higher order polynomials, are also investigated, which generally provide smaller value of 𝑅₂, which is a measure of how good the curve fit is with the data and is bounded between zero and unity. A higher value of 𝑅₂ always provides a better fit of the regression to the data set. Finally, as shown in Figure 1(c), 𝐶_𝑣.3 fits a 2nd-order polynomial of α quite well.

(a)

(b)

(c)

Figure 1 

𝐶𝑣,𝑘 versus the total void (drop) fraction α: (a) 𝐶𝑣,1, (b) 𝐶𝑣,2, and (c) 𝐶𝑣,3.

As a result of the above analysis, three assumptions are made as follows:(1)𝐶_𝑣,1(2)𝐶≃1.5;_𝑣,2=0.8824𝛼^−0.747;(3)𝐶_𝑣,3=−0.9028𝛼²+0.0552𝛼+1.1555. (4.1) Here, α stands for the total void fraction of Therminol 59 drops. It should be noted that for a more general analysis, the momentum flux parameters may be better correlated with additional parameters.

In addition, it is reasonable to assume that the velocities of the different phases are the same, that is, 𝑣₁=𝑣₂=𝑣₃=𝑣, because of the similar densities between water and Therminol 59. (the density difference between the two fluids is about 2.7% at 20°C.) Since the velocity for each phase follows similar turbulent velocity profile along the radial direction of the pipe, this assumption could also be partially verified from comparisons of the area-averaged velocities between group-1 and group-2 drops, as shown in Table 1. In most runs, 𝑣_𝑔1 is approximately equal to 𝑣_𝑔2, except for three runs, that is, Runs 16, 20, and 27. In these three runs, it is shown that the void fraction of group-2 droplets is significantly smaller than that of group-1 droplets. Also, the average velocity of the dispersed phase is found to be close to that of group-1 drops. These indicate that in Runs 16, 20, and 27, the numbers of group-2 drops measured are quite small, and the flows fall into bubbly-to-cap bubbly transition region. The zero relative motion assumption reduces the two-fluid-IATE model into a nonslip two-fluid model, similar to the homogeneous equilibrium model (HEM).

Furnishing the above assumptions, that is, (4.1) and equal phase velocities, the characteristic function 𝑔(𝜆) is simplified to a function of the void fraction of group-1 𝛼₁, void fraction of group-2 𝛼₂, characteristic root 𝜆, and velocity 𝑣, only. In addition, if one introduces a new parameter, 𝑤 as the dimensionless velocity, and defines it as 𝑤=𝜆/𝑣, then the characteristic equation 𝑔(𝜆)=0 converts to 𝑔(𝑤)=0. The roots of 𝑔(𝑤)=0 have the same mathematical properties as those of 𝑔(𝜆)=0, that is, if 𝑔(𝑤)=0 has four real roots, so does the original equation 𝑔(𝜆)=0.

Figures 2, 3, 4, 5, 6, 7, and 8 provide plots of 𝑔(𝑤) versus 𝑤 for 17 flow conditions. In these figures, subplot (a) is a plot in a relatively large range of 𝑤, while subplot (b) shows a plot in a sufficiently narrow range of 𝑤 for the benefit of observing the two roots near which 𝑔(𝑤) has small variations in the real axis. Figures 2–8 show clearly that 𝑔(𝑤) has four positive roots, which are consistent with the physical meanings. As illustrated in the plots, the typical shape of 𝑔(𝑤) is that 𝑔(𝑤) has one pair of roots that are close to each other and are approaching a less-than-unity 𝑤 value, whereas the other pair intersects at the values of 𝑤 in the order of unity. Attention should be paid to the first pair of roots, which demonstrates an essential role on the stability. For the flow conditions investigated, it is found out that the modified two-fluid-IATE model satisfies the necessary condition of the well-posedness, which indicates that the flows tend to be stable. Therefore, it is mathematically consistent with the HEM.

(a)

(b)

Figure 2 

𝑔(𝑤) for Runs 13, 17, and 18 over (a) a wide range of 𝜔 from 0 to 4 and (b) a narrow range of 𝜔.

(a)

(b)

Figure 3 

𝑔(𝑤) for Runs 14, 15, and 19 over (a) a wide range of 𝜔 from 0 to 3 and (b) a narrow range of 𝜔.

(a)

(b)

Figure 4 

𝑔(𝑤) for Runs 16, 21, and 29 over (a) a wide range of 𝜔 from 0 to 6 and (b) a narrow range of 𝜔.

(a)

(b)

Figure 5 

𝑔(𝑤) for Runs 20 and 27 over (a) a wide range of 𝜔 from 0 to 10 and (b) a narrow range of 𝜔.

(a)

(b)

Figure 6 

𝑔(𝑤) for Runs 22 and 26 over (a) a wide range of 𝜔 from 0 to 6 and (b) a narrow range of 𝜔.

(a)

(b)

Figure 7 

𝑔(𝑤) for Run 23 over (a) a wide range of 𝜔 from 0 to 3.5 and (b) a narrow range of 𝜔.

(a)

(b)

Figure 8 

𝑔(𝑤) for Runs 24, 25, and 28 over (a) a wide range of 𝜔 from 0 to 8 and (b) a narrow range of 𝜔.

The advantage of the two-fluid-IATE model with the momentum flux parameters over the conventional two-fluid model can be shown if the relative velocity between the phases is considered. Based on a force balance, the relative velocity between the dispersed phase and the continuous phase, 𝑈_𝑟, is obtained as [21]𝑈_𝑟𝑈_𝑟=43𝐷_𝑑Δ𝜌𝑔(1−𝛼)𝐶_𝐷𝜌₃, (4.2) where 𝐷_𝑑, 𝑔, 𝐶_𝐷, and 𝜌₃ are the bubble diameter, gravitational acceleration, drag force coefficient, and density of the continuous phase (field-3), respectively. It should be noted that the bubbles/drops are assumed spherical in (4.2). The drag force coefficient for distorted bubbles/drops is employed for field-1 as [22]𝐶_𝐷=2𝐷_𝑑3(1−𝛼)𝑔Δ𝜌𝜎, (4.3) where 𝜎 is the surface tension. Assuming that the diameter of the slug bubbles is 80% of the pipe inner diameter, the drag force coefficient for slug bubbles/drops is used for field-2 as [22]𝐶_𝐷=9.8(1−𝛼)³. (4.4) As a result, the relative velocity is a function of the total void fraction α as shown in Figure 9.

Figure 9 

Relative velocity versus the total void fraction α for the distorted and slug drops.

Our study shows that the eigenvalue 𝜆 becomes complex for the one-dimensional two-fluid three-field model with the unity momentum flux parameters if the relative velocities given by (4.2) are taken into account, indicating that the model is ill-posed (no hyperbolic region). In contrast, if the two-fluid-IATE model with the momentum flux parameters is considered, then the stability criteria, which ensures the hyperbolicity of the model, is shown in Figure 10 as a map with the continuous phase velocity 𝜈₃ as the 𝑦-axis and group-1 void fraction 𝛼₁ as the 𝑥-axis. When the total void fraction is smaller than 0.45, the model is always hyperbolic no matter what the continuous phase velocity is. When the total void fraction is between 0.45 and 0.65, the continuous phase velocity must increase to ensure the model hyperbolicity as the group-1 void fraction increases, indicating that the hyperbolic region becomes smaller. When the total void fraction is larger than 0.65, the model is unstable, which indicates the same mathematical property of the conventional two-fluid model.

Figure 10 

Stability criteria of the one-dimensional two-fluid-IATE model with momentum flux parameters.

Furthermore, the stability criteria can be represented using the slip ratios for field-1 and field-2 dispersed phase with respect to the continuous phase velocity, which are, respectively, defined as𝑆₁=𝑣₁𝑣₃,𝑆₂=𝑣₂𝑣₃. (4.5) The stability criteria for 𝑆₁ and 𝑆₂ are plotted in Figure 11. In the current case study for the upward liquid-liquid flow, the dispersed phase is lighter than the continuous phase; therefore, the slip ratios defined by (4.5) are always larger than unity. It is clearly seen that, except the case when the total void fraction is 0.45, the region of hyperbolicity decreases as the group-1 void fraction decreases. When the total void fraction reaches 0.65, the slip ratios slightly larger than 1, which is close to the HEM in order to ensure the model hyperbolicity. When the total void fraction is beyond 0.65, there is no hyperbolicity region if the momentum flux parameters given by (4.1) are adopted. Note that, to ensure model hyperbolicity, the criteria for 𝑆₁ and 𝑆₂ shown in Figure 11 should be satisfied simultaneously.

(a)

(b)

Figure 11 

Stability criteria of the one-dimensional two-fluid-IATE model with momentum flux parameters for (a) 𝑆1 and (b) 𝑆2.

5. Conclusions

This study examines the hyperbolicity property of the one-dimensional incompressible two-fluid-IATE model. In this study, momentum flux parameters are applied to take into account the effect of the void fraction and velocity profiles. These momentum flux parameters have physical meaning based on the local flow distributions and do not introduce any artificial terms. In addition, the proposed two-fluid-IATE model can be applied dynamically over a wide range of flow regimes, up to the churn-turbulent flow, including the flow transition regions. As a case study, a characteristic analysis has been performed to obtain the necessary conditions of the well-posedness of the one-dimensional incompressible two-fluid-IATE model for liquid-liquid slug flows in a circular pipe. Based on the available experimental data, the constitutive correlations of the momentum flux parameters were approximated as functions of the total void fraction of the dispersed phase. It was shown that the inclusion of the momentum flux parameters in the two-fluid-IATE model helps the model maintain hyperbolicity. The model’s hyperbolicity depends heavily on local flow parameters, such as the void fraction and phase velocities. When the relative velocity is considered, stability criteria pose restrictions on the continuous phase velocity to achieve model hyperbolicity.

This analysis only provides the necessary condition of the well-posedness of the two-fluid-IATE model. On the other hand, the question remains whether the instability growth rate for perturbation is unbounded as the amplitude of the perturbation decreases to zero. This requires a study of the dynamic behavior using the linear stability theory. The sufficient condition could be provided by using a wake equation with an accompanying infinitesimal perturbation [11, 12, 23] but has not been conducted for the current liquid-liquid applications. In addition, it remains questionable as how to generalize the results of incompressible model for the applications of gas-liquid two-phase flows, where the gas compressibility needs to be considered. The effect of the compressibility of the gas phase is essential for the validity of two-fluid-IATE applications in general two-phase flows [24]. Despite the aforementioned considerations, the introduction of the momentum flux parameters appears promising to the regularization of two-fluid-IATE

Appendix

 

The detailed matrices of [𝐴], [𝐵], and [𝐶] are given as follows:[𝐴]=𝜌1000000001000000−1−1000000₁𝑣₁0𝛼₁𝜌₁000000𝜌₂𝑣₂0𝛼₂𝜌₂0000−𝜌₃𝑣₃−𝜌₃𝑣₃001−𝛼₁−𝛼₂𝜌₃−20003𝐷−𝜒^𝑐𝐷_𝑠𝑚,1²𝑎_𝑖,1𝛼₁𝐷0000010−𝜒^𝑐𝐷_𝑠𝑚,1²𝑎_𝑖,1𝛼₁−23𝑎_𝑖,2𝛼₂[𝐵]=𝑣000001₁,0𝛼₁000000𝑣₂0𝛼₂0000−𝑣₃−𝑣₃𝜌00𝒵000₁𝑣₁𝐶_𝑣,1𝑣₁02𝛼₁𝜌₁𝐶_𝑣,1𝑣₁00𝛼₁000𝜌₂𝑣₂𝐶_𝑣,2𝑣₂02𝛼₂𝜌₂𝐶_𝑣,2𝑣₂0𝛼₂00−𝜌₃𝑣₃𝐶_𝑣,3𝑣₃−𝜌₃𝑣₃𝐶_𝑣,3𝑣₃00𝒰𝒵00𝒫0ℛ000𝑣₁02𝒬−3𝑎_𝑖,2𝛼₂𝑣₂1𝒮−3𝑎_𝑖,2000𝑣₂,[𝐶]=−Δ𝑚̇₁₂𝜌₁Δ𝑚̇₁₂𝜌₂0𝑀^∗_𝑖,1𝑀^∗_𝑖,2𝑀^∗_𝑖,3Φ₁Φ₂, (A.1)(1) where 𝒫 denotes −[2/3−𝜒(𝐷_𝑐/𝐷_𝑠𝑚,1)²](𝑎_𝑖,1/𝛼₁)𝑣₁, 𝒬 denotes −𝜒(𝐷_𝑐/𝐷_𝑠𝑚,1)²(𝑎_𝑖,1/𝛼₁)𝑣₁, ℛ denotes [1/3+𝜒(𝐷_𝑐/𝐷_𝑠𝑚,1)²], 𝒮 denotes −𝜒(𝐷_𝑐/𝐷_𝑠𝑚,1)²𝑎_𝑖,1, 𝒰 denotes 2(1−𝛼₁−𝛼₂)𝜌₃𝐶_𝑣,3𝑣₃, and 𝒵 denotes (1−𝛼₁−𝛼₂).

Acknowledgment

This work was performed under the auspices of US Nuclear Regulatory Commission (NRC) through a Faculty Development Grant (Award no.: NRC-38-08-945). The authors would like to thank the NRC for the support.