- Email: [email protected]
Numerical solution of hyperbolic two-fluid two-phase flow model with non-reflecting boundary conditions
International Journal of Engineering Science 40 (2002) 789–803 www.elsevier.com/locate/ijengsci
Numerical solution of hyperbolic two-fluid two-phase flow model with non-reflecting boundary conditions Moon-Sun Chung a
a,*
, Keun-Shik Chang b, Sung-Jae Lee
a
Thermal-Hydraulics Safety Research Team, Korea Atomic Energy Research Institute, 150 Dukjin-Dong, Yusong-Gu, Taejon 305-353, South Korea b Department of Aerospace Engineering, Korea Advanced Institute of Science and Technology, 373-1 Kusong-Dong, Yusong-Gu, Taejon 305-701, South Korea Received 21 July 2000; received in revised form 11 May 2001; accepted 9 July 2001 (Communicated by J.T. ODEN)
Abstract Flux vector splitting method is applied to the two-fluid six-equation model of two-phase flow, which takes account of surface tension effect via the interfacial pressure jump terms in the momentum equations. The latter terms using the concept of finite-thickness interface are derived as a simple function of fluid bulk moduli. We proved that the governing equation system is hyperbolic with real eigenvalues in the bubbly, slug, and annular flow regimes. The governing equations do not need any conventional artificial stabilizing terms like the virtual mass terms. Sonic speeds obtained through characteristic analysis show excellent agreement with the existing experimental data. Edwards pipe problem is solved as a benchmark test of the present two-phase equation model. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Two-fluid model; Surface tension; Hyperbolic system; Sonic speed; Flux vector splitting; Edwards pipe
1. Introduction The non-equilibrium exchange of mass, momentum, and energy between the gas and the liquid phases can be conveniently accounted for by using the two-fluid models. However, the common form of two-fluid models assuming single pressure across the interface gives complex eigenvalues
*
Corresponding author. Tel.: +82-42-868-2895. E-mail address: [email protected] (M.-S. Chung).
0020-7225/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 2 2 5 ( 0 1 ) 0 0 0 9 2 - 1
790
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
Nomenclature A interfacial area A; B coefficient matrices C speed of sound E source matrix ¼ A1 C F flux vector h enthalpy per unit mass J; K Jacobian matrices L fluid bulk modulus p pressure Q volumetric source term for heat R averaged bubble radius t time u internal energy per unit mass U primitive variable vector v flow velocity V control volume x space coordinate Greek symbols a volumetric phase concentration b relative surface thickness d surface thickness e eigenvalue switching factor / source term k eigenvalue K eigenvalue matrix q fluid density r surface tension Subscripts and superscripts g gas phase i interface j index for coordinate k index for each fluid l liquid phase m mth component n nth time-step, current time level p constant pressure s saturation condition w on the wall )1 inverse of a matrix first-order flux
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
791
and consequent numerical instability because the initial value problem is ill-posed, see [1–3]. In contrast, the models assuming two pressures across the interface present hyperbolic type governing equations with real characteristics. However, the existing two-pressure model represents only particular flow type like a stratified flow, see [4]. It has been known that mathematical property of the governing equations is improved by introducing physical terms in the governing equations. Ramshaw and Trapp [2], for example, added the surface tension equation to the governing equation system and Travis et al. [5] took account of the viscous stresses in the momentum equations. Stuhmiller [6] argued that when the interfacial pressure force terms were taken into account, the characteristic speed of the void wave became partially real. He analyzed that to assure real characteristics for the flow having unequal phase velocities, the interfacial pressure had to be lower than the bulk pressure by an amount proportional to the square of the relative velocity. Rousseau and Ferch [7], on the other hand, derived a condition for the hyperbolic governing equations in terms of static pressure difference between phases. They showed that the condition is identical to that of Kelvin–Helmholtz instability against the long wavelengths. Ardron [8] examined the one-dimensional two-fluid equations for a stratified flow, ignoring viscosity but retaining gravity and surface tension effect. The equations were observed stable over a realistic range of conditions, producing wave velocities in good agreement with the exact solution. A new promising approach of eliminating the instability was proposed by the present authors and their coworkers: The interfacial pressure jump terms are introduced based on the surface tension in the two-fluid momentum equations, see [9–11]. The system of equations manifested real characteristics in all the bubbly, slug, and annular flow regimes when the interfacial pressure jump terms are expressed as a product of effective bulk moduli and the gradient of interfacial area density. To avoid the difficulty of mathematically treating the interfacial area transport equation, we developed an expression of the interfacial pressure jump terms dependent only upon the fluid bulk moduli. In Section 2, we will recapitulate how the surface tension terms are brought about in the governing equations. Finite-volume flux vector splitting (FVS) method is formulated in order to solve the system of equations numerically in Section 3. Unlike the conventional methods of twophase flow, it is less dissipative and captures sharp gradient of flows very well without oscillation. In Section 4, we discuss the numerical solution of Edwards pipe problem as a benchmark test of transient two-phase flow.
2. Hyperbolic two-fluid model The conservation equation system consists of mass, momentum, and energy conservation equations by which the area-averaged phasic properties are related one another. The one-dimensional six equations take the following form: Continuity: oðak qk Þ oðak qk mk Þ þ ¼ /c;k ; ot ox
ð1Þ
792
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
Momentum: oðak qk mk Þ oðak qk m2k Þ opk oak þ ðpk pi Þ ¼ /c;k vi þ Fi;k þ Fw;k þ Bk ¼ /e;k ; þ þ ak ox ox ot ox
ð2Þ
Internal energy: oðak qk uk Þ oðak qk vk uk Þ oak oðak vk Þ þ þ pk ¼ Qi;k þ Fi;k vi þ /c;k hk þ Fw;k vk ¼ /e;k ; þ ot ox ox ot
ð3Þ
where ak ; qk ; pk ; vk ; uk , and hk denote void fraction, density, pressure, velocity, internal energy, and enthalpy of the phase k, respectively, and Bk is a gravitational body force. The source functions, /c;k ; /m;k , and /e;k , are made of algebraic constitutive relations, as shall be mentioned at the end of this section. We assign k ¼ g for the gas and k ¼ l for the liquid. The interfacial pressure jump term, ðpk pi Þðoak =oxÞ, which will be expressed as the product of fluid bulk modulus and an infinitesimal variable, is very small relative to other terms in the momentum equation (2). It shall be shown in Section 4 by the characteristic analysis that these terms make the equation system hyperbolic. 2.1. Interfacial pressure jump terms Young and Laplace proposed a well-known surface tension: pg pl ¼
2r ; R
ð4Þ
where r is the surface tension and R is the averaged bubble radius. We now assume a finite interfacial thickness d between the two radii Rg and Rl shown in Fig. 1. For an imaginary sphere having a radius Ri , where Rg þ Rl ¼ 2Ri , we rewrite Eq. (4) as
Fig. 1. Hypothetical sphere at Ri (solid circle) with an infinitesimal thickness d.
pg pl ¼
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
793
r 2d 2d r ¼ : Rg þ d=2 d Rl d=2 d
ð5Þ
For small dilation of the bubble marked by the radial increment DRi and surface area change DAi , we can expect corresponding volume change interior and exterior to the film, namely, DVg and DVl . From the bubble geometry, we can easily derive the following relations:
d=2 ; Rg þ d=2 Rl DAi d=2 ¼ 1þ : Rl d=2 2 DVl Rg DAi ¼ 2 DVg
1
ð6Þ ð7Þ
We can rewrite Eq. (5) as 4r Rg DAi 4r Rl DAi ¼ : 1 1þ pg pl ¼ d d 2 DVg 2 DVl
ð8Þ
It is recognized that the relation, L ¼ cr=d, is used in the physical chemistry and statistical mechanics, where L is a bulk modulus and c is a constant, see [12,13]. We assume that the surface tension stress, 4r=d ¼ L, which plays a role of Lagrangian multiplier as explained by Aubin and Ekeland [14], can be replaced by the phasic bulk moduli added together: 4r ¼ Lg þ Ll ¼ qg Cg2 þ ql Cl2 ; d
ð9Þ
where Ck is the sonic speed. The pressure jump is split into the phasic components as shown in [9–11]: Rg DAi Rl DAi pg pl ¼ ðpg pi Þ þ ðpi pl Þ ¼ Lg 1 Ll 1 þ : ð10Þ 2 DVg 2 DVl Here, Lg and Ll are the bulk modulus of the gas and the liquid phase, respectively, and pi is the interfacial pressure on the imaginary sphere Ri . The pressure jump interior and exterior to the sphere Ri has the following fraction as derived from Eqs. (8)–(10), pg pi Lg ¼ ; pg pl Lg þ Ll
pl pi Ll ¼ : p g p l Lg þ Ll
ð11Þ
Consequently, the interfacial pressure jump terms in the momentum equation (2) have the expression ðpk pi Þ
oak 2r Lk oak oak ¼ ð1Þn ¼ ð1Þn bLk ; Ri Lg þ Ll ox ox ox
ð12Þ
794
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
where the exponent n stands for the liquid if n ¼ 1 and for the gas if n ¼ 2. The parameter b named ‘relative surface thickness’ becomes by using Eq. (9), b
2r 1 d ¼ : Ri Lg þ Ll 2Ri
ð13Þ
Therefore, the interfacial pressure jump term is represented by the product of bulk modulus Lk and the relative surface thickness b which is a function of temperature and pressure. For a bubbly flow assumed of a perfect mixture, the mixture bulk modulus is Lm ¼ V
dp dp dp ¼V : ¼ V dV dVg þ dVl Vg dp=Lg þ Vl dp=Ll
ð14Þ
Since it holds that Lg Ll , Eq. (14) yields Lm Lg =ag , which is true in the void fraction range ag P Lg =Ll . We assume here that the mixture bulk modulus is equal to the bulk modulus of the gas by taking ag ¼ Oð1Þ, namely, Lm ¼ qg Cg2 . The resultant surface thickness evaluated by Eq. (9) becomes d ¼ 2r=Lg , which is of the order of magnitude 106 –107 m as shown in Fig. 2. It is equivalent to the initial bubble radius in homogeneous nucleation, see [15]. For the slug flow, the effective bulk modulus Lk can be obtained from a simplified physical model using the assumption of no elastic interaction between the phases. Because a wave traveling in one phase is not disturbed by the presence of the other phase in this model, the effective bulk modulus of one phase is same as the bulk modulus of its sole phase, namely, Lg ¼ qg Cg2 and Ll ¼ ql Cl2 . In the annular flow, the pressure wave in the gas phase is not transmitted to the liquid phase but is mostly either reflected back or changed into capillary waves on the liquid surface. We can then assume Lg ¼ qg Cg2 and Ll ¼ 0. It is suggested by Van Stralen [16] that the minimum radius of the bubble nuclei can be expressed in the nucleation process as a function of saturation temperature, latent heat, and
Fig. 2. Bulk modulus and surface thickness of the water–vapor mixture.
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
795
superheating temperature. It is suggested that the minimum radius of bubble nuclei be of the order of magnitude 106 –107 m for many liquids. Therefore we can deduce that the radius of curvature Ri has the same order of magnitude as the surface thickness d. Eq. (5) suggests that for Rg P 0 and Rl Rg ¼ d, we must have the inequality 0 < b 6 1. For the bubbly flow, since the surface thickness and the bubble diameter have the same order of magnitude in the limit, we are led to the approximation, b ¼ Oð1Þ. In contrast, for the annular flow with an interface having a large radius of curvature, b becomes relatively small. On this reason, we can make interpolation of b as a function of void fraction. For the convenience, we treat this parameter as a constant in the characteristic analysis. 2.2. Transformation of governing equations Using the definition Xk ðoqk =opk Þuk , Yk ðoqk =ouk Þpk , and the identity opg =ox ¼ opl =ox which is derived from Eq. (4) for the equilibrium states of the bubble, we can write out Eqs. (1)–(3) as follows: Continuity: oag opg oag opg ovg oug oug þ ag Xg þ qg vg þ ag Xg vg þ qg ag þ ag Yg þ vg ¼ /c;g ; qg ot ot ox ox ox ot ox oal opg oal opg ovl oul oul ql þ al Xl þ ql vl þ al Xl vl þ ql al þ al Yl þ vl ¼ /c;l : ot ot ox ox ox ot ox
ð15Þ ð16Þ
Momentum: ag qg al ql
ovg opg ovg oag þ ag þ ag qg vg þ bLg ¼ /m;g ; ot ox ox ox
ovl opg ovl oal þ al þ al ql vl bLl ¼ /m;l : ot ox ox ox
ð17Þ ð18Þ
Internal energy: ag qg al ql
oug oug oag oag ovg þ ag qg v g þ pg þ pg v g þ pg ag ¼ /e;g ; ot ox ot ox ox
oul oul oal oal ovl þ al ql vl þ pl þ pl vl þ pl al ¼ /e;l : ot ox ot ox ox
ð19Þ ð20Þ
It holds that ag þ al ¼ 1. If we use an auxiliary thermodynamic equation dqk ¼
oq opk
dpk þ
uk
oq ouk
duk ¼ Xk dpk þ Yk duk ; Pk
ð21Þ
796
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
then a matrix equation holds A
oU oU þB ¼ C; ot ox
ð22Þ
where U ¼ ½ ag
pg
vg
vl
ug
ul T ;
ð23Þ
T C ¼ /c;g /c;l /m;g /m;l /e;g /e;l ; 3 2 0 0 ag Yg 0 qg ag Xg 6 ql al Xl 0 0 0 al Yl 7 7 6 7 6 0 0 a q 0 0 0 g g 7; A¼6 6 0 0 0 7 0 0 al ql 7 6 4 pg 0 0 0 ag qg 0 5 0 0 0 0 al ql pl 2 qg vg ag vg Xg ag qg 0 ag Yg vg 6 ql vl al vl Xl 0 al ql 0 6 6 bLg a a q v 0 0 g g g g B¼6 6 bLl a 0 a q v 0 l l l l 6 4 pg vg 0 ag pg 0 ag qg vg 0 0 al pl 0 pl vl
ð24Þ
ð25Þ
3 0 al Yl vl 7 7 0 7 7: 0 7 7 0 5 al ql vl
ð26Þ
The source vector C in Eq. (22) requires constitutive equations related with the transport of mass, momentum, and energy between phases. To maintain the consistency with the interfacial transfer relations, we prefer the same constitutive equations as were adopted in RELAP5/MOD3 code for the bubbly, slug, and annular flow regimes. In practice, we identify flow transition regimes like the bubble-to-slug and the slug-to-annular flow regimes by interpolation as a function of void fraction. Eq. (22) can be transformed into oU oU þG ¼ E; ot ox
ð27Þ
where G ¼ A1 B and E ¼ A1 C. The governing equations are then straightforward: oV oV þH ¼ M; ot ox
ð28Þ
where H ¼ J G J 1 , M ¼ J E, and J ¼ oV=oU. The conservation vector V is V ¼ ag qg
al ql
ag qg v g
al ql vl
ag qg ug
al ql ul
T
:
ð29Þ
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
We define the flux vector F by 3 2 ag qg vg 7 6 al ql vl 7 6 2 6 ag qg vg þ pg ag þ bLg ag 7 7; 6 F¼6 2 7 6 al ql vl þ pl al þ bLl al 7 5 4 ag qg vg ug al ql vl ul
797
ð30Þ
where the relative surface thickness b is a relatively infinitesimal variable. An alternative form of Eq. (28) is oV oF oF þ þW ¼ M; ot ox ox
ð31Þ
where W ¼ ðJ G KÞ K 1 and K is another Jacobian matrix oF=oU. To obtain the flux vector F, Eq. (31) is further transformed to oF oF þR ¼ N; ot ox where R ¼ K G 2 qg 6 ql 6 6 qg vg J ¼6 6 ql vl 6 4 qg ug ql ul
ð32Þ
K 1 , N ¼ K E, and the Jacobian matrices J, K are 3 ag Xg 0 0 ag Yg 0 7 al Xl 0 0 0 al Yl 7 7 ag vg Xg ag qg 0 ag vg Yg 0 7; 7 al vl Xl 0 al ql 0 al vl Yl 7 5 ag ug Xg 0 0 ag ðug Yg þ qg Þ 0 al ul Xl 0 0 0 al ðul Yl þ ql Þ
2
qg vg 6 ql vl 6 6 qg v2g þ pg þ bLg K ¼6 6 q v2 pl bLl 6 l l 4 qg vg ug ql vl ul
ag vg Xg al vl Xl ag v2g Xg þ ag al v2l Xl þ al ag vg ug Xg al vl ul Xl
ag qg 0 2ag qg vg 0 ag qg ug 0
0 al ql 0 2al ql vl 0 al ql ul
ag v g Yg 0 ag v2g Yg 0 ag vg ðug Yg þ qg Þ 0
ð33Þ
3 0 7 al vl Yl 7 7 0 7: 2 7 al vl Yl 7 5 0 al vl ðul Yl þ ql Þ ð34Þ
3. Flux vector splitting The fluxes normal to the cell interfaces are calculated from the Riemann problem, which is shown in the work of Stadtke et al. [17]. For this purpose, the governing equation (27) is projected normal to the cell interface:
798
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
oU oU þ Gn ¼ En: ot on
ð35Þ
Corresponding to the six real eigenvalues, a linearly independent set of eigenvectors exists which is a necessary condition for a hyperbolic system of equations. A similarity transformation exists by which the characteristic form of the equations is obtained: L1 n
oU oU þ Kn L1 ¼ L1 n n En; ot on
ð36Þ
where the eigenvalue matrix is Kn ¼ L1 n G n Ln . The columns of matrix Ln are the right eigenvectors of the coefficient matrix G n and the rows of L1 n are the left eigenvectors of G n . Using the individual eigenvalues, the coefficient matrix G n can be split into the elementary parts Gn ¼
6 X
G n;m ;
ð37Þ
m¼1
with G n;m ¼ Ln Kn;m L1 n . The diagonal matrix Kn;m consists of the mth eigenvalue only. To obtain flux vectors at the cell interfaces, Eq. (32) is used as oF oF þ Rn ¼ N n; ot on
ð38Þ
where the coefficient matrix Rn ¼ K n G n K 1 n with K n ¼ oFn =oU. The similarity transformation does not change the eigenvalues of the governing equation system. Therefore, the characteristic form of the governing equations, equivalent to Eq. (36), is L01 n
oF oF þ Kn L01 ¼ L01 n n N n; ot on
ð39Þ
where the transformed right eigenvectors are L01 ¼ K n Ln . n If we linearize Eq. (39), then the numerical flux at the cell interface is composed of the positive and negative contributions depending on the sign of eigenvalues as X Rm X Rm Fj þ Fjþ1 ; ð40Þ Fjþ1=2 ¼ km j km jþ1 m;km >0 m;km <0 where km denotes the mth component of eigenvalue matrix. The resultant form of the governing equations based on the finite volume method is obtained from the conservative equation (31) as Vjnþ1 ¼ Vjn
Dt Dt Fjþ1=2 Fj1=2 Wjn Fjþ1=2 Fj1=2 þ DtMjn : Dx Dx
ð41Þ
A non-physical discontinuity might appear at the sonic transition of the computed results. It can be corrected by a simple modification similar to that used by Buning and Steger [18]. It redefines the eigenvalues k m using a small parameter em , namely,
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
k m ¼ km
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 þ ðem =km Þ2 2
799
ð42Þ
:
We assign a very small number to em , for example, of the order ðem =km Þ2 6 Oð104 Þ, to make modification of the eigenvalues truly small. 4. System eigenvalues and numerical example 4.1. Characteristic analysis The eigenvalues of the equation system represent the propagation speeds of small-amplitude, short-wavelength perturbations as indicated by Whitham [19]. Whereas both source and dispersion terms of the source vector play an important role for long-wavelength disturbances and the nonlinear wave interaction causes dominant effect on the large-amplitude disturbances. If the eigenvalues are all real, the equation system is hyperbolic and its solutions are stable against small disturbances. The eigenvalues of coefficient matrix G in Eq. (27) are obtained by DetðG kIÞ ¼ 0:
ð43Þ
Here, we assume that the ‘relative surface thickness’ is a constant, 1.0, in order to derive the eigenvalues analytically from Eq. (43). As a result, we can derive a sixth-order polynomial equation as follows: P6 ðkÞ ¼ ðk vg Þðk vl Þðk4 þ Z1 k3 þ Z2 k2 þ Z3 k þ Z4 Þ ¼ 0;
ð44Þ
where the coefficients are Z1 ¼ 2ðvg þ vl Þ, Z2 ¼
Z3 ¼
1 þ ag Cl2 ql # ( " ) L L 2 2 g l al Cg2 qg vg þ vl þ 2vg vl Cl2 þ ag Cl2 ql vg þ vl þ 2vg vl Cg2 ; qg ql al Cg2 qg
2 þ ag Cl2 ql !# ( " ) L L g l v2g þ ag Cl2 ql vl Cg2 v2g þ vg al Cg2 qg vg Cl2 v2l þ vl v2l ; qg ql al Cg2 qg
1 Z4 ¼ 2 al Cg qg þ ag Cl2 ql
( al Cg2 qg
"
Lg v2g Cl2 v2l qg
!# þ ag Cl2 ql
Cg2 v2g
L
l
ql
) v2l :
800
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
Table 1 Three sets of six eigenvalues Eigenvalues Bubbly regime
Slug regime Annular regime
k1;2 ¼ vg ; vl ;
k3;4 ¼ vg Cg ;
k5;6 ¼ vl Cl
k1;2 ¼ vg ; vl ;
k3;4 ¼ vg Cg ;
k5;6 ¼ vl Cl
k1;2 ¼ vg ; vl ;
k3;4 ¼ vg Cg ;
k5;6 ¼ vl Cl
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qg Cg2 2 al Cg qg þ ag Cl2 ql
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi al qg Cg2 al Cg2 qg þ ag Cl2 ql
The analytical solution of the characteristic equation (44) gives three sets of six real eigenvalues as listed in Table 1, which represent the bubbly, slug, and annular flows. The first two eigenvalues, k1;2 ¼ vg ; vl , represent the convection velocity of the gas and the liquid phases, respectively. The other two eigenvalues, k3 and k5 , represent approximately the sonic speeds in the gas and the liquid phases. For the bubbly and slug flows, the total sonic speed can be obtained by the void fraction weighting [10,11]: Ct ¼
k3 k5 : al k3 þ ag k5
ð45Þ
For the annular flow, the individual phasic sonic speeds, k3 and k5 , shall be used. In Fig. 3, we compared the computed total and phasic sonic speeds with the experimental data produced by Henry et al. [20]. For the bubbly flow, the total sonic speed agrees reasonably well with the experimental data in the void fraction range 0 < ag < 0:2 as shown in Fig. 3(a). The increasing deviation in the range ag > 0:2 is probably caused by transition of the flow regimes. Fig. 3(b) shows that the sonic speed of the water–air slug flow is in good agreement between the computed and experimental data in the entire void fraction range 0 < ag < 1. The sonic speed of the gas phase of the annular flow agrees well in Fig. 3(c) with the experimental data [20]. For the liquid phase, unfortunately, experimental data do not exist. The computed result shows that sonic speed of the liquid phase is subject to a rapid initial decrease for the very low the void fraction range, due to the effect of increasing elasticity at the interface. 4.2. Edwards pipe problem This problem is the benchmark test for Edwards and O’Brien’s [21] blowdown experiment. It represents a loss-of-coolant accident (LOCA) problem in the pressurized water reactor. The pipe is horizontal, 4-m long with the cross-section area 0:00456 m2 . It is filled with subcooled water under the initial pressure 7.0 MPa. Discharge of the flow is abruptly initiated by suddenly opening one end of the pipe that has a narrow cross-section area 0:00397 m2 . An expansion wave sweeps backward from the open end of the pipe, accompanied by flashing of the subcooled water due to the severe depressurization.
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
801
Fig. 3. Calculated sonic speeds compared with the experimental data: (a) Bubbly flow (pg ¼ 283 kPa). (b) Slug flow (pg ¼ 172 kPa). (c) Annular flow (pg ¼ 100 kPa).
In Fig. 4(a), the pressure is predicted at the open end of the pipe; the pipe is composed of 20 cells. The calculated time-dependent vapor pressure lies close to the measured pressure data. After the expansion wave is reflected at the closed end of the pipe, which is observed by the slight undershoot of the pressure in Fig. 4(a), the pressure is maintained near 2.8 MPa for about 0.15 s. The result of the present model shows that the pressure undershooting is less and the second depressurization occurs somewhat earlier in the calculated data than in the experimental data. However, the overall comparison appears as good as RELAP5/MOD3 code does. Fig. 4(b) shows the calculated time-dependent void fraction at the midsection of the pipe. Significant discrepancy which might have been originated from the inaccurate constitutive equations is apparent in the initial stage of the flow development, up to t ¼ 0:10 s. This deviation of the void fraction is also evident with the RELAP5/MOD3 code. Nevertheless, comparison of the present results with the experiment is significantly improved in later times. In contrast,
802
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
(a)
(b)
Fig. 4. Edwards pipe blowdown problem: (a) pressure history at the open end of the pipe; (b) void-fraction history at the midsection of the Edwards pipe.
RELAP5/MOD3 shows irregular void-fraction change up to t ¼ 0:25 s. The present result is also smoother than the earlier data of the present authors [10]. 5. Conclusions A hyperbolic system of two-phase two-fluid conservation laws has been derived using the interfacial pressure jump terms expressed ultimately as a function of fluid bulk moduli. Owing to its hyperbolicity, the equation system can be solved by an upwind numerical method like the flux vector splitting method. All the eigenvalues of the equation system turned out to be real which must be a remarkable improvement over the earlier two-phase formulations. The computed sonic speeds gave good comparison with the measured data. It has been demonstrated that the present model can be efficiently applied to the wave-dominant two-phase flows having strong expansion and shock waves with abrupt depressurization and flashing phenomena such as Edwards pipe problem.
Acknowledgements The authors are grateful to Dr. Won-Jae Lee in KAERI for the helpful technical discussions he offered while we prepare this paper. References [1] R.W. Lyczkowski, D. Gidaspow, C.W. Solbrig, E.D. Hughes, Characteristics and stability analyses of transient one-dimensional two-phase flow equations and their finite difference approximations, Nucl. Sci. Eng. 66 (1978) 378–396.
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
803
[2] J.D. Ramshaw, J.A. Trapp, Characteristics, stability, and short-wave length phenomena in two-phase flow equation systems, Nucl. Sci. Eng. 66 (1978) 93–102. [3] H.B. Stewart, Stability of two-phase flow calculation using two-fluid models, J. Comput. Phys. 33 (1979) 259–270. [4] V.H. Ransom, D.L. Hicks, Hyperbolic two-pressure models for two-phase flow, J. Comput. Phys. 53 (1984) 124–151. [5] J.R. Travis, F.H. Harlow, A.A. Amsden, Numerical calculation of two-phase flows, Nucl. Sci. Eng. 61 (1976) 1–10. [6] J.H. Stuhmiller, The influence of interfacial pressure forces on the character of two-phase flow model equations, Int. J. Multiphase Flow 3 (1977) 551–560. [7] J.C. Rousseau, R.L. Ferch, A note on two-phase separated flow models, Int. J. Multiphase Flow 5 (1979) 489–493. [8] K.H. Ardron, One-dimensional two-fluid equations for horizontal stratified two-phase flow, Int. J. Multiphase Flow 6 (1980) 295–304. [9] S.J. Lee, K.S. Chang, S.J. Kim, Surface tension effect in the two-fluid equation system, Int. J. Heat Mass Transfer 41 (1998) 2821–2826. [10] M.S. Chung, K.S. Chang, S.J. Lee, Wave propagation in two-phase flow based on a new hyperbolic two-fluid model, Numer. Heat Transfer A 38 (2000) 169–191. [11] M.S. Chung, S.J. Lee, K.S. Chang, Effect of interfacial pressure jump and virtual mass terms on sound wave propagation in the two-phase flow, J. Sound Vib. 244 (2001) 717–728. [12] P.A. Egelstaff, B. Widom, Liquid surface tension near the triple point, J. Chem. Phys. 53 (1970) 2667–2669. [13] R.D.J. Present, On the product of surface tension and compressibility of the liquids, J. Chem. Phys. 61 (1974) 4267–4269. [14] J.P. Aubin, I. Ekeland, Applied Nonlinear Analysis, Wiley, New York, 1984. [15] M. Blander, J.L. Katz, Bubble nucleation in liquids, AIChE J. 21 (1975) 833–848. [16] S.J.D. VanStralen, The mechanism of nucleate boiling in pure liquids and in binary mixtures – Part I, Int. J. Heat Mass Transfer 9 (1966) 995–1020. [17] H. Stadtke, G. Franchello, B. Worth, Numerical simulation of multi-dimensional two-phase flow based on hyperbolic flow equations, in: The 31st meeting of the European Two-Phase Flow Group, Piacenza, Italy, June 6– 8, 1994. [18] P.G. Buning, J.L. Steger, Solution of the two-dimensional Euler equations with generalized coordinate transforms using flux vector splitting, AIAA Paper 82-0971, 1982. [19] G.B. Whitham, Linear and Nonlinear Waves, Wiley, New York, 1974. [20] R.E. Henry, M.A. Grolmes, K.H. Fauske, Pressure pulse propagation in two-phase one- and two-component mixtures, ANL-7792, 1971. [21] A.R. Edwards, O’Brien, Studies of phenomena connected with the depressurization of water reactors, J. Br. Nucl. Energy Soc. 9 (1970) 125–135.
Numerical solution of hyperbolic two-fluid two-phase flow model with non-reflecting boundary conditions Moon-Sun Chung a
a,*
, Keun-Shik Chang b, Sung-Jae Lee
a
Thermal-Hydraulics Safety Research Team, Korea Atomic Energy Research Institute, 150 Dukjin-Dong, Yusong-Gu, Taejon 305-353, South Korea b Department of Aerospace Engineering, Korea Advanced Institute of Science and Technology, 373-1 Kusong-Dong, Yusong-Gu, Taejon 305-701, South Korea Received 21 July 2000; received in revised form 11 May 2001; accepted 9 July 2001 (Communicated by J.T. ODEN)
Abstract Flux vector splitting method is applied to the two-fluid six-equation model of two-phase flow, which takes account of surface tension effect via the interfacial pressure jump terms in the momentum equations. The latter terms using the concept of finite-thickness interface are derived as a simple function of fluid bulk moduli. We proved that the governing equation system is hyperbolic with real eigenvalues in the bubbly, slug, and annular flow regimes. The governing equations do not need any conventional artificial stabilizing terms like the virtual mass terms. Sonic speeds obtained through characteristic analysis show excellent agreement with the existing experimental data. Edwards pipe problem is solved as a benchmark test of the present two-phase equation model. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Two-fluid model; Surface tension; Hyperbolic system; Sonic speed; Flux vector splitting; Edwards pipe
1. Introduction The non-equilibrium exchange of mass, momentum, and energy between the gas and the liquid phases can be conveniently accounted for by using the two-fluid models. However, the common form of two-fluid models assuming single pressure across the interface gives complex eigenvalues
*
Corresponding author. Tel.: +82-42-868-2895. E-mail address: [email protected] (M.-S. Chung).
0020-7225/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 2 2 5 ( 0 1 ) 0 0 0 9 2 - 1
790
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
Nomenclature A interfacial area A; B coefficient matrices C speed of sound E source matrix ¼ A1 C F flux vector h enthalpy per unit mass J; K Jacobian matrices L fluid bulk modulus p pressure Q volumetric source term for heat R averaged bubble radius t time u internal energy per unit mass U primitive variable vector v flow velocity V control volume x space coordinate Greek symbols a volumetric phase concentration b relative surface thickness d surface thickness e eigenvalue switching factor / source term k eigenvalue K eigenvalue matrix q fluid density r surface tension Subscripts and superscripts g gas phase i interface j index for coordinate k index for each fluid l liquid phase m mth component n nth time-step, current time level p constant pressure s saturation condition w on the wall )1 inverse of a matrix first-order flux
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
791
and consequent numerical instability because the initial value problem is ill-posed, see [1–3]. In contrast, the models assuming two pressures across the interface present hyperbolic type governing equations with real characteristics. However, the existing two-pressure model represents only particular flow type like a stratified flow, see [4]. It has been known that mathematical property of the governing equations is improved by introducing physical terms in the governing equations. Ramshaw and Trapp [2], for example, added the surface tension equation to the governing equation system and Travis et al. [5] took account of the viscous stresses in the momentum equations. Stuhmiller [6] argued that when the interfacial pressure force terms were taken into account, the characteristic speed of the void wave became partially real. He analyzed that to assure real characteristics for the flow having unequal phase velocities, the interfacial pressure had to be lower than the bulk pressure by an amount proportional to the square of the relative velocity. Rousseau and Ferch [7], on the other hand, derived a condition for the hyperbolic governing equations in terms of static pressure difference between phases. They showed that the condition is identical to that of Kelvin–Helmholtz instability against the long wavelengths. Ardron [8] examined the one-dimensional two-fluid equations for a stratified flow, ignoring viscosity but retaining gravity and surface tension effect. The equations were observed stable over a realistic range of conditions, producing wave velocities in good agreement with the exact solution. A new promising approach of eliminating the instability was proposed by the present authors and their coworkers: The interfacial pressure jump terms are introduced based on the surface tension in the two-fluid momentum equations, see [9–11]. The system of equations manifested real characteristics in all the bubbly, slug, and annular flow regimes when the interfacial pressure jump terms are expressed as a product of effective bulk moduli and the gradient of interfacial area density. To avoid the difficulty of mathematically treating the interfacial area transport equation, we developed an expression of the interfacial pressure jump terms dependent only upon the fluid bulk moduli. In Section 2, we will recapitulate how the surface tension terms are brought about in the governing equations. Finite-volume flux vector splitting (FVS) method is formulated in order to solve the system of equations numerically in Section 3. Unlike the conventional methods of twophase flow, it is less dissipative and captures sharp gradient of flows very well without oscillation. In Section 4, we discuss the numerical solution of Edwards pipe problem as a benchmark test of transient two-phase flow.
2. Hyperbolic two-fluid model The conservation equation system consists of mass, momentum, and energy conservation equations by which the area-averaged phasic properties are related one another. The one-dimensional six equations take the following form: Continuity: oðak qk Þ oðak qk mk Þ þ ¼ /c;k ; ot ox
ð1Þ
792
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
Momentum: oðak qk mk Þ oðak qk m2k Þ opk oak þ ðpk pi Þ ¼ /c;k vi þ Fi;k þ Fw;k þ Bk ¼ /e;k ; þ þ ak ox ox ot ox
ð2Þ
Internal energy: oðak qk uk Þ oðak qk vk uk Þ oak oðak vk Þ þ þ pk ¼ Qi;k þ Fi;k vi þ /c;k hk þ Fw;k vk ¼ /e;k ; þ ot ox ox ot
ð3Þ
where ak ; qk ; pk ; vk ; uk , and hk denote void fraction, density, pressure, velocity, internal energy, and enthalpy of the phase k, respectively, and Bk is a gravitational body force. The source functions, /c;k ; /m;k , and /e;k , are made of algebraic constitutive relations, as shall be mentioned at the end of this section. We assign k ¼ g for the gas and k ¼ l for the liquid. The interfacial pressure jump term, ðpk pi Þðoak =oxÞ, which will be expressed as the product of fluid bulk modulus and an infinitesimal variable, is very small relative to other terms in the momentum equation (2). It shall be shown in Section 4 by the characteristic analysis that these terms make the equation system hyperbolic. 2.1. Interfacial pressure jump terms Young and Laplace proposed a well-known surface tension: pg pl ¼
2r ; R
ð4Þ
where r is the surface tension and R is the averaged bubble radius. We now assume a finite interfacial thickness d between the two radii Rg and Rl shown in Fig. 1. For an imaginary sphere having a radius Ri , where Rg þ Rl ¼ 2Ri , we rewrite Eq. (4) as
Fig. 1. Hypothetical sphere at Ri (solid circle) with an infinitesimal thickness d.
pg pl ¼
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
793
r 2d 2d r ¼ : Rg þ d=2 d Rl d=2 d
ð5Þ
For small dilation of the bubble marked by the radial increment DRi and surface area change DAi , we can expect corresponding volume change interior and exterior to the film, namely, DVg and DVl . From the bubble geometry, we can easily derive the following relations:
d=2 ; Rg þ d=2 Rl DAi d=2 ¼ 1þ : Rl d=2 2 DVl Rg DAi ¼ 2 DVg
1
ð6Þ ð7Þ
We can rewrite Eq. (5) as 4r Rg DAi 4r Rl DAi ¼ : 1 1þ pg pl ¼ d d 2 DVg 2 DVl
ð8Þ
It is recognized that the relation, L ¼ cr=d, is used in the physical chemistry and statistical mechanics, where L is a bulk modulus and c is a constant, see [12,13]. We assume that the surface tension stress, 4r=d ¼ L, which plays a role of Lagrangian multiplier as explained by Aubin and Ekeland [14], can be replaced by the phasic bulk moduli added together: 4r ¼ Lg þ Ll ¼ qg Cg2 þ ql Cl2 ; d
ð9Þ
where Ck is the sonic speed. The pressure jump is split into the phasic components as shown in [9–11]: Rg DAi Rl DAi pg pl ¼ ðpg pi Þ þ ðpi pl Þ ¼ Lg 1 Ll 1 þ : ð10Þ 2 DVg 2 DVl Here, Lg and Ll are the bulk modulus of the gas and the liquid phase, respectively, and pi is the interfacial pressure on the imaginary sphere Ri . The pressure jump interior and exterior to the sphere Ri has the following fraction as derived from Eqs. (8)–(10), pg pi Lg ¼ ; pg pl Lg þ Ll
pl pi Ll ¼ : p g p l Lg þ Ll
ð11Þ
Consequently, the interfacial pressure jump terms in the momentum equation (2) have the expression ðpk pi Þ
oak 2r Lk oak oak ¼ ð1Þn ¼ ð1Þn bLk ; Ri Lg þ Ll ox ox ox
ð12Þ
794
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
where the exponent n stands for the liquid if n ¼ 1 and for the gas if n ¼ 2. The parameter b named ‘relative surface thickness’ becomes by using Eq. (9), b
2r 1 d ¼ : Ri Lg þ Ll 2Ri
ð13Þ
Therefore, the interfacial pressure jump term is represented by the product of bulk modulus Lk and the relative surface thickness b which is a function of temperature and pressure. For a bubbly flow assumed of a perfect mixture, the mixture bulk modulus is Lm ¼ V
dp dp dp ¼V : ¼ V dV dVg þ dVl Vg dp=Lg þ Vl dp=Ll
ð14Þ
Since it holds that Lg Ll , Eq. (14) yields Lm Lg =ag , which is true in the void fraction range ag P Lg =Ll . We assume here that the mixture bulk modulus is equal to the bulk modulus of the gas by taking ag ¼ Oð1Þ, namely, Lm ¼ qg Cg2 . The resultant surface thickness evaluated by Eq. (9) becomes d ¼ 2r=Lg , which is of the order of magnitude 106 –107 m as shown in Fig. 2. It is equivalent to the initial bubble radius in homogeneous nucleation, see [15]. For the slug flow, the effective bulk modulus Lk can be obtained from a simplified physical model using the assumption of no elastic interaction between the phases. Because a wave traveling in one phase is not disturbed by the presence of the other phase in this model, the effective bulk modulus of one phase is same as the bulk modulus of its sole phase, namely, Lg ¼ qg Cg2 and Ll ¼ ql Cl2 . In the annular flow, the pressure wave in the gas phase is not transmitted to the liquid phase but is mostly either reflected back or changed into capillary waves on the liquid surface. We can then assume Lg ¼ qg Cg2 and Ll ¼ 0. It is suggested by Van Stralen [16] that the minimum radius of the bubble nuclei can be expressed in the nucleation process as a function of saturation temperature, latent heat, and
Fig. 2. Bulk modulus and surface thickness of the water–vapor mixture.
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
795
superheating temperature. It is suggested that the minimum radius of bubble nuclei be of the order of magnitude 106 –107 m for many liquids. Therefore we can deduce that the radius of curvature Ri has the same order of magnitude as the surface thickness d. Eq. (5) suggests that for Rg P 0 and Rl Rg ¼ d, we must have the inequality 0 < b 6 1. For the bubbly flow, since the surface thickness and the bubble diameter have the same order of magnitude in the limit, we are led to the approximation, b ¼ Oð1Þ. In contrast, for the annular flow with an interface having a large radius of curvature, b becomes relatively small. On this reason, we can make interpolation of b as a function of void fraction. For the convenience, we treat this parameter as a constant in the characteristic analysis. 2.2. Transformation of governing equations Using the definition Xk ðoqk =opk Þuk , Yk ðoqk =ouk Þpk , and the identity opg =ox ¼ opl =ox which is derived from Eq. (4) for the equilibrium states of the bubble, we can write out Eqs. (1)–(3) as follows: Continuity: oag opg oag opg ovg oug oug þ ag Xg þ qg vg þ ag Xg vg þ qg ag þ ag Yg þ vg ¼ /c;g ; qg ot ot ox ox ox ot ox oal opg oal opg ovl oul oul ql þ al Xl þ ql vl þ al Xl vl þ ql al þ al Yl þ vl ¼ /c;l : ot ot ox ox ox ot ox
ð15Þ ð16Þ
Momentum: ag qg al ql
ovg opg ovg oag þ ag þ ag qg vg þ bLg ¼ /m;g ; ot ox ox ox
ovl opg ovl oal þ al þ al ql vl bLl ¼ /m;l : ot ox ox ox
ð17Þ ð18Þ
Internal energy: ag qg al ql
oug oug oag oag ovg þ ag qg v g þ pg þ pg v g þ pg ag ¼ /e;g ; ot ox ot ox ox
oul oul oal oal ovl þ al ql vl þ pl þ pl vl þ pl al ¼ /e;l : ot ox ot ox ox
ð19Þ ð20Þ
It holds that ag þ al ¼ 1. If we use an auxiliary thermodynamic equation dqk ¼
oq opk
dpk þ
uk
oq ouk
duk ¼ Xk dpk þ Yk duk ; Pk
ð21Þ
796
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
then a matrix equation holds A
oU oU þB ¼ C; ot ox
ð22Þ
where U ¼ ½ ag
pg
vg
vl
ug
ul T ;
ð23Þ
T C ¼ /c;g /c;l /m;g /m;l /e;g /e;l ; 3 2 0 0 ag Yg 0 qg ag Xg 6 ql al Xl 0 0 0 al Yl 7 7 6 7 6 0 0 a q 0 0 0 g g 7; A¼6 6 0 0 0 7 0 0 al ql 7 6 4 pg 0 0 0 ag qg 0 5 0 0 0 0 al ql pl 2 qg vg ag vg Xg ag qg 0 ag Yg vg 6 ql vl al vl Xl 0 al ql 0 6 6 bLg a a q v 0 0 g g g g B¼6 6 bLl a 0 a q v 0 l l l l 6 4 pg vg 0 ag pg 0 ag qg vg 0 0 al pl 0 pl vl
ð24Þ
ð25Þ
3 0 al Yl vl 7 7 0 7 7: 0 7 7 0 5 al ql vl
ð26Þ
The source vector C in Eq. (22) requires constitutive equations related with the transport of mass, momentum, and energy between phases. To maintain the consistency with the interfacial transfer relations, we prefer the same constitutive equations as were adopted in RELAP5/MOD3 code for the bubbly, slug, and annular flow regimes. In practice, we identify flow transition regimes like the bubble-to-slug and the slug-to-annular flow regimes by interpolation as a function of void fraction. Eq. (22) can be transformed into oU oU þG ¼ E; ot ox
ð27Þ
where G ¼ A1 B and E ¼ A1 C. The governing equations are then straightforward: oV oV þH ¼ M; ot ox
ð28Þ
where H ¼ J G J 1 , M ¼ J E, and J ¼ oV=oU. The conservation vector V is V ¼ ag qg
al ql
ag qg v g
al ql vl
ag qg ug
al ql ul
T
:
ð29Þ
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
We define the flux vector F by 3 2 ag qg vg 7 6 al ql vl 7 6 2 6 ag qg vg þ pg ag þ bLg ag 7 7; 6 F¼6 2 7 6 al ql vl þ pl al þ bLl al 7 5 4 ag qg vg ug al ql vl ul
797
ð30Þ
where the relative surface thickness b is a relatively infinitesimal variable. An alternative form of Eq. (28) is oV oF oF þ þW ¼ M; ot ox ox
ð31Þ
where W ¼ ðJ G KÞ K 1 and K is another Jacobian matrix oF=oU. To obtain the flux vector F, Eq. (31) is further transformed to oF oF þR ¼ N; ot ox where R ¼ K G 2 qg 6 ql 6 6 qg vg J ¼6 6 ql vl 6 4 qg ug ql ul
ð32Þ
K 1 , N ¼ K E, and the Jacobian matrices J, K are 3 ag Xg 0 0 ag Yg 0 7 al Xl 0 0 0 al Yl 7 7 ag vg Xg ag qg 0 ag vg Yg 0 7; 7 al vl Xl 0 al ql 0 al vl Yl 7 5 ag ug Xg 0 0 ag ðug Yg þ qg Þ 0 al ul Xl 0 0 0 al ðul Yl þ ql Þ
2
qg vg 6 ql vl 6 6 qg v2g þ pg þ bLg K ¼6 6 q v2 pl bLl 6 l l 4 qg vg ug ql vl ul
ag vg Xg al vl Xl ag v2g Xg þ ag al v2l Xl þ al ag vg ug Xg al vl ul Xl
ag qg 0 2ag qg vg 0 ag qg ug 0
0 al ql 0 2al ql vl 0 al ql ul
ag v g Yg 0 ag v2g Yg 0 ag vg ðug Yg þ qg Þ 0
ð33Þ
3 0 7 al vl Yl 7 7 0 7: 2 7 al vl Yl 7 5 0 al vl ðul Yl þ ql Þ ð34Þ
3. Flux vector splitting The fluxes normal to the cell interfaces are calculated from the Riemann problem, which is shown in the work of Stadtke et al. [17]. For this purpose, the governing equation (27) is projected normal to the cell interface:
798
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
oU oU þ Gn ¼ En: ot on
ð35Þ
Corresponding to the six real eigenvalues, a linearly independent set of eigenvectors exists which is a necessary condition for a hyperbolic system of equations. A similarity transformation exists by which the characteristic form of the equations is obtained: L1 n
oU oU þ Kn L1 ¼ L1 n n En; ot on
ð36Þ
where the eigenvalue matrix is Kn ¼ L1 n G n Ln . The columns of matrix Ln are the right eigenvectors of the coefficient matrix G n and the rows of L1 n are the left eigenvectors of G n . Using the individual eigenvalues, the coefficient matrix G n can be split into the elementary parts Gn ¼
6 X
G n;m ;
ð37Þ
m¼1
with G n;m ¼ Ln Kn;m L1 n . The diagonal matrix Kn;m consists of the mth eigenvalue only. To obtain flux vectors at the cell interfaces, Eq. (32) is used as oF oF þ Rn ¼ N n; ot on
ð38Þ
where the coefficient matrix Rn ¼ K n G n K 1 n with K n ¼ oFn =oU. The similarity transformation does not change the eigenvalues of the governing equation system. Therefore, the characteristic form of the governing equations, equivalent to Eq. (36), is L01 n
oF oF þ Kn L01 ¼ L01 n n N n; ot on
ð39Þ
where the transformed right eigenvectors are L01 ¼ K n Ln . n If we linearize Eq. (39), then the numerical flux at the cell interface is composed of the positive and negative contributions depending on the sign of eigenvalues as X Rm X Rm Fj þ Fjþ1 ; ð40Þ Fjþ1=2 ¼ km j km jþ1 m;km >0 m;km <0 where km denotes the mth component of eigenvalue matrix. The resultant form of the governing equations based on the finite volume method is obtained from the conservative equation (31) as Vjnþ1 ¼ Vjn
Dt Dt Fjþ1=2 Fj1=2 Wjn Fjþ1=2 Fj1=2 þ DtMjn : Dx Dx
ð41Þ
A non-physical discontinuity might appear at the sonic transition of the computed results. It can be corrected by a simple modification similar to that used by Buning and Steger [18]. It redefines the eigenvalues k m using a small parameter em , namely,
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
k m ¼ km
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 þ ðem =km Þ2 2
799
ð42Þ
:
We assign a very small number to em , for example, of the order ðem =km Þ2 6 Oð104 Þ, to make modification of the eigenvalues truly small. 4. System eigenvalues and numerical example 4.1. Characteristic analysis The eigenvalues of the equation system represent the propagation speeds of small-amplitude, short-wavelength perturbations as indicated by Whitham [19]. Whereas both source and dispersion terms of the source vector play an important role for long-wavelength disturbances and the nonlinear wave interaction causes dominant effect on the large-amplitude disturbances. If the eigenvalues are all real, the equation system is hyperbolic and its solutions are stable against small disturbances. The eigenvalues of coefficient matrix G in Eq. (27) are obtained by DetðG kIÞ ¼ 0:
ð43Þ
Here, we assume that the ‘relative surface thickness’ is a constant, 1.0, in order to derive the eigenvalues analytically from Eq. (43). As a result, we can derive a sixth-order polynomial equation as follows: P6 ðkÞ ¼ ðk vg Þðk vl Þðk4 þ Z1 k3 þ Z2 k2 þ Z3 k þ Z4 Þ ¼ 0;
ð44Þ
where the coefficients are Z1 ¼ 2ðvg þ vl Þ, Z2 ¼
Z3 ¼
1 þ ag Cl2 ql # ( " ) L L 2 2 g l al Cg2 qg vg þ vl þ 2vg vl Cl2 þ ag Cl2 ql vg þ vl þ 2vg vl Cg2 ; qg ql al Cg2 qg
2 þ ag Cl2 ql !# ( " ) L L g l v2g þ ag Cl2 ql vl Cg2 v2g þ vg al Cg2 qg vg Cl2 v2l þ vl v2l ; qg ql al Cg2 qg
1 Z4 ¼ 2 al Cg qg þ ag Cl2 ql
( al Cg2 qg
"
Lg v2g Cl2 v2l qg
!# þ ag Cl2 ql
Cg2 v2g
L
l
ql
) v2l :
800
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
Table 1 Three sets of six eigenvalues Eigenvalues Bubbly regime
Slug regime Annular regime
k1;2 ¼ vg ; vl ;
k3;4 ¼ vg Cg ;
k5;6 ¼ vl Cl
k1;2 ¼ vg ; vl ;
k3;4 ¼ vg Cg ;
k5;6 ¼ vl Cl
k1;2 ¼ vg ; vl ;
k3;4 ¼ vg Cg ;
k5;6 ¼ vl Cl
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qg Cg2 2 al Cg qg þ ag Cl2 ql
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi al qg Cg2 al Cg2 qg þ ag Cl2 ql
The analytical solution of the characteristic equation (44) gives three sets of six real eigenvalues as listed in Table 1, which represent the bubbly, slug, and annular flows. The first two eigenvalues, k1;2 ¼ vg ; vl , represent the convection velocity of the gas and the liquid phases, respectively. The other two eigenvalues, k3 and k5 , represent approximately the sonic speeds in the gas and the liquid phases. For the bubbly and slug flows, the total sonic speed can be obtained by the void fraction weighting [10,11]: Ct ¼
k3 k5 : al k3 þ ag k5
ð45Þ
For the annular flow, the individual phasic sonic speeds, k3 and k5 , shall be used. In Fig. 3, we compared the computed total and phasic sonic speeds with the experimental data produced by Henry et al. [20]. For the bubbly flow, the total sonic speed agrees reasonably well with the experimental data in the void fraction range 0 < ag < 0:2 as shown in Fig. 3(a). The increasing deviation in the range ag > 0:2 is probably caused by transition of the flow regimes. Fig. 3(b) shows that the sonic speed of the water–air slug flow is in good agreement between the computed and experimental data in the entire void fraction range 0 < ag < 1. The sonic speed of the gas phase of the annular flow agrees well in Fig. 3(c) with the experimental data [20]. For the liquid phase, unfortunately, experimental data do not exist. The computed result shows that sonic speed of the liquid phase is subject to a rapid initial decrease for the very low the void fraction range, due to the effect of increasing elasticity at the interface. 4.2. Edwards pipe problem This problem is the benchmark test for Edwards and O’Brien’s [21] blowdown experiment. It represents a loss-of-coolant accident (LOCA) problem in the pressurized water reactor. The pipe is horizontal, 4-m long with the cross-section area 0:00456 m2 . It is filled with subcooled water under the initial pressure 7.0 MPa. Discharge of the flow is abruptly initiated by suddenly opening one end of the pipe that has a narrow cross-section area 0:00397 m2 . An expansion wave sweeps backward from the open end of the pipe, accompanied by flashing of the subcooled water due to the severe depressurization.
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
801
Fig. 3. Calculated sonic speeds compared with the experimental data: (a) Bubbly flow (pg ¼ 283 kPa). (b) Slug flow (pg ¼ 172 kPa). (c) Annular flow (pg ¼ 100 kPa).
In Fig. 4(a), the pressure is predicted at the open end of the pipe; the pipe is composed of 20 cells. The calculated time-dependent vapor pressure lies close to the measured pressure data. After the expansion wave is reflected at the closed end of the pipe, which is observed by the slight undershoot of the pressure in Fig. 4(a), the pressure is maintained near 2.8 MPa for about 0.15 s. The result of the present model shows that the pressure undershooting is less and the second depressurization occurs somewhat earlier in the calculated data than in the experimental data. However, the overall comparison appears as good as RELAP5/MOD3 code does. Fig. 4(b) shows the calculated time-dependent void fraction at the midsection of the pipe. Significant discrepancy which might have been originated from the inaccurate constitutive equations is apparent in the initial stage of the flow development, up to t ¼ 0:10 s. This deviation of the void fraction is also evident with the RELAP5/MOD3 code. Nevertheless, comparison of the present results with the experiment is significantly improved in later times. In contrast,
802
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
(a)
(b)
Fig. 4. Edwards pipe blowdown problem: (a) pressure history at the open end of the pipe; (b) void-fraction history at the midsection of the Edwards pipe.
RELAP5/MOD3 shows irregular void-fraction change up to t ¼ 0:25 s. The present result is also smoother than the earlier data of the present authors [10]. 5. Conclusions A hyperbolic system of two-phase two-fluid conservation laws has been derived using the interfacial pressure jump terms expressed ultimately as a function of fluid bulk moduli. Owing to its hyperbolicity, the equation system can be solved by an upwind numerical method like the flux vector splitting method. All the eigenvalues of the equation system turned out to be real which must be a remarkable improvement over the earlier two-phase formulations. The computed sonic speeds gave good comparison with the measured data. It has been demonstrated that the present model can be efficiently applied to the wave-dominant two-phase flows having strong expansion and shock waves with abrupt depressurization and flashing phenomena such as Edwards pipe problem.
Acknowledgements The authors are grateful to Dr. Won-Jae Lee in KAERI for the helpful technical discussions he offered while we prepare this paper. References [1] R.W. Lyczkowski, D. Gidaspow, C.W. Solbrig, E.D. Hughes, Characteristics and stability analyses of transient one-dimensional two-phase flow equations and their finite difference approximations, Nucl. Sci. Eng. 66 (1978) 378–396.
M.-S. Chung et al. / International Journal of Engineering Science 40 (2002) 789–803
803
[2] J.D. Ramshaw, J.A. Trapp, Characteristics, stability, and short-wave length phenomena in two-phase flow equation systems, Nucl. Sci. Eng. 66 (1978) 93–102. [3] H.B. Stewart, Stability of two-phase flow calculation using two-fluid models, J. Comput. Phys. 33 (1979) 259–270. [4] V.H. Ransom, D.L. Hicks, Hyperbolic two-pressure models for two-phase flow, J. Comput. Phys. 53 (1984) 124–151. [5] J.R. Travis, F.H. Harlow, A.A. Amsden, Numerical calculation of two-phase flows, Nucl. Sci. Eng. 61 (1976) 1–10. [6] J.H. Stuhmiller, The influence of interfacial pressure forces on the character of two-phase flow model equations, Int. J. Multiphase Flow 3 (1977) 551–560. [7] J.C. Rousseau, R.L. Ferch, A note on two-phase separated flow models, Int. J. Multiphase Flow 5 (1979) 489–493. [8] K.H. Ardron, One-dimensional two-fluid equations for horizontal stratified two-phase flow, Int. J. Multiphase Flow 6 (1980) 295–304. [9] S.J. Lee, K.S. Chang, S.J. Kim, Surface tension effect in the two-fluid equation system, Int. J. Heat Mass Transfer 41 (1998) 2821–2826. [10] M.S. Chung, K.S. Chang, S.J. Lee, Wave propagation in two-phase flow based on a new hyperbolic two-fluid model, Numer. Heat Transfer A 38 (2000) 169–191. [11] M.S. Chung, S.J. Lee, K.S. Chang, Effect of interfacial pressure jump and virtual mass terms on sound wave propagation in the two-phase flow, J. Sound Vib. 244 (2001) 717–728. [12] P.A. Egelstaff, B. Widom, Liquid surface tension near the triple point, J. Chem. Phys. 53 (1970) 2667–2669. [13] R.D.J. Present, On the product of surface tension and compressibility of the liquids, J. Chem. Phys. 61 (1974) 4267–4269. [14] J.P. Aubin, I. Ekeland, Applied Nonlinear Analysis, Wiley, New York, 1984. [15] M. Blander, J.L. Katz, Bubble nucleation in liquids, AIChE J. 21 (1975) 833–848. [16] S.J.D. VanStralen, The mechanism of nucleate boiling in pure liquids and in binary mixtures – Part I, Int. J. Heat Mass Transfer 9 (1966) 995–1020. [17] H. Stadtke, G. Franchello, B. Worth, Numerical simulation of multi-dimensional two-phase flow based on hyperbolic flow equations, in: The 31st meeting of the European Two-Phase Flow Group, Piacenza, Italy, June 6– 8, 1994. [18] P.G. Buning, J.L. Steger, Solution of the two-dimensional Euler equations with generalized coordinate transforms using flux vector splitting, AIAA Paper 82-0971, 1982. [19] G.B. Whitham, Linear and Nonlinear Waves, Wiley, New York, 1974. [20] R.E. Henry, M.A. Grolmes, K.H. Fauske, Pressure pulse propagation in two-phase one- and two-component mixtures, ANL-7792, 1971. [21] A.R. Edwards, O’Brien, Studies of phenomena connected with the depressurization of water reactors, J. Br. Nucl. Energy Soc. 9 (1970) 125–135.