- Email: [email protected]
Analysis of density wave instability in a boiling flow using a characteristic method
Nonlinear
Analysis,
Theory,
Mehods
Pergamon
& Application, Vol. 30, No. 5. pp. 2707-2795. 1997 Proc. 2nd World Congress of Nonlinear Analysrs 0 1997 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0362-546x197 $17.00 + 0.00
PII: so362-546X(97)003&-4
Analysis of Density Wave Instability in a Boiling Flow Using a Characteristic Method MOTOAKI
OKAZAKI
Department of Reactor Engineedng, Japan Atomic Energy Research Institute, Tokai-mura, Naka-gun, Ibamki-ken 319-11 Key words
and phases:
Density
wave instability,
two-phase
flow, characteristic
Japan method,
sound speed, numerical
method
1. Introduction Density wave instability in a boiling flow is known to occur in a heated channel with a constant pressure difference between inlet and outlet. ln order to make analysis of the instability, physical aspects of unsteady boiling flow must be described in an essentially correct manner. The characteristic method for derivation of ordinary differential equations from partial differential equations of two-phase flow is very appropriate to meet this objective. Because, this method is only one which describes the propagation of unsteady flow change in a channel in conformity with its physical meaning. However this method is not fully developed for the practical use although some paper@‘@) were presented until now. We presented the numerical method (4) of analyzing unsteady steam or water single flow and mixing two-phase flow in a heated channel, separately. In this paper, we present the analyzing method of combined flow of subcooled water and boiling flow in a heated channel and some calculated results of density wave instability.
2. Basic Equations The basic equations Mass conservation Momentum balance Energy conservation
for homogeneous %
:
saturated two-phase
+g(p*),=
flow are as follows.
0
(2.1)
au +,au+Lap+Lap,=o
: :
at
az
P, h
P, az
3. Derivation of Characteristic Equations (1) Characteristic
Equation Change the mass conservation
of State Change equation (2.1) as follows.
2787
(2.2)
2788
Second World
Congress
of Nonlinear
Analysts
(3.1)
ChangeEq.(2.3), using the following relations p,=aP,+(l-ah,
,
h =xh,+(l-x)h, m
) c=h-$
, x=Qp, , w=p#i hn
(34
and combining Eqs.(Zl) and (2.2), we get (3.3)
We can consider that the left hand side of Eq.(3.3) expressesthe entropy change by thermodynamicalstate changealong the flow path and therefore the right hand side of Eq(3.2) does the heat contributing to entropy increase. Consequently Eq(3.3) can be transformed into the characteristicequation(3.4) as follows. along
Where
(3.5)
Equation (3.4) is the ordinary differential equation describingthe thermodynamicalstate change occurredin a flow. For the saturatedtwo-phase flow where phasechangeoccurs receiving external heat, the equationof phasechangeis derived from Eq.(3.4). By using the following relations x2%
1-x=
.
(I-alp,
pill
Pm ’
for the Eq.(3.4), we get the phasechangeequation as follows. (3.6) (3.7)
(2) Characteristic equation of momentum balance (2.1) Water single flow If we assumeas P =f(pJ)
, then
% =(2J,+(%),$
(3.8)
From the general expressionsof thermodynamics,we get
(~),=~(~),5 (ig=: ’ ($go)p($);-l If we use
(3.9)
and substitutingEqs.(3.9) for Eq.(3.8), then Eq(3.8) becomesas follows.
Second
World
Congress
of Nonlinear
Analysts
2789
(3.8)’ Substituting Eq.(3.8)’ into JZq.(3.1), we get
Further, Eq.(3.10) is changed to Eq.(3.11) using Eqs.(3.5). (3.11) Now, Eqs.(3.11) and (2.2) are considered as partial differential equations for to
t
and
z , respectively.
P and u with respect
From these two equations, we can get the characteristic equations as
follows. dz --=U+a, dt
along
dz -=u-a
along
1dP du --+p--&=---&
DS -, ~---$
a&
,
dr
PT
1dP ---+p-=-ad2
du a?
PT aC”K
DS Dt-
dp1 a2
(3.13)
Where (2.2)
Saturated two-phase
flow
As for the sonic velocity in the two-phase
flow, we consider dx=O in addition to &=O
, in
the pressure propagation process, because phase change does not occur in the very fast transient. lf we assume P, =.RpA Let V.,,=l/P,
To consider
, then
(G-i$={(z$+[2q,~}~
(3.14)
and using Eqs.(3.2), for the first term of the right hand side of E!q.(3.14), we get
h=O
in each phase in the pressure propagation process, we use the relation for gas or
liquid.
Change the above equation, using the general expressions of thermodynamics, we get
SecondWorld Congress of NonlinearAnalysts
2790
(3.16)
Substituting Eq.(3.16) into Eq.(3.1.5), we can get as follows. (3.17)
Where,
(3.18)
(3.19)
The second term of the right hand side of Eq(3.14)
is expressed as follows. (321)
Substituting Eqs.(3.14) and (3.17) into Eq.(3.1), we get
The secondterm in the Eq(3.22) expressesconcerningthe entropy changeof eachphaseand the third doesconcerning the quality changealong the flow path, respectively. Therefore, Eq(3.22) becomesas follows, using Eq.(3.21).
Now, Eqs.(3.23)and (2.2) are consideredas partial differential equationsfor P and u with respect to t and z , respectively. From these two equations,we can get the characteristic equationsas follows. along
along
dz -=u+a,,
1 dP du 0, DP --+pm-=--+a,L---2Dt ax dt a? ah
dt
1dP
i&-, dt
*’
-a,-;i;+%-
du -----aa, DP a; Dt
V -VI
2
Dx Dt
aP,
az
(3.24)
vs-v~ .Dx -3 = ,,; Dt &
3. Solution method by using difference equations Characteristicsgrid method is used for the numerical calculation along the each characteristic curve as an initial value problem. In Fig.3.1 and Fig.3.2, solid lines show the characteristiccurvesfor
Second World Congress of Nonlinear Analysts
2791
the changes of pressureand velocity in the flow causedby the pressurepropagation,
and dotted lines
showthe onesfor the changeof enthalpy or quality in the flow causedby the fluid flow. Initial value of enthalpy or quality for the solid lines are given by the solutionsobtainedon the characteristiccurve of dotted lines and the initial value of pressureand velocity for the dotted lines are given by the solutionsobtained on the characteristiccurve of solid lines, for the next step calculation, respectively. When characteristicscurves deal only liquid single flow or only two-phase flow, numerical calculationscan be made in a very stable manner. But, in the case that boiling initiation occurs between the two characteristicscurve of RQ and SQ in Fig.2, numerical calculation becomesvery difficult to obtain solutionsin a stable mannerdue to the vapor generationand the suddendecreaseof sonic velocity. The sudden change of sonic velocity caused by boiling initiation makes the characteristicscurve steepand the time interval larger, and it makeslarger amountof boiling and wider area of boiling causinga larger changeof density. All thesemattersare very sensitiveto the changesof pressureand velocity, which makesdifficult to obtain an exact solution. Further, a remarkabledistortion of characteristicscurve are made due to suddendecreaseof sonic velocity, which meansa large difference of time levels betweenadjacentcurve is calculated. If it is small, adjacenttime level approachesto about the samelevel in making progress of calculationby being made smallerdistancebetweenpoints R and S in Fig.3.2. But if it is greater than a certain limited value, the difference of time levels cannot be made smaller any more. According as the distancebetweenR and S becomessmall, the pitch of characteristicline of cizldt=U (we call this u-line) is neededto be smaller, becauseby using these lines the average values of density and quality are estimatedin the region between R and S. For this purpose, the numberof u-line are increasedby using the interpolation methodjust before boiling initiation point (we call boiling front) comesto the region due to the increaseof heat input. The numerical calculation becomesmore accurate with the increaseof number of u-line. The spatial average value of fluid condition is determinedby the numberof u-line and the averagevalue of time is determinedfor each u-line betweenpoint 1 to 2 in Fig.3.2. The averagevalues of time and spaceare evaluated for the fluid of each region of RQ and SQ separately. In the region of R-S, for the u-line with no evaporation, the velocity u and density p are determinedas the values of upper streampoint R, and for the u-line with evaporation, quality change is determinedby using the velocity u which is determinedinterpolatedly betweenthe boiling initiating point and the downstreampoint S. Using the averagevalues of time and spacefor each of regions RQ and SQ+ respectively, the solutions for point Q are obtained. And if the iterative solutions converged,we take them as the final solutions.
‘Figs. 1
Characteristics
grid
. Fig*3*2
Characteristics on 2-t plane
lines
Fig.4.1
Test pipe
Second World
2792
Congress
of Nonlinear
Analysts
4. Some examples of numerical calculation results The two cases in the boiling flow from subcooled water to saturated two-phase flow are presented. One is the steady boiling flow and the another is an oscillatory flow due to density wave instability which is caused by a little larger heating rate than in the case of steady boiling flow. Unsteady flow calculations are made for the test pipe shown in Fig.4.1, where both ends are opened and at one end pressure is kept constant 3.5 MPa and at the another end pressure is decreased gradually until 2.0 MPa. For the case of steady flow, boundary conditions of pressure and heating rate given for test pipe are shown in Fig.4.2 and Fig.4.3, respectively. The inlet temperature of water is 215’C, which means 27.5-C subcooling to the inlet pressure 3.5 MPa. As for the characteristics grid, the number of calculation points is 41 including both ends, which means the test pipe is divided into N=40, and the number of u-lines is set as NxS initially. Calculation results for the case of steady boiling flow at grid No.=30 are shown in Figs.4.4-4.8. The steady boiling flow is attained after the reach of constant heat input. We can observe it in each of pressure and velocity staying about a constant value although a very small high frequent It is not clear that this oscillation means the oscillation is contained, as shown in Figs.4.6-4.7. phenomena of boiling flow itself or a kind of numerical instability occurred in the calculation of boiling initiation region. Anyway, boiling front stays at about 3.4 m as shown in Fig.4.8. The change of void fraction and acoustic velocity are shown in Fig.4.6 and Fig.4.7, respectively. Calculation results for the case of density wave instability in a boiling flow from a subcooled water are shown in Figs.4.9-4.16. First, the heat flux change is shown in Fig.4.9. The move of characteristics grid point caused by the sonic velocity change due to boiling is shown in Fig.4.lO(a),@). The density wave instability behaviors calculated in the subcooled water region (grid No.=4) are shown in Fig.4.11 and Fig.4.12, for pressure and velocity, respectively. And the behaviors in the saturated boiling region (grid No.=30) are shown in Figs.4.13-4.15, for pressure, velocity, void fraction and move of boiling front, respectively. From these results, the followings are shown. The oscillation phases of density wave instability in the subcooled region and boiling region are almost in the same phase, except that the phase of velocity in boiling region show a small phase lag. The phase of pressure shows a 180 degree difference compared to that of velocity. The pressure profile in the flow channel shows a monotonous pressure decrease because pressure oscillation is in phase through the flow channel.
0.0
Fig.4.2
Pressure
change
at exit
Fig.4.3
3.0
External heat case of steady
6.0
flux change flow
12.0
9.0
in the
15.c
Second World
Congress
of Nonlinear
Analysts
2793
3.2
2.4
Pressure
cig.4.4
change at Grid No.=30
Fig.45 1600~’ 1200
Velocity
change at Grid No.=30
’ ’ ’ I ’ ’ ’ ’ I ’ ’ ’
I
I4--
I ’ ’ ’ ’ ! ’ ’ ’
1
.-.--I ’_....__.___ _.______ -j_____- ._______ i ______.____ 1 ; / 1
..-.-. ._..i-.... --.---i-.-_.__ j...__-._._ .i___..._.__ kc i”=(:i _ .._..__ 1 4! 1.-.-. - .i...-.-.-.--.i-.____. -.+..-~ .___ -i--‘0
Fig.45
‘Void fraction change
at Grid No.=30
orI I I '1 0.0 3.0 6.0 9.0 Fig.4.7 Acoustic velocity change
12.0
15.0
at Grid No.=30
5 4 3 2 1 0 0.0
6.0
3.0
Fig.4.8 M&e 5_,11~,II~~,I~I~,~~~~_ LT.-
.-._._ j-.---
-.-.--
15.0
Fig.4.9 External
of Boiling front
.-.-. J-. I
case s 0
I
I 4
12.0
9.0
.-.- -.j”
--.-
-..-1
.-.--
Move of grid points
(No.=5,6,7,8)
I I +,q-e * L-L tx.-~c,~ .+-a---v-
,I ---+-
u- -w-e_ t’+)’ t
.-_
4
0
Fig.e.lO(a)
5
heat flux change in ,the of density wave instability flow
.-s 4
1 0 0I.0
FIgA.lO$r
12.0
16.0
fvloves~f grid points (No.=1 3,20,30,38)
Second World
2794
Congress
of Nonlinear
Analysts
12.0 16.0 a.0 Velocity change at Grid No.=4. (Subcooled water region) 4.0
0.0
Pressure change at Grid No.=4. (Subcooled water region)
Fig.4.11
Fig.4.12 40
I!!.I,III1
I"',"" I
I I
I
5z 301____. I I 1 --.-.__ .___ i---..--...-.c-.-.----.i._-. 1 20 .--.-.-.-.-j ,o
.-
_....
__.. _._._
--
__..___
---+
___._...__
r
.__.
-._
__._
____-_,
__!
_
.-
-.
.-....
-.-.-.-.-.
Io~o
Fig.4.13
16.0' 12.0 a.0 Pressure change at Grid No.=30. (Boiling region)
0.0
;
4.0
a.0
Velocity change (Boiling region)
Fig.4.14
at Grid N0.=30.
II ! ,/’-J--l 0.5-....-..‘.-.-...., j+L-._..-.. +. ._.._._______ c_____________: ;:I !! I 0.25 -----.c-.i---.------~.~.-.--.----~.~-....-.~.~-~ c ,/ 0 - I ,, d 0.0
Fig.4.15
j I 4 , 4.0
Time (s) c , I ,z,
a.0
Void fraction change (Boiling region)
1 ,* 12.0
16.0
at Grid No=30.
0.0
4.0
Fig.4.16
8.0
12.0
Move of boiling front
16.0
SecondWorld Congress of NonlinearAnalysts
2795
5. Nomenclature a : Sonic velocity m/s r : Latent heat J&z a, : Sonic velocity for saturated two-phase flow s : Specific entropy Jikg-K A : Flow area m2 t : Time set C, : $/& used in Fig.5. T : Temperature K C, : Specific isobaric heat capacity u : Fluid flow velocity m/s J/kg-K Cy : Specific isochoric heat capacity J/kg-K V : Specific volume m3/kg d : Diameter m W : Mass flow rate kg/s e : Specific internal energy Jk X : Flow quality h : Specific enthalpy J&z z : Distance m H : Defined by Eq(3.7) a : Void fraction e : Length m p : Coefficient of thermal expansion l/K Lp : Peripheral length kg/m3 p : Density N : Number of division of test pipe 4 : Heat flux w/m2 P : Pressure Pa K : isothermal compressibility l/Pa ¶ : Heat contributing to entropy increase of unit mass Jh (subscript) E : External heat Q : Point Q f : Friction R : Point R g : Gas S : Point S Q : Liquid TP:Two-Phase mixture m : g, !Z or mixture 6. Conclusion Numerical calculation method using characteristics grid for analyzing unsteady boling flow from subcooled water is presented. (2) Caluculation results of two types of boiling flow are presented. One type is steady boiling flow which is attained after heat flux becomes constant. Another type is unsteady boiling flow with density wave instability even after constant heat flux is given. (1)
References (1) (2) (3)
(4)
Banejee, S. and Hancox, W.T., “On the Development of Methods for Analysing Transient Eh,+ Boiling,” Int. J. Multiphase Flow., Vo1.4, ~~437-460 (1978) Gidaspow, D., Rasouli, F. and Shin, Y.W., “An Unequal Velocity Model for Transient Two-Phase Flow by the Method of Characteristics,” Nucl. Sci. and Engg., 84, pp179-195 (1983) Lyczkowski, R.W., Gidaspow, D., Solbrig, C.W. and Hughes, E.D, “Characteristics and Stability Analysis of Transient One-Dimensional Two-Phase Flow Equations and Their Finite Difference Approximations,” Nucl. Sci. and Engg., 66, pp378-396 (1978) Okazaki, M., “Analysis of Density Wave Instability in a Boiling Flow Using a Characteristic Method “ASME Conference. PVP-Vo1.270, 65-73 (1994)
Analysis,
Theory,
Mehods
Pergamon
& Application, Vol. 30, No. 5. pp. 2707-2795. 1997 Proc. 2nd World Congress of Nonlinear Analysrs 0 1997 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0362-546x197 $17.00 + 0.00
PII: so362-546X(97)003&-4
Analysis of Density Wave Instability in a Boiling Flow Using a Characteristic Method MOTOAKI
OKAZAKI
Department of Reactor Engineedng, Japan Atomic Energy Research Institute, Tokai-mura, Naka-gun, Ibamki-ken 319-11 Key words
and phases:
Density
wave instability,
two-phase
flow, characteristic
Japan method,
sound speed, numerical
method
1. Introduction Density wave instability in a boiling flow is known to occur in a heated channel with a constant pressure difference between inlet and outlet. ln order to make analysis of the instability, physical aspects of unsteady boiling flow must be described in an essentially correct manner. The characteristic method for derivation of ordinary differential equations from partial differential equations of two-phase flow is very appropriate to meet this objective. Because, this method is only one which describes the propagation of unsteady flow change in a channel in conformity with its physical meaning. However this method is not fully developed for the practical use although some paper@‘@) were presented until now. We presented the numerical method (4) of analyzing unsteady steam or water single flow and mixing two-phase flow in a heated channel, separately. In this paper, we present the analyzing method of combined flow of subcooled water and boiling flow in a heated channel and some calculated results of density wave instability.
2. Basic Equations The basic equations Mass conservation Momentum balance Energy conservation
for homogeneous %
:
saturated two-phase
+g(p*),=
flow are as follows.
0
(2.1)
au +,au+Lap+Lap,=o
: :
at
az
P, h
P, az
3. Derivation of Characteristic Equations (1) Characteristic
Equation Change the mass conservation
of State Change equation (2.1) as follows.
2787
(2.2)
2788
Second World
Congress
of Nonlinear
Analysts
(3.1)
ChangeEq.(2.3), using the following relations p,=aP,+(l-ah,
,
h =xh,+(l-x)h, m
) c=h-$
, x=Qp, , w=p#i hn
(34
and combining Eqs.(Zl) and (2.2), we get (3.3)
We can consider that the left hand side of Eq.(3.3) expressesthe entropy change by thermodynamicalstate changealong the flow path and therefore the right hand side of Eq(3.2) does the heat contributing to entropy increase. Consequently Eq(3.3) can be transformed into the characteristicequation(3.4) as follows. along
Where
(3.5)
Equation (3.4) is the ordinary differential equation describingthe thermodynamicalstate change occurredin a flow. For the saturatedtwo-phase flow where phasechangeoccurs receiving external heat, the equationof phasechangeis derived from Eq.(3.4). By using the following relations x2%
1-x=
.
(I-alp,
pill
Pm ’
for the Eq.(3.4), we get the phasechangeequation as follows. (3.6) (3.7)
(2) Characteristic equation of momentum balance (2.1) Water single flow If we assumeas P =f(pJ)
, then
% =(2J,+(%),$
(3.8)
From the general expressionsof thermodynamics,we get
(~),=~(~),5 (ig=: ’ ($go)p($);-l If we use
(3.9)
and substitutingEqs.(3.9) for Eq.(3.8), then Eq(3.8) becomesas follows.
Second
World
Congress
of Nonlinear
Analysts
2789
(3.8)’ Substituting Eq.(3.8)’ into JZq.(3.1), we get
Further, Eq.(3.10) is changed to Eq.(3.11) using Eqs.(3.5). (3.11) Now, Eqs.(3.11) and (2.2) are considered as partial differential equations for to
t
and
z , respectively.
P and u with respect
From these two equations, we can get the characteristic equations as
follows. dz --=U+a, dt
along
dz -=u-a
along
1dP du --+p--&=---&
DS -, ~---$
a&
,
dr
PT
1dP ---+p-=-ad2
du a?
PT aC”K
DS Dt-
dp1 a2
(3.13)
Where (2.2)
Saturated two-phase
flow
As for the sonic velocity in the two-phase
flow, we consider dx=O in addition to &=O
, in
the pressure propagation process, because phase change does not occur in the very fast transient. lf we assume P, =.RpA Let V.,,=l/P,
To consider
, then
(G-i$={(z$+[2q,~}~
(3.14)
and using Eqs.(3.2), for the first term of the right hand side of E!q.(3.14), we get
h=O
in each phase in the pressure propagation process, we use the relation for gas or
liquid.
Change the above equation, using the general expressions of thermodynamics, we get
SecondWorld Congress of NonlinearAnalysts
2790
(3.16)
Substituting Eq.(3.16) into Eq.(3.1.5), we can get as follows. (3.17)
Where,
(3.18)
(3.19)
The second term of the right hand side of Eq(3.14)
is expressed as follows. (321)
Substituting Eqs.(3.14) and (3.17) into Eq.(3.1), we get
The secondterm in the Eq(3.22) expressesconcerningthe entropy changeof eachphaseand the third doesconcerning the quality changealong the flow path, respectively. Therefore, Eq(3.22) becomesas follows, using Eq.(3.21).
Now, Eqs.(3.23)and (2.2) are consideredas partial differential equationsfor P and u with respect to t and z , respectively. From these two equations,we can get the characteristic equationsas follows. along
along
dz -=u+a,,
1 dP du 0, DP --+pm-=--+a,L---2Dt ax dt a? ah
dt
1dP
i&-, dt
*’
-a,-;i;+%-
du -----aa, DP a; Dt
V -VI
2
Dx Dt
aP,
az
(3.24)
vs-v~ .Dx -3 = ,,; Dt &
3. Solution method by using difference equations Characteristicsgrid method is used for the numerical calculation along the each characteristic curve as an initial value problem. In Fig.3.1 and Fig.3.2, solid lines show the characteristiccurvesfor
Second World Congress of Nonlinear Analysts
2791
the changes of pressureand velocity in the flow causedby the pressurepropagation,
and dotted lines
showthe onesfor the changeof enthalpy or quality in the flow causedby the fluid flow. Initial value of enthalpy or quality for the solid lines are given by the solutionsobtainedon the characteristiccurve of dotted lines and the initial value of pressureand velocity for the dotted lines are given by the solutionsobtained on the characteristiccurve of solid lines, for the next step calculation, respectively. When characteristicscurves deal only liquid single flow or only two-phase flow, numerical calculationscan be made in a very stable manner. But, in the case that boiling initiation occurs between the two characteristicscurve of RQ and SQ in Fig.2, numerical calculation becomesvery difficult to obtain solutionsin a stable mannerdue to the vapor generationand the suddendecreaseof sonic velocity. The sudden change of sonic velocity caused by boiling initiation makes the characteristicscurve steepand the time interval larger, and it makeslarger amountof boiling and wider area of boiling causinga larger changeof density. All thesemattersare very sensitiveto the changesof pressureand velocity, which makesdifficult to obtain an exact solution. Further, a remarkabledistortion of characteristicscurve are made due to suddendecreaseof sonic velocity, which meansa large difference of time levels betweenadjacentcurve is calculated. If it is small, adjacenttime level approachesto about the samelevel in making progress of calculationby being made smallerdistancebetweenpoints R and S in Fig.3.2. But if it is greater than a certain limited value, the difference of time levels cannot be made smaller any more. According as the distancebetweenR and S becomessmall, the pitch of characteristicline of cizldt=U (we call this u-line) is neededto be smaller, becauseby using these lines the average values of density and quality are estimatedin the region between R and S. For this purpose, the numberof u-line are increasedby using the interpolation methodjust before boiling initiation point (we call boiling front) comesto the region due to the increaseof heat input. The numerical calculation becomesmore accurate with the increaseof number of u-line. The spatial average value of fluid condition is determinedby the numberof u-line and the averagevalue of time is determinedfor each u-line betweenpoint 1 to 2 in Fig.3.2. The averagevalues of time and spaceare evaluated for the fluid of each region of RQ and SQ separately. In the region of R-S, for the u-line with no evaporation, the velocity u and density p are determinedas the values of upper streampoint R, and for the u-line with evaporation, quality change is determinedby using the velocity u which is determinedinterpolatedly betweenthe boiling initiating point and the downstreampoint S. Using the averagevalues of time and spacefor each of regions RQ and SQ+ respectively, the solutions for point Q are obtained. And if the iterative solutions converged,we take them as the final solutions.
‘Figs. 1
Characteristics
grid
. Fig*3*2
Characteristics on 2-t plane
lines
Fig.4.1
Test pipe
Second World
2792
Congress
of Nonlinear
Analysts
4. Some examples of numerical calculation results The two cases in the boiling flow from subcooled water to saturated two-phase flow are presented. One is the steady boiling flow and the another is an oscillatory flow due to density wave instability which is caused by a little larger heating rate than in the case of steady boiling flow. Unsteady flow calculations are made for the test pipe shown in Fig.4.1, where both ends are opened and at one end pressure is kept constant 3.5 MPa and at the another end pressure is decreased gradually until 2.0 MPa. For the case of steady flow, boundary conditions of pressure and heating rate given for test pipe are shown in Fig.4.2 and Fig.4.3, respectively. The inlet temperature of water is 215’C, which means 27.5-C subcooling to the inlet pressure 3.5 MPa. As for the characteristics grid, the number of calculation points is 41 including both ends, which means the test pipe is divided into N=40, and the number of u-lines is set as NxS initially. Calculation results for the case of steady boiling flow at grid No.=30 are shown in Figs.4.4-4.8. The steady boiling flow is attained after the reach of constant heat input. We can observe it in each of pressure and velocity staying about a constant value although a very small high frequent It is not clear that this oscillation means the oscillation is contained, as shown in Figs.4.6-4.7. phenomena of boiling flow itself or a kind of numerical instability occurred in the calculation of boiling initiation region. Anyway, boiling front stays at about 3.4 m as shown in Fig.4.8. The change of void fraction and acoustic velocity are shown in Fig.4.6 and Fig.4.7, respectively. Calculation results for the case of density wave instability in a boiling flow from a subcooled water are shown in Figs.4.9-4.16. First, the heat flux change is shown in Fig.4.9. The move of characteristics grid point caused by the sonic velocity change due to boiling is shown in Fig.4.lO(a),@). The density wave instability behaviors calculated in the subcooled water region (grid No.=4) are shown in Fig.4.11 and Fig.4.12, for pressure and velocity, respectively. And the behaviors in the saturated boiling region (grid No.=30) are shown in Figs.4.13-4.15, for pressure, velocity, void fraction and move of boiling front, respectively. From these results, the followings are shown. The oscillation phases of density wave instability in the subcooled region and boiling region are almost in the same phase, except that the phase of velocity in boiling region show a small phase lag. The phase of pressure shows a 180 degree difference compared to that of velocity. The pressure profile in the flow channel shows a monotonous pressure decrease because pressure oscillation is in phase through the flow channel.
0.0
Fig.4.2
Pressure
change
at exit
Fig.4.3
3.0
External heat case of steady
6.0
flux change flow
12.0
9.0
in the
15.c
Second World
Congress
of Nonlinear
Analysts
2793
3.2
2.4
Pressure
cig.4.4
change at Grid No.=30
Fig.45 1600~’ 1200
Velocity
change at Grid No.=30
’ ’ ’ I ’ ’ ’ ’ I ’ ’ ’
I
I4--
I ’ ’ ’ ’ ! ’ ’ ’
1
.-.--I ’_....__.___ _.______ -j_____- ._______ i ______.____ 1 ; / 1
..-.-. ._..i-.... --.---i-.-_.__ j...__-._._ .i___..._.__ kc i”=(:i _ .._..__ 1 4! 1.-.-. - .i...-.-.-.--.i-.____. -.+..-~ .___ -i--‘0
Fig.45
‘Void fraction change
at Grid No.=30
orI I I '1 0.0 3.0 6.0 9.0 Fig.4.7 Acoustic velocity change
12.0
15.0
at Grid No.=30
5 4 3 2 1 0 0.0
6.0
3.0
Fig.4.8 M&e 5_,11~,II~~,I~I~,~~~~_ LT.-
.-._._ j-.---
-.-.--
15.0
Fig.4.9 External
of Boiling front
.-.-. J-. I
case s 0
I
I 4
12.0
9.0
.-.- -.j”
--.-
-..-1
.-.--
Move of grid points
(No.=5,6,7,8)
I I +,q-e * L-L tx.-~c,~ .+-a---v-
,I ---+-
u- -w-e_ t’+)’ t
.-_
4
0
Fig.e.lO(a)
5
heat flux change in ,the of density wave instability flow
.-s 4
1 0 0I.0
FIgA.lO$r
12.0
16.0
fvloves~f grid points (No.=1 3,20,30,38)
Second World
2794
Congress
of Nonlinear
Analysts
12.0 16.0 a.0 Velocity change at Grid No.=4. (Subcooled water region) 4.0
0.0
Pressure change at Grid No.=4. (Subcooled water region)
Fig.4.11
Fig.4.12 40
I!!.I,III1
I"',"" I
I I
I
5z 301____. I I 1 --.-.__ .___ i---..--...-.c-.-.----.i._-. 1 20 .--.-.-.-.-j ,o
.-
_....
__.. _._._
--
__..___
---+
___._...__
r
.__.
-._
__._
____-_,
__!
_
.-
-.
.-....
-.-.-.-.-.
Io~o
Fig.4.13
16.0' 12.0 a.0 Pressure change at Grid No.=30. (Boiling region)
0.0
;
4.0
a.0
Velocity change (Boiling region)
Fig.4.14
at Grid N0.=30.
II ! ,/’-J--l 0.5-....-..‘.-.-...., j+L-._..-.. +. ._.._._______ c_____________: ;:I !! I 0.25 -----.c-.i---.------~.~.-.--.----~.~-....-.~.~-~ c ,/ 0 - I ,, d 0.0
Fig.4.15
j I 4 , 4.0
Time (s) c , I ,z,
a.0
Void fraction change (Boiling region)
1 ,* 12.0
16.0
at Grid No=30.
0.0
4.0
Fig.4.16
8.0
12.0
Move of boiling front
16.0
SecondWorld Congress of NonlinearAnalysts
2795
5. Nomenclature a : Sonic velocity m/s r : Latent heat J&z a, : Sonic velocity for saturated two-phase flow s : Specific entropy Jikg-K A : Flow area m2 t : Time set C, : $/& used in Fig.5. T : Temperature K C, : Specific isobaric heat capacity u : Fluid flow velocity m/s J/kg-K Cy : Specific isochoric heat capacity J/kg-K V : Specific volume m3/kg d : Diameter m W : Mass flow rate kg/s e : Specific internal energy Jk X : Flow quality h : Specific enthalpy J&z z : Distance m H : Defined by Eq(3.7) a : Void fraction e : Length m p : Coefficient of thermal expansion l/K Lp : Peripheral length kg/m3 p : Density N : Number of division of test pipe 4 : Heat flux w/m2 P : Pressure Pa K : isothermal compressibility l/Pa ¶ : Heat contributing to entropy increase of unit mass Jh (subscript) E : External heat Q : Point Q f : Friction R : Point R g : Gas S : Point S Q : Liquid TP:Two-Phase mixture m : g, !Z or mixture 6. Conclusion Numerical calculation method using characteristics grid for analyzing unsteady boling flow from subcooled water is presented. (2) Caluculation results of two types of boiling flow are presented. One type is steady boiling flow which is attained after heat flux becomes constant. Another type is unsteady boiling flow with density wave instability even after constant heat flux is given. (1)
References (1) (2) (3)
(4)
Banejee, S. and Hancox, W.T., “On the Development of Methods for Analysing Transient Eh,+ Boiling,” Int. J. Multiphase Flow., Vo1.4, ~~437-460 (1978) Gidaspow, D., Rasouli, F. and Shin, Y.W., “An Unequal Velocity Model for Transient Two-Phase Flow by the Method of Characteristics,” Nucl. Sci. and Engg., 84, pp179-195 (1983) Lyczkowski, R.W., Gidaspow, D., Solbrig, C.W. and Hughes, E.D, “Characteristics and Stability Analysis of Transient One-Dimensional Two-Phase Flow Equations and Their Finite Difference Approximations,” Nucl. Sci. and Engg., 66, pp378-396 (1978) Okazaki, M., “Analysis of Density Wave Instability in a Boiling Flow Using a Characteristic Method “ASME Conference. PVP-Vo1.270, 65-73 (1994)