- Email: [email protected]
Surface tension effect on stability of two-phase stratified flow
Fluid Dynamics Research 39 (2007) 279 – 291
Surface tension effect on stability of two-phase stratified flow M.R. Ansari∗ , A. Eskandari Sani Mechanical Engineering Department, Tarbiat Modarres University, P.O. Box 14115-143, Tehran, Islamic Republic of Iran Received 19 March 2006; received in revised form 27 June 2006; accepted 29 June 2006 Communicated by A.D. Gilbert
Abstract In this paper, the effect of wavelength on the disturbed interface of a two-phase stratified regime will be discussed. It is assumed that the two phases are completely separated and are flowing parallel to each other along the duct. In order to analyze the stability of the two-phase stratified flow, a small disturbance is considered at the interface. Growth or decline of the interface disturbance, relative velocity and surface tension in time will be discussed analytically and numerically. Numerical analysis is applied using the finite difference method. The results of the numerical method are compared with the analytical results of prior investigators. The results of the present paper show that the non-linear conservation equations model discussed is capable of predicting two-phase stratified flow behavior. The surface tension increases the physics of the model and has a positive effect on the convergence of the numerical result of the two-phase flow model. The convergence of the numerical results was obtained and is presented in the results section of the paper. This conclusion is unique to this paper and the Ransom and Hicks model [Ransom, V.H., Hicks, D.L., 1984. Hyperbolic two-pressure models for two-phase flow. J. Comput. Phys. 53, 124–151]. © 2006 The Japan Society of Fluid Mechanics and Elsevier B.V. All rights reserved. Keywords: Two-phase stratified regime; Hydrodynamic instability; Interface; Surface tension
1. Introduction There are two views regarding two-phase flow modeling, from the perspectives of analytical and numerical methods. From the analytical framework, because of the complexity of the two-phase flow ∗ Corresponding author. Tel.: +98 21 8801 1001; fax: +98 21 8800 5040.
E-mail addresses: [email protected] (M.R. Ansari), [email protected] (A.E. Sani). 0169-5983/$32.00 © 2006 The Japan Society of Fluid Mechanics and Elsevier B.V. All rights reserved. doi:10.1016/j.fluiddyn.2006.06.002
280
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
governing equations for slug or mixed regimes, the solutions of the equations are not possible except in simple conditions. The detailed solution of the model for phase velocities, pressures and densities cannot be obtained, though a general or averaged condition can be calculated for the regime. The second method common in two-phase flow analysis is, instead of considering two fluids with different specification, a single phase proposed where every specification is the average of the related specifications of the fluids. Hence, the parameters used in the model are the average values of the two phases. This method is used when the general values of the whole system are required, as explained by Nguyen (1994). To overcome these difficulties, numerical methods are applied in order to solve the equations. Here, the flow parameters are studied on a section vertical to flow direction by Nguyen (1994), and Ramshaw et al. (1978). In this article, a two-phase stratified flow was studied in a horizontal duct. The two fluids, air and water are flowing parallel to each other in the duct. The effect of wavelength of the disturbance which is generated on the interface of fluids, relative velocity and surface tension were taken into consideration. The limit for initiation of instability on the interface is analyzed and compared with the results of previous investigators. The case is analyzed for the Kelvin–Helmholtz instability condition. The solutions of the non-linear conservation equations, including the surface tension effect, are demonstrated in this article. The present model predicts more physics with a higher stability limit in comparison with prior results. 2. Numerical analysis 2.1. Introduction A horizontal duct is considered with air and water that are flowing separately and concurrently. See Fig. 1. The condition is considered as a two-fluid model. In order to find the instability initiation from the stratified regime, a disturbance in the form of a sinusoidal wave with a very small amplitude was applied on the interface. Any wave growth (or damping) will be discussed with the effect of physical parameters. In order to examine the wave growth on the duct scale, the flow regime is considered to be onedimensional flow and parallel to the duct length. The duct height is very small in comparison with the duct length so the one-dimensional flow assumption could be established. The gravitational effect, viscosity, and mass transfer between phases were neglected. The liquid density was assumed to be constant and the duct scale is a function of wavelength. Any variation of the amplitude (increase or decrease) with respect to time was considered as wave growth (unstable condition) or wave damping (stable condition).
Fig. 1. Two-phase two fluid flow in a horizontal duct.
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
281
2.2. Governing equations The mass and momentum conservation equations for each fluid are considered. Because of the assumption of constant temperature and the non-phase change condition, the energy conservation equation was not taken into account. In this case, the mass equations for both fluids have been demonstrated by Ramshaw et al. (1978) j(1 1 )/jt + j(1 1 u1 )/jx = 0,
(1)
j(2 2 )/jt + j(2 2 u2 )/jx = 0.
(2)
The momentum equations for both fluids are as follows: j(1 1 u1 )/jt + j(1 1 u1 u1 )/jx + 1 jP1 /jx = 0,
(3)
j(2 2 u2 )/jt + j(2 2 u2 u2 )/jx + 2 jP2 /jx = 0.
(4)
In the above equations, is the void fraction, is the density, u is the velocity and P is the pressure where index 1 is for liquid and 2 for gas. Other equations that could be included are as follows: 1 + 2 = 1.
(5)
This is written in dimensionless form with respect to the duct height P1 − P2 = −/R,
(6)
where is the surface tension and R is the curvature radius at any point of the two-phase interface, White (2000). See Fig. 2. The value of R with respect to void fraction can be written as 1/R = −H (j2 1 /jx 2 ){1 + H 2 [j1 /jx]2 }−3/2 ,
Fig. 2. Pressure variation around curvature of interface.
(7)
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
282
Fig. 3. Mesh points configuration and its labeling.
where H is the duct height. As was assumed in the one-dimensional case, the value of 1/R can be approximated as [H j1 /jx]2 >1.
(8)
2.3. Solution of the governing equations In order to solve the governing equations, the finite difference mesh method was applied. A fully implicit numerical scheme was used to increase the numerical stability; see Ransom and Trapp (1983). The spatial points are fixed, but staggered in arrangement. The fluids’ properties, density, pressure and void fraction were calculated at the center of the mesh, and the fluids’ velocities were found on the mesh boundaries. The differential equations (1)–(4) are differentiated with respect to time and space using the above-mentioned method. The staggered mesh method application is shown in Fig. 3. The flux terms in the basic equations are differenced using donor-cell methods. The parameter SU is calculated at point j + 1/2 and defined as follows (Ansari et al., 1988): SU J +1/2 = 1/2[SJ +1 (UJ +1/2 − |UJ +1/2 |) + SJ (UJ +1/2 + |UJ +1/2 |)],
(9)
where S is a scalar and U is a vector. The derivatives are differentiated as follows: ∇J SU = [SU J +1/2 − SU J −1/2 ]/x.
(10)
3. The stability analysis The stability analysis of the two-phase flow equations was discussed by Ramshaw et al. (1978). For the stability condition, the critical wave number, kc , was defined as kc = |u1 − u2 |[1 2 /H (1 2 + 2 1 )]1/2 .
(11)
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
283
5 Sigma=0.05N/m Sigma=0.073N/m Sigma=0.1N/m Sigma=0.5N/m Sigma=1.0N/m
Critical Wave Length (m)
4
3
2
1
0
0
0.5
1
1.5
2
Relative Velocity (m/s)
Fig. 4. Critical wavelength change with relative velocity of two phases for different values of the surface tension.
From Eq. (11), it becomes clear that for wavelengths longer than 2/kc , the system is unstable, and stable for the shorter values. The instability at large wavelengths is called Kelvin–Helmholtz instability and is discussed by Kundu et al. (2002) and Mata and Pereyra (2002). See Fig. 4. It shows the variation of the critical wavelengths with relative velocity of two phases for different values of the surface tension. If k >kc (this is similar to = 0), for those wavelengths greater than 2/kc , the physical instability of Kelvin–Helmholtz type can be estimated from those equations with and without the inclusion of the surface tension parameter. Thus, a system of equations that contains complex characteristics would be capable of predicting physical instabilities of long waves. However, the result from short waves > 0 differs from = 0. For wavelengths shorter than 2/kc , the system is stable and for larger than 2/kc the system is unstable. 4. The initial and boundary conditions Periodic boundary conditions were used. This means that any value of parameters at outlet x = L were fed in at the inlet at x = 0. The duct scale is assumed as one wavelength. In other words, there are N mesh points on the duct scale. The condition at mesh point N + 1 is the same as at mesh point 1. See Fig. 3 for a description of the mesh. The initial value of velocities for both phases was changed in order to find the values for initiation of instability. A small amplitude perturbation was imposed on the interface as follows: 1 (x) = 10 + 0.005 sin(2x/),
(12)
where is the initial value of the wavelength and 10 is the liquid volume fraction at initial condition. The value of was varied from 0.1 to 1 m in order to study the wavelength effect on the stability or instability.
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
284
The initial values of phase velocities were defined as u1 (x) = u10 10 /1 (x),
(13)
u2 (x) = u20 (1 − 10 )/[1 − 1 (x)],
(14)
where u10 and u20 are the initial phase velocities.
5. Results 5.1. Effect of wavelength and initial relative velocity In order to study the wavelength effect on the interface between two phases, first a wavelength of 1 m length was assumed. The initial and boundary conditions are as follows: water = 1000 kg/m3 , gas = 1.3 kg/m3 , the surface tension between two phases = 0.0728 N/m, the initial uwater = 1 m/s and ugas = 3 m/s. The time duration for study = 0.2 s (the stability or instability initiation condition is considered), the initial value of void fraction = 0.5, the wavelength = 1 m, the duct length = 1 m and the duct height = 0.1 m. Fig. 4 is obtained analytically and shows the wavelength effect versus relative velocity with respect to surface tension. It shows that for the relative velocity of 2 m/s, the critical wavelength is 0.16 m. Thus, a wavelength higher than this value will grow with time. Hence, the growth of 1 m wavelength is certain and has been shown in Fig. 5. The result in Fig. 5 has been obtained numerically for relative velocity of 2 m/s with 1 m wavelength. 0.58 Time=0.0 s Time=0.2 s
0.56
Interface Shape (-)
0.54
0.52
0.5
0.48
0.46
0.44 0
0.25
0.5
0.75
1
Duct Length (m)
Fig. 5. Wave growth at the interface of stratified two-phase flow (relative velocity = 2 m/s and wavelength = 1 m).
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
285
0.58 Time = 0.0 sec Time = 0.2 sec
0.56
Interface Shape (-)
0.54
0.52
0.5
0.48
0.46
0.44 0
0.25
0.5
0.75
1
Duct Length (m)
Fig. 6. Wave damping at the interface of stratified two-phase flow (relative velocity = 0.35 m/s and wavelength = 1 m).
In the second step, all the parameters were kept as same as before, except that the relative velocity of the two-phase flow is decreased. The result indicates that growth was also decreased as expected. See Fig. 6. In the third step, in order to present the effect of wavelength, the parameters are kept as same as before except that the wavelength was decreased to 0.5 m and the relative velocity was taken into account at 10 m/s and 1.3 m/s, respectively. See Fig. 7. In the fourth step, the effect of the relative velocity was considered for a 1 m wavelength. The wave growth was evaluated as H (maximum − minimum ). See Fig. 8. The effect of different wavelengths on the growth factor is shown on Fig. 9 with respect to 2 m/s relative velocity (the growth factor definition is outlined in Section 5.4). 5.2. Surface tension effect In this section, the results of the two-phase flow equations, including the surface tension parameter, were compared with the same set of equations without the inclusion of the surface tension parameter. The time duration for wave growth was considered for 0.2 s. See Fig. 10 for comparison of the results that have been obtained with the same method (the differentiation of the equations and grid generation are same for both methods). In the first model, the surface tension equation was used as a coupling between two phases. In the second model, the hydrostatic pressure assumption declared by Ardron (1980) was applied as follows: P1 = Pi + 1/21 1 gH ,
(15)
P2 = Pi − 1/2(1 − 1 )2 gH ,
(16)
286
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291 0.6
Time = 0.0s Time = 0.2s, Relative velocity = 10m/s Time = 0.2s, Relative velocity = 0.3m/s
0.58
Interface Shape (-)
0.56
0.54
0.52
0.5
0.48
0.46
0.44 0
0.25
0.5
0.75
1
DuctLength(m)
Fig. 7. Wave growth or damping at the interface of stratified two-phase flow (wavelength = 0.5 m).
where P1 and P2 are the phase’s pressure, Pi is the interface pressure and H is the duct height. The initial and boundary condition for both models are the same. As is obvious from the figure, the growth for the model with surface tension is lower. This result is compared with the findings of previous investigators in Fig. 11. It should be noted that the numerical method which has been used by Ansari (2004) for the 5E2P model is the same as in the present work (the 5E2P model is explained by Ransom and Hicks (1984)). The data for DLY model was obtained from Physical Benchmarking Exercise (1987). 5.3. Mesh independency In order to verify the results obtained by the numerical method, the total number of mesh points was increased on the same duct scale. The results are presented in Fig. 12. As is obvious, when the number of mesh points was increased, convergence was obtained. The difference between the results that were obtained for total mesh point numbers 90 and 120 is negligible. 5.4. Void fraction effect One of the parameters which affects the two-phase flow behavior is the void fraction. In order to study this effect for a constant wavelength and constant relative velocity, the void fraction was changed during different iterations of computer runs. See Fig. 13. It demonstrates that as the void fraction increases, the growth (instability) decreases (the required pressure difference in order to lift the liquid from interface
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
287
16 Relative velocity = 0.01m/s Relative velocity = 0.1 m/s Relative velocity = 0.2 m/s Relative velocity = 2.0 m/s
14
Wave Growth (mm)
12
10
8
6
4
2 0
0.05
0.1
0.15
0.2
Time(s)
Fig. 8. Effect of relative velocity on wave growth or damping at the interface of stratified two-phase flow (wavelength = 1 m).
0.14 Wave length = 0.5m Wave length = 0.75m Wave length = 1m Wave length = 1.5m
Growth Factor (-)
0.13
0.12
0.11
0.1
0
0.025
0.05
0.075
0.1
Time (s)
Fig. 9. Effect of two-phase flow interfacial wavelength on growth factor with respect to time at relative velocity of 2 m/s.
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291 18 Sigma Model Interface Pressure Model
Wave Growth (mm)
16
14
12
10
8
0
0.05
0.1
0.15
0.2
Time (s)
Fig. 10. Comparison of two different models using same solution.
Present Model 5E2P Model DLY Model
40
Wave Growth (mm)
288
30
20
10
0
0
0.05
0.1
0.15
0.2
Time (s)
Fig. 11. Comparison of wave growth obtained by different two-phase flow models.
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291 12.5 Grid Number = 30 Grid Number = 60 Grid Number = 90 Grid Number = 120
Wave Growth (mm)
12
11.5
11
10.5
10 0
0.025
0.05
0.075
0.1
Time (s)
Fig. 12. Effect of total grid number on sigma model.
Wave Growth (mm)
11.3
11.2
11.1
0
0.2
0.4
0.6
0.8
1
Void Fraction (-)
Fig. 13. Effect of void fraction on wave growth at interface of stratified two-phase flow.
289
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
290
decreases). This was anticipated, and can be explained by the analytical analysis of the growth factor that varies as: Im() = 0
if k kc ,
Im() = k[1 2 (kc2 − k 2 )/(1 2 + 2 1 )]
if k < kc .
(17)
With respect to this equation, when k < kc , the wave growth will increase. On the other hand, as the void fraction 2 increases, because of higher density of the liquid in comparison with gas density, the growth factor will decrease. 6. Conclusion The non-linear conservation equations of two-phase stratified two-fluid flow including the surface tension parameter were solved by the finite difference technique. The governing equations were considered in transient form with one space dimension using a computer code created for this purpose. A wide range of results were generated in the horizontal duct, using air and water. The results indicated that: 1. As the initial wavelength was increased, the wave growth increased. 2. If the relative velocity increases from some particular value, the initial wave will become unstable. 3. The surface tension parameter has a positive effect on stabilizing flow instability (the result of other models without including the surface tension parameter were compared and presented in the results section). 4. As the liquid volume fraction decreases, the flow becomes more stable and the wave growth decreases. 5. The convergence of the results obtained from two-phase flow models based on hydrostatic pressure assumption (single-pressure) is not possible. The two-phase flow model including the surface tension parameter not only has more predictive capability, but also establishes convergent results. The results obtained by such a model could lead to better modeling of the behavior of two-phase flow systems and hence, to better system design and control. 6. The results obtained by two-pressure model (Ransom & Hicks, 5E2P) and the present work are consistent with one another. The single-pressure model (hydrostatic pressure assumption) result is over-predicted. References Ansari, M.R., 2004. Effect of pressure on two-phase stratified flow modeling. J. Nucl. Sci. Technol. 41 (7), 709–714. Ansari, M.R., Nariai, H., Hirano, M., Watanabe, T., 1988. Numerical analysis on slugging of air–water stratified flow in horizontal duct using MINCS code. Proceedings of the Third International Topical Meeting on Nuclear Power Thermal Hydraulics and Operations, Seoul, Korea, A1-115-122. Ardron, K.H., 1980. One-dimensional two-fluid equations for horizontal stratified two-phase flow. Int. J. Multiphase Flow 6, 295–304. Kundu, P.K., Cohen, I.M., 2002. Fluid Mechanics. Academic Press, New York. Mata, C., Pereyra, E., 2002. Stability of stratified gas–liquid flows. Development of Aerospace engineering and mechanics. University of Minnesota, AEM, 107 Akerman Hall, MN, USA. Nguyen, V.T., 1994. Two-phase, Gas–liquid Flow in Pipes—Mixture Properties and Equation of Channel. The Mitre Co., McLean, VA, USA.
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
291
Physical Benchmarking Exercise, 1987. DOE/EPRI, Second International Workshop of Two-phase Flow Fundamentals. Rensselar Polytechnic Institute, Troy, NY, USA. Ramshaw, J.D., Trapp, J.A., 1978. Characteristics, stability and short wavelength phenomena in two-phase flow equation systems. Nucl. Sci. Eng. 66, 93–102. Ransom, V.H., Hicks, D.L., 1984. Hyperbolic two-pressure models for two-phase flow. J. Comput. Phys. 53, 124–151. Ransom, V.H., Trapp, J.A., 1983. Applied mathematical methods in nuclear thermal hydraulic. In: Thermal Hydraulics of Nuclear Reactors, A.N.S, vol. 1, p. 99. White, F.M., 2000. Fluid Mechanics. Mc Graw-Hill, Inc., New York.
Surface tension effect on stability of two-phase stratified flow M.R. Ansari∗ , A. Eskandari Sani Mechanical Engineering Department, Tarbiat Modarres University, P.O. Box 14115-143, Tehran, Islamic Republic of Iran Received 19 March 2006; received in revised form 27 June 2006; accepted 29 June 2006 Communicated by A.D. Gilbert
Abstract In this paper, the effect of wavelength on the disturbed interface of a two-phase stratified regime will be discussed. It is assumed that the two phases are completely separated and are flowing parallel to each other along the duct. In order to analyze the stability of the two-phase stratified flow, a small disturbance is considered at the interface. Growth or decline of the interface disturbance, relative velocity and surface tension in time will be discussed analytically and numerically. Numerical analysis is applied using the finite difference method. The results of the numerical method are compared with the analytical results of prior investigators. The results of the present paper show that the non-linear conservation equations model discussed is capable of predicting two-phase stratified flow behavior. The surface tension increases the physics of the model and has a positive effect on the convergence of the numerical result of the two-phase flow model. The convergence of the numerical results was obtained and is presented in the results section of the paper. This conclusion is unique to this paper and the Ransom and Hicks model [Ransom, V.H., Hicks, D.L., 1984. Hyperbolic two-pressure models for two-phase flow. J. Comput. Phys. 53, 124–151]. © 2006 The Japan Society of Fluid Mechanics and Elsevier B.V. All rights reserved. Keywords: Two-phase stratified regime; Hydrodynamic instability; Interface; Surface tension
1. Introduction There are two views regarding two-phase flow modeling, from the perspectives of analytical and numerical methods. From the analytical framework, because of the complexity of the two-phase flow ∗ Corresponding author. Tel.: +98 21 8801 1001; fax: +98 21 8800 5040.
E-mail addresses: [email protected] (M.R. Ansari), [email protected] (A.E. Sani). 0169-5983/$32.00 © 2006 The Japan Society of Fluid Mechanics and Elsevier B.V. All rights reserved. doi:10.1016/j.fluiddyn.2006.06.002
280
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
governing equations for slug or mixed regimes, the solutions of the equations are not possible except in simple conditions. The detailed solution of the model for phase velocities, pressures and densities cannot be obtained, though a general or averaged condition can be calculated for the regime. The second method common in two-phase flow analysis is, instead of considering two fluids with different specification, a single phase proposed where every specification is the average of the related specifications of the fluids. Hence, the parameters used in the model are the average values of the two phases. This method is used when the general values of the whole system are required, as explained by Nguyen (1994). To overcome these difficulties, numerical methods are applied in order to solve the equations. Here, the flow parameters are studied on a section vertical to flow direction by Nguyen (1994), and Ramshaw et al. (1978). In this article, a two-phase stratified flow was studied in a horizontal duct. The two fluids, air and water are flowing parallel to each other in the duct. The effect of wavelength of the disturbance which is generated on the interface of fluids, relative velocity and surface tension were taken into consideration. The limit for initiation of instability on the interface is analyzed and compared with the results of previous investigators. The case is analyzed for the Kelvin–Helmholtz instability condition. The solutions of the non-linear conservation equations, including the surface tension effect, are demonstrated in this article. The present model predicts more physics with a higher stability limit in comparison with prior results. 2. Numerical analysis 2.1. Introduction A horizontal duct is considered with air and water that are flowing separately and concurrently. See Fig. 1. The condition is considered as a two-fluid model. In order to find the instability initiation from the stratified regime, a disturbance in the form of a sinusoidal wave with a very small amplitude was applied on the interface. Any wave growth (or damping) will be discussed with the effect of physical parameters. In order to examine the wave growth on the duct scale, the flow regime is considered to be onedimensional flow and parallel to the duct length. The duct height is very small in comparison with the duct length so the one-dimensional flow assumption could be established. The gravitational effect, viscosity, and mass transfer between phases were neglected. The liquid density was assumed to be constant and the duct scale is a function of wavelength. Any variation of the amplitude (increase or decrease) with respect to time was considered as wave growth (unstable condition) or wave damping (stable condition).
Fig. 1. Two-phase two fluid flow in a horizontal duct.
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
281
2.2. Governing equations The mass and momentum conservation equations for each fluid are considered. Because of the assumption of constant temperature and the non-phase change condition, the energy conservation equation was not taken into account. In this case, the mass equations for both fluids have been demonstrated by Ramshaw et al. (1978) j(1 1 )/jt + j(1 1 u1 )/jx = 0,
(1)
j(2 2 )/jt + j(2 2 u2 )/jx = 0.
(2)
The momentum equations for both fluids are as follows: j(1 1 u1 )/jt + j(1 1 u1 u1 )/jx + 1 jP1 /jx = 0,
(3)
j(2 2 u2 )/jt + j(2 2 u2 u2 )/jx + 2 jP2 /jx = 0.
(4)
In the above equations, is the void fraction, is the density, u is the velocity and P is the pressure where index 1 is for liquid and 2 for gas. Other equations that could be included are as follows: 1 + 2 = 1.
(5)
This is written in dimensionless form with respect to the duct height P1 − P2 = −/R,
(6)
where is the surface tension and R is the curvature radius at any point of the two-phase interface, White (2000). See Fig. 2. The value of R with respect to void fraction can be written as 1/R = −H (j2 1 /jx 2 ){1 + H 2 [j1 /jx]2 }−3/2 ,
Fig. 2. Pressure variation around curvature of interface.
(7)
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
282
Fig. 3. Mesh points configuration and its labeling.
where H is the duct height. As was assumed in the one-dimensional case, the value of 1/R can be approximated as [H j1 /jx]2 >1.
(8)
2.3. Solution of the governing equations In order to solve the governing equations, the finite difference mesh method was applied. A fully implicit numerical scheme was used to increase the numerical stability; see Ransom and Trapp (1983). The spatial points are fixed, but staggered in arrangement. The fluids’ properties, density, pressure and void fraction were calculated at the center of the mesh, and the fluids’ velocities were found on the mesh boundaries. The differential equations (1)–(4) are differentiated with respect to time and space using the above-mentioned method. The staggered mesh method application is shown in Fig. 3. The flux terms in the basic equations are differenced using donor-cell methods. The parameter SU is calculated at point j + 1/2 and defined as follows (Ansari et al., 1988): SU J +1/2 = 1/2[SJ +1 (UJ +1/2 − |UJ +1/2 |) + SJ (UJ +1/2 + |UJ +1/2 |)],
(9)
where S is a scalar and U is a vector. The derivatives are differentiated as follows: ∇J SU = [SU J +1/2 − SU J −1/2 ]/x.
(10)
3. The stability analysis The stability analysis of the two-phase flow equations was discussed by Ramshaw et al. (1978). For the stability condition, the critical wave number, kc , was defined as kc = |u1 − u2 |[1 2 /H (1 2 + 2 1 )]1/2 .
(11)
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
283
5 Sigma=0.05N/m Sigma=0.073N/m Sigma=0.1N/m Sigma=0.5N/m Sigma=1.0N/m
Critical Wave Length (m)
4
3
2
1
0
0
0.5
1
1.5
2
Relative Velocity (m/s)
Fig. 4. Critical wavelength change with relative velocity of two phases for different values of the surface tension.
From Eq. (11), it becomes clear that for wavelengths longer than 2/kc , the system is unstable, and stable for the shorter values. The instability at large wavelengths is called Kelvin–Helmholtz instability and is discussed by Kundu et al. (2002) and Mata and Pereyra (2002). See Fig. 4. It shows the variation of the critical wavelengths with relative velocity of two phases for different values of the surface tension. If k >kc (this is similar to = 0), for those wavelengths greater than 2/kc , the physical instability of Kelvin–Helmholtz type can be estimated from those equations with and without the inclusion of the surface tension parameter. Thus, a system of equations that contains complex characteristics would be capable of predicting physical instabilities of long waves. However, the result from short waves > 0 differs from = 0. For wavelengths shorter than 2/kc , the system is stable and for larger than 2/kc the system is unstable. 4. The initial and boundary conditions Periodic boundary conditions were used. This means that any value of parameters at outlet x = L were fed in at the inlet at x = 0. The duct scale is assumed as one wavelength. In other words, there are N mesh points on the duct scale. The condition at mesh point N + 1 is the same as at mesh point 1. See Fig. 3 for a description of the mesh. The initial value of velocities for both phases was changed in order to find the values for initiation of instability. A small amplitude perturbation was imposed on the interface as follows: 1 (x) = 10 + 0.005 sin(2x/),
(12)
where is the initial value of the wavelength and 10 is the liquid volume fraction at initial condition. The value of was varied from 0.1 to 1 m in order to study the wavelength effect on the stability or instability.
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
284
The initial values of phase velocities were defined as u1 (x) = u10 10 /1 (x),
(13)
u2 (x) = u20 (1 − 10 )/[1 − 1 (x)],
(14)
where u10 and u20 are the initial phase velocities.
5. Results 5.1. Effect of wavelength and initial relative velocity In order to study the wavelength effect on the interface between two phases, first a wavelength of 1 m length was assumed. The initial and boundary conditions are as follows: water = 1000 kg/m3 , gas = 1.3 kg/m3 , the surface tension between two phases = 0.0728 N/m, the initial uwater = 1 m/s and ugas = 3 m/s. The time duration for study = 0.2 s (the stability or instability initiation condition is considered), the initial value of void fraction = 0.5, the wavelength = 1 m, the duct length = 1 m and the duct height = 0.1 m. Fig. 4 is obtained analytically and shows the wavelength effect versus relative velocity with respect to surface tension. It shows that for the relative velocity of 2 m/s, the critical wavelength is 0.16 m. Thus, a wavelength higher than this value will grow with time. Hence, the growth of 1 m wavelength is certain and has been shown in Fig. 5. The result in Fig. 5 has been obtained numerically for relative velocity of 2 m/s with 1 m wavelength. 0.58 Time=0.0 s Time=0.2 s
0.56
Interface Shape (-)
0.54
0.52
0.5
0.48
0.46
0.44 0
0.25
0.5
0.75
1
Duct Length (m)
Fig. 5. Wave growth at the interface of stratified two-phase flow (relative velocity = 2 m/s and wavelength = 1 m).
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
285
0.58 Time = 0.0 sec Time = 0.2 sec
0.56
Interface Shape (-)
0.54
0.52
0.5
0.48
0.46
0.44 0
0.25
0.5
0.75
1
Duct Length (m)
Fig. 6. Wave damping at the interface of stratified two-phase flow (relative velocity = 0.35 m/s and wavelength = 1 m).
In the second step, all the parameters were kept as same as before, except that the relative velocity of the two-phase flow is decreased. The result indicates that growth was also decreased as expected. See Fig. 6. In the third step, in order to present the effect of wavelength, the parameters are kept as same as before except that the wavelength was decreased to 0.5 m and the relative velocity was taken into account at 10 m/s and 1.3 m/s, respectively. See Fig. 7. In the fourth step, the effect of the relative velocity was considered for a 1 m wavelength. The wave growth was evaluated as H (maximum − minimum ). See Fig. 8. The effect of different wavelengths on the growth factor is shown on Fig. 9 with respect to 2 m/s relative velocity (the growth factor definition is outlined in Section 5.4). 5.2. Surface tension effect In this section, the results of the two-phase flow equations, including the surface tension parameter, were compared with the same set of equations without the inclusion of the surface tension parameter. The time duration for wave growth was considered for 0.2 s. See Fig. 10 for comparison of the results that have been obtained with the same method (the differentiation of the equations and grid generation are same for both methods). In the first model, the surface tension equation was used as a coupling between two phases. In the second model, the hydrostatic pressure assumption declared by Ardron (1980) was applied as follows: P1 = Pi + 1/21 1 gH ,
(15)
P2 = Pi − 1/2(1 − 1 )2 gH ,
(16)
286
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291 0.6
Time = 0.0s Time = 0.2s, Relative velocity = 10m/s Time = 0.2s, Relative velocity = 0.3m/s
0.58
Interface Shape (-)
0.56
0.54
0.52
0.5
0.48
0.46
0.44 0
0.25
0.5
0.75
1
DuctLength(m)
Fig. 7. Wave growth or damping at the interface of stratified two-phase flow (wavelength = 0.5 m).
where P1 and P2 are the phase’s pressure, Pi is the interface pressure and H is the duct height. The initial and boundary condition for both models are the same. As is obvious from the figure, the growth for the model with surface tension is lower. This result is compared with the findings of previous investigators in Fig. 11. It should be noted that the numerical method which has been used by Ansari (2004) for the 5E2P model is the same as in the present work (the 5E2P model is explained by Ransom and Hicks (1984)). The data for DLY model was obtained from Physical Benchmarking Exercise (1987). 5.3. Mesh independency In order to verify the results obtained by the numerical method, the total number of mesh points was increased on the same duct scale. The results are presented in Fig. 12. As is obvious, when the number of mesh points was increased, convergence was obtained. The difference between the results that were obtained for total mesh point numbers 90 and 120 is negligible. 5.4. Void fraction effect One of the parameters which affects the two-phase flow behavior is the void fraction. In order to study this effect for a constant wavelength and constant relative velocity, the void fraction was changed during different iterations of computer runs. See Fig. 13. It demonstrates that as the void fraction increases, the growth (instability) decreases (the required pressure difference in order to lift the liquid from interface
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
287
16 Relative velocity = 0.01m/s Relative velocity = 0.1 m/s Relative velocity = 0.2 m/s Relative velocity = 2.0 m/s
14
Wave Growth (mm)
12
10
8
6
4
2 0
0.05
0.1
0.15
0.2
Time(s)
Fig. 8. Effect of relative velocity on wave growth or damping at the interface of stratified two-phase flow (wavelength = 1 m).
0.14 Wave length = 0.5m Wave length = 0.75m Wave length = 1m Wave length = 1.5m
Growth Factor (-)
0.13
0.12
0.11
0.1
0
0.025
0.05
0.075
0.1
Time (s)
Fig. 9. Effect of two-phase flow interfacial wavelength on growth factor with respect to time at relative velocity of 2 m/s.
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291 18 Sigma Model Interface Pressure Model
Wave Growth (mm)
16
14
12
10
8
0
0.05
0.1
0.15
0.2
Time (s)
Fig. 10. Comparison of two different models using same solution.
Present Model 5E2P Model DLY Model
40
Wave Growth (mm)
288
30
20
10
0
0
0.05
0.1
0.15
0.2
Time (s)
Fig. 11. Comparison of wave growth obtained by different two-phase flow models.
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291 12.5 Grid Number = 30 Grid Number = 60 Grid Number = 90 Grid Number = 120
Wave Growth (mm)
12
11.5
11
10.5
10 0
0.025
0.05
0.075
0.1
Time (s)
Fig. 12. Effect of total grid number on sigma model.
Wave Growth (mm)
11.3
11.2
11.1
0
0.2
0.4
0.6
0.8
1
Void Fraction (-)
Fig. 13. Effect of void fraction on wave growth at interface of stratified two-phase flow.
289
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
290
decreases). This was anticipated, and can be explained by the analytical analysis of the growth factor that varies as: Im() = 0
if k kc ,
Im() = k[1 2 (kc2 − k 2 )/(1 2 + 2 1 )]
if k < kc .
(17)
With respect to this equation, when k < kc , the wave growth will increase. On the other hand, as the void fraction 2 increases, because of higher density of the liquid in comparison with gas density, the growth factor will decrease. 6. Conclusion The non-linear conservation equations of two-phase stratified two-fluid flow including the surface tension parameter were solved by the finite difference technique. The governing equations were considered in transient form with one space dimension using a computer code created for this purpose. A wide range of results were generated in the horizontal duct, using air and water. The results indicated that: 1. As the initial wavelength was increased, the wave growth increased. 2. If the relative velocity increases from some particular value, the initial wave will become unstable. 3. The surface tension parameter has a positive effect on stabilizing flow instability (the result of other models without including the surface tension parameter were compared and presented in the results section). 4. As the liquid volume fraction decreases, the flow becomes more stable and the wave growth decreases. 5. The convergence of the results obtained from two-phase flow models based on hydrostatic pressure assumption (single-pressure) is not possible. The two-phase flow model including the surface tension parameter not only has more predictive capability, but also establishes convergent results. The results obtained by such a model could lead to better modeling of the behavior of two-phase flow systems and hence, to better system design and control. 6. The results obtained by two-pressure model (Ransom & Hicks, 5E2P) and the present work are consistent with one another. The single-pressure model (hydrostatic pressure assumption) result is over-predicted. References Ansari, M.R., 2004. Effect of pressure on two-phase stratified flow modeling. J. Nucl. Sci. Technol. 41 (7), 709–714. Ansari, M.R., Nariai, H., Hirano, M., Watanabe, T., 1988. Numerical analysis on slugging of air–water stratified flow in horizontal duct using MINCS code. Proceedings of the Third International Topical Meeting on Nuclear Power Thermal Hydraulics and Operations, Seoul, Korea, A1-115-122. Ardron, K.H., 1980. One-dimensional two-fluid equations for horizontal stratified two-phase flow. Int. J. Multiphase Flow 6, 295–304. Kundu, P.K., Cohen, I.M., 2002. Fluid Mechanics. Academic Press, New York. Mata, C., Pereyra, E., 2002. Stability of stratified gas–liquid flows. Development of Aerospace engineering and mechanics. University of Minnesota, AEM, 107 Akerman Hall, MN, USA. Nguyen, V.T., 1994. Two-phase, Gas–liquid Flow in Pipes—Mixture Properties and Equation of Channel. The Mitre Co., McLean, VA, USA.
M.R. Ansari, A.E. Sani / Fluid Dynamics Research 39 (2007) 279 – 291
291
Physical Benchmarking Exercise, 1987. DOE/EPRI, Second International Workshop of Two-phase Flow Fundamentals. Rensselar Polytechnic Institute, Troy, NY, USA. Ramshaw, J.D., Trapp, J.A., 1978. Characteristics, stability and short wavelength phenomena in two-phase flow equation systems. Nucl. Sci. Eng. 66, 93–102. Ransom, V.H., Hicks, D.L., 1984. Hyperbolic two-pressure models for two-phase flow. J. Comput. Phys. 53, 124–151. Ransom, V.H., Trapp, J.A., 1983. Applied mathematical methods in nuclear thermal hydraulic. In: Thermal Hydraulics of Nuclear Reactors, A.N.S, vol. 1, p. 99. White, F.M., 2000. Fluid Mechanics. Mc Graw-Hill, Inc., New York.