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

Authors: Xia Wang; Xiaodong Sun
Journal: Journal of Applied Mathematics
 (2010)

---

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


---

## Body

## 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 parameterCv,k is defined for field-k [20]:(2.1)Cv,k=〈αkvk2〉(〈αk〉〈〈vk〉〉2),
where vk 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 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 (2.2)∂α1∂t+v1∂α1∂z+α1∂v1∂z=-Δṁ12ρ1,∂α2∂t+v2∂α2∂z+α2∂v2∂z=Δṁ12ρ2,-∂α1∂t-v3∂α1∂z-∂α2∂t-v3∂α2∂z+(1-α1-α2)∂v3∂z=0,ρ1v1∂α1∂t+ρ1Cv,1v12∂α1∂z+α1ρ1∂v1∂t+2α1ρ1Cv,1v1∂v1∂z+α1∂p∂z=Mi,1*,ρ2v2∂α2∂t+ρ2Cv,2v22∂α2∂z+α2ρ2∂v2∂t+2α2ρ2Cv,2v2∂v2∂z+α2∂p∂z=Mi,2*,-ρ3v3∂α1∂t-ρ3Cv,3v32∂α1∂z-ρ3v3∂α2∂t-ρ3Cv,3v32∂α2∂z+(1-α1-α2)ρ3∂v3∂t+2(1-α1-α2)ρ3Cv,3v3∂v3∂z+(1-α1-α2)∂p∂z=Mi,3*,-[23-χ(DcDsm,1)2](ai,1α1)∂α1∂t-[23-χ(DcDsm,1)2](ai,1α1)v1∂α1∂z+[13+χ(DcDsm,1)2]ai,1∂v1∂z+∂ai,1∂t+v1∂ai,1∂z=Φ1,-χ(DcDsm,1)2(ai,1α1)∂α1∂t-χ(DcDsm,1)2(ai,1α1)v1∂α1∂z-χ(DcDsm,1)2ai,1∂v1∂z-(2ai,23α2)∂α2∂t-(2ai,23α2)v2∂α2∂z+13ai,2∂v2∂z+∂ai,2∂t+v2∂ai,2∂z=Φ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, p,ρk,Δṁ12,Mi,k*,χ,Dc,Dsm, 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 ρ1is 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 ofMi,k* to the regularization of two-fluid model have attracted attention of many researchers [15, 16]. It was suggested that the inclusion of Mi,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 Mi,k* that could be applicable in any circumstances remained a challenge. Noting the argument that Mi,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 Mi,k* since the contribution of Mi,k* 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 vectorx=(α1,α2,v1,v2,v3,p,ai,1,ai,2) is defined, then the governing equations of the system can be rewritten in the form of matrix as(3.1)[A]∂x∂t+[B]∂x∂z=[C],
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 (3.2)g(λ)=g1(λ)-g2(λ),
where the corresponding polynomials are given as (3.3)g1(λ)=ρ3g11(λ)g12(λ),g11(λ)=λ2-2Cv,3v3λ+Cv,3v32,g12(λ)=(α1+α2)λ2-2(α1Cv,2v2+α2Cv,1v1)λ+(α1Cv,2v22+α2Cv,1v12),g2(λ)=-(1-α1-α2)ρ1g21(λ)g22(λ),g21(λ)=λ2-2Cv,1v1λ+Cv,1v12,g22(λ)=λ2-2Cv,2v2λ+Cv,2v22.
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λ3 (let λ1>λ2>λ3);(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 ing(λ). 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.


## 4. 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/m3 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/m3 at 20°C. Local flow variables in 17 different flow conditions were measured at three measurement elevations along the flow directions with L/Dh=13, 67, and 121, respectively [9]. Table 1 lists the cross-sectional area-averaged data for each flow condition at L/Dh=13.Table 1
Area-averaged data for Runs 13–29 atL/Dh=13 [9]*.
Run no.α1 (%)α2 (%)vg1 (m/s)vg2 (m/s)vg,avg** (m/s)1324.965.980.360.370.371438.938.160.380.330.371532.6914.780.490.460.491620.410.110.450.16***0.451724.156.870.520.500.521826.1410.950.550.570.551919.9020.120.640.640.642015.840.140.620.07***0.612122.301.900.570.510.552220.173.390.750.770.762316.9917.390.850.820.842417.641.050.710.550.712514.371.780.840.900.852622.230.320.790.810.792713.130.141.580.41***1.582811.870.911.561.591.562919.272.291.691.841.70*: The run numbers are consistent with Vasavada [9].**: The averaged velocity was void fraction-weighted velocity of Therminol 59 provided by Vasavada [9].***: Due to the small sampling numbers of group-2 droplets, these velocity data are considered not statistically sound.The momentum flux parameterCv,k (k=1,2, and 3) are calculated according to the experimental data [9] when the total void fraction of the organic drops, α(α=α1+α2), is in the range between 0.1 and 0.6. The results are given in Figure 1, where Cv,k are plotted against α. From the plots, the relations between Cv,k and α could be determined. Figure 1(a) illustrates that Cv,1 is a relatively weak function of α. Most of Cv,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 Cv,1≃1.5. In Figure 1(b), although it shows a scattered distribution when α is small, Cv,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 Cv,2, such as a linear function or higher order polynomials, are also investigated, which generally provide smaller value of R2, 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 R2 always provides a better fit of the regression to the data set. Finally, as shown in Figure 1(c), Cv.3 fits a 2nd-order polynomial of α quite well.Cv,k versus the total void (drop) fraction α: (a) Cv,1, (b) Cv,2, and (c) Cv,3.
(a)(b)(c)As a result of the above analysis, three assumptions are made as follows:(4.1)(1)Cv,1≃1.5;(2)Cv,2=0.8824α-0.747;(3)Cv,3=-0.9028α2+0.0552α+1.1555.
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,v1=v2=v3=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, vg1 is approximately equal to vg2, 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.Figures2, 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.g(w) for Runs 13, 17, and 18 over (a) a wide range of ω from 0 to 4 and (b) a narrow range of ω.
(a)(b)g(w) for Runs 14, 15, and 19 over (a) a wide range of ω from 0 to 3 and (b) a narrow range of ω.
(a)(b)g(w) for Runs 16, 21, and 29 over (a) a wide range of ω from 0 to 6 and (b) a narrow range of ω.
(a)(b)g(w) for Runs 20 and 27 over (a) a wide range of ω from 0 to 10 and (b) a narrow range of ω.
(a)(b)g(w) for Runs 22 and 26 over (a) a wide range of ω from 0 to 6 and (b) a narrow range of ω.
(a)(b)g(w) for Run 23 over (a) a wide range of ω from 0 to 3.5 and (b) a narrow range of ω.
(a)(b)g(w) for Runs 24, 25, and 28 over (a) a wide range of ω from 0 to 8 and (b) a narrow range of ω.
(a)(b)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,Ur, is obtained as [21](4.2)Ur|Ur|=43DdΔρg(1-α)CDρ3,
where Dd, g, CD, 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](4.3)CD=2Dd3(1-α)gΔρσ,
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](4.4)CD=9.8(1-α)3.
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 ν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.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)S1=v1v3,S2=v2v3.
The stability criteria for S1 and S2 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 S1 and S2 shown in Figure 11 should be satisfied simultaneously.Stability criteria of the one-dimensional two-fluid-IATE model with momentum flux parameters for (a)S1 and (b) S2.
(a)(b)


## 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



---

*Source: 476839-2010-12-12.xml*