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

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 twofluid 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.


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 twofluid 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 twofluid model with IATEs is named "two-fluid-IATE model" hereafter.fraction and velocity distributions, the momentum flux parameter C v,k is defined for field-k 20 : where v k is the velocity of field-k.The signs and are introduced to represent area averaging and phase fraction-weighted area averaging, respectively, which are defined for a parameter ψ as ψ A ψdA /A and ψ α k ψ k / α k , respectively.Here, A is the flow cross-sectional area, and α k is the void fraction of field-k.In the conventional one-dimensional two-fluid model, the momentum flux parameter is assumed to be unity, resulting in an illposed 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
Here, p, ρ k , Δ ṁ12 , M * i,k , χ, D c , D sm , and Φ k are the pressure, density of field-k, net intergroup mass transfer rate from group-1 to group-2 fluid particles, generalized drag force of field-k, 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-k, 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 ρ 1 is identical to ρ 2 , 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 M * i,k to the regularization of two-fluid model have attracted attention of many researchers 15, 16 .It was suggested that the inclusion of M * i,k , 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 M * i,k that could be applicable in any circumstances remained a challenge.Noting the argument that M * i,k usually tends to stabilize the two-fluid model, the analysis of the model hyperbolicity would be conservative if one adopts an algebraic form for M * i,k since the contribution of M * i,k to stabilizing the model will then not be credited in the characteristic analysis.

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 x α 1 , α 2 , v 1 , v 2 , v 3 , p, a i,1 , a i,2 is defined, then the governing equations of the system can be rewritten in the form of matrix as where A and B are the coefficient matrices.The detailed matrices of A , B , and C 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 | A λ − B | F λ 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 , F λ is rearranged and simplified to where the corresponding polynomials are given as

3.3
Note that F λ is replaced by g λ since four out of the eight roots for F λ 0 are guaranteed to be real roots.
Considering the fact that ρ 1 and ρ 3 are positive and the constrains of 1 − α 1 − α 2 being between zero and unity to maintain its physical meaning, we recognize that the coefficient of the highest order of g λ , that is, the 4th order, should be always positive.All roots of g λ are real if the following two requirements are satisfied: 1 g λ 0 has three real roots λ 1 , λ 2 , and 2 g λ 2 > 0, g λ 1 < 0, and g λ 3 < 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 g λ .The plots of g λ could therefore be given, from which it becomes straightforward to study graphically the hyperbolicity property of g λ .If there are four intersections between g λ and the real axis, then the characteristic equation g λ 0 will have four real roots for λ.More details are given in the next section.

Results and Discussions: A Case Study
To demonstrate the feasibility of the studies discussed above, g λ 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 3 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 3 at 20 • C. Local flow variables in 17 different flow conditions were measured at three  A higher value of R 2 always provides a better fit of the regression to the data set.Finally, as shown in Figure 1 c , C v.3 fits a 2nd-order polynomial of α quite well.As a result of the above analysis, three assumptions are made as follows:

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, v 1 v 2 v 3 v, 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, v g1 is approximately equal to v g2 , 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 g λ is simplified to a function of the void fraction of group-1 α 1 , void fraction of group-2 α 2 , characteristic root λ, and velocity v, only.In addition, if one introduces a new parameter, w as the dimensionless velocity, and defines it as w λ/v, then the characteristic equation g λ 0 converts to g w 0. The roots of g w 0 have the same mathematical properties as those of g λ 0, that is, if g w 0 has four real roots, so does the original equation g λ 0. Figures 2, 3, 4, 5, 6, 7, and 8 provide plots of g w versus w for 17 flow conditions.In these figures, subplot a is a plot in a relatively large range of w, while subplot b shows a plot in a sufficiently narrow range of w for the benefit of observing the two roots near which g w has small variations in the real axis.Figures 2-8 show clearly that g w has four positive roots, which are consistent with the physical meanings.As illustrated in the plots, the typical shape of g w is that g w has one pair of roots that are close to each other and are approaching a less-than-unity w value, whereas the other pair intersects at the values of w 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.
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, U r , is obtained as 21 where D d , g, C D , and ρ 3 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 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 As a result, the relative velocity is a function of the total void fraction αas shown in Figure 9.
Our study shows that the eigenvalue λ becomes complex for the one-dimensional twofluid 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 ν 3 as the y-axis and group-1 void fraction α 1 as the x-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.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 The stability criteria for S 1 and S 2 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 S 1 and S 2 shown in Figure 11 should be satisfied simultaneously.This analysis only provides the necessary condition of the well-posedness of the twofluid-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.

Conclusions
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 A , B , and C are given as follows:

Figure 2 :
Figure 2: g w for Runs 13, 17, and 18 over a a wide range of ω from 0 to 4 and b a narrow range of ω.

Figure 3 :
Figure 3: g w for Runs 14, 15, and 19 over a a wide range of ω from 0 to 3 and b a narrow range of ω.

Figure 4 :
Figure 4: g w for Runs 16, 21, and 29 over a a wide range of ω from 0 to 6 and b a narrow range of ω.

Figure 5 :
Figure 5: g w for Runs 20 and 27 over a a wide range of ω from 0 to 10 and b a narrow range of ω.

Figure 6 :
Figure 6: g w for Runs 22 and 26 over a a wide range of ω from 0 to 6 and b a narrow range of ω.

Figure 7 :
Figure 7: g w for Run 23 over a a wide range of ω from 0 to 3.5 and b a narrow range of ω.

Figure 8 :
Figure 8: g w for Runs 24, 25, and 28 over a a wide range of ω from 0 to 8 and b a narrow range of ω.

Figure 9 :Figure 10 :
Figure 9: Relative velocity versus the total void fraction α for the distorted and slug drops.

1 S 2 bFigure 11 :
Figure 11: Stability criteria of the one-dimensional two-fluid-IATE model with momentum flux parameters for a S 1 and b S 2 .