- Email: [email protected]
An analytical method for thermal-hydraulic transients in piping networks
441
AN ANALYTICALMETHOD FOR THERMALHYDRAULXTRANSIEIYTS Nmwcmcs Ryuhei ICAWABE EtuqpRawa&&ikmtwyHitachi,
IN PIPING
Ltd.., 1168 Monyama-cho, Hitachi -shi, Ibamkl- ken, 316 Japan
Raccivod 1 May 1982
To Predict the thermal-hydraulictransients,an analytical method has been developed for single and two-phase flow in arbitrary piping networks. In this method the piping network i8 representedbv vesselsand flow channels. The thermal-hydraulic transients in the channel are described by partial differential quationr derived from mau: momentum and energy cotutervation 18~. The partial differential quotions are solved implicitly, simultaneously for the whole nt?twork, with the ordinary differential equation8 that da&be the change of vessctpressuresand enthalpies. Numerical calculation error is evaluatedin the implicit method for the integration of partial differential equations of channel flow. In the numerical calculation M artificial diffusion appears with a diffusion coefficient Ar X2/2, where Ar is a time step and A denotes the propagation velocity of the perturbation
1. InThe cooling system ?f a nuclear power plant consists of a complex piping ncttwork, in which water, steam or their twephasc mixture transport the heat from the reactor core to iF.e turbine. To estimate thermal shocks on the cooling system components, hydraulic loads on the piping system, and temperature rises of the fuel rods during and after a power change, a reactor scram and other large accidents, computer programs are necessary that can predict thermal-hydraulic transients in the piping network, It is desirable that the computer following two features:
programs have the
(i) arbitrary time step: the time constant of the transients in the piping network can range widely, i.e. from millisecond8 in water hammer phenomena to several minute8 in power changes in the plant, The time step width of the calculation should not be limited by the stability requirement in order to avoid the expensive calculation for plow long term tranoient, and if a very
small time step is adopted, the computer program bhould give the exact solution even for the water hammer phenomena; (ii) arbitrary flow channel network: the modelhag of the system depends on the requirement to the accuracy of the results. For example the effects of parallel chanI
nel should be analyzed in the core heat up calculation, and the core may be rcprescnted by single averaged channel in whole-plant dynamic analysis,
T.A. Porsching, J.H. Murphy and J.rL Redfield (I] show a method that io applicable to an a;rbitrary hydraulic network with an implicit method in numerical integration, i.e. unlimited time step. In their method the flow network is represented by nodes and links, where mass and energy are defined in the node8 and the links represent the flow path between nodes, ThesGfore the dynamics of a system are described by ordinary differential equations. The fluid in pipe is continuous and it8 behavior is described by partial differential quations. When the piping system forms a network, these partial differential quations for each single channel mu8t be solved simultaneously. In thi8 paper, an implicit method is presented that is applicable to arbitrary networks, and estimation of the error in the numerical calculation is dascribcd. 2. Auumptions For the calculation, the hydrsulic network is divided into vessels and ‘flow channels. In each vessel the pres8ure is considered uniform. A one-dimensional expression is applied to the channel flow and thermal quilibrium &St8 between gas and liquid phases.
0029-5493/82/0000-0000/$02.75 Q 1982North-Holland
442
R. Kawab¢ / Thermai.ky&'aulic tramienu be piping ewtworks
The ¢ag,rgy conservation law is
3. Basic equ2tioas
ap(e
3.1. Pressure m a vessel
u'/2)
a~,(e
at
The pressure in a vessel is calculated by an ordinary dJferential equation and is expressed as a function of total mass and internal energy in the vessel P= F(M. U).
(l)
The change rates of mass and energy are related to the flog and the fluid enthalpy in the channels, as follows
.~ = E w,-E w,.
(2)
~' = Eh,W. -- E h , w , .
(3)
a
+
u2/2)
0uP - - - ~ - + q + psu cos O.
(7)
For two phase flow, the pressure, void fraction, and velocities of liquid and gas are the key dependent variables. Hence the mass conservation law is
(op, + (i - , , ) p , ) =
- ~(,,p,u, + ( l--
Q)PlMI}
•
(s)
The momentum conservation law for gas phase is
aap,u,
I
at
where i means the channels the fluid of which flows towards the vessel, and j means the channels where the fluid flow out from the vessel. Differentiation of eq. (I) and substitution of eqs. (2) and (3). give the ordinary differential equation to calculate the vessel pressure.
+
&
=
ae~o=u2= &
aP =~
+
=p,s cos s - F.,
(9)
and that for liquid phase is
a(I-e)p,u~_(l_a)aP
a(i-=)p,u,= at
-
az
az"
+ (! - = ) p , x c o s e - r,.
(io)
The energy conservation law is /
t
/
Steam fraction in out-flowing fluid affects the pres,,ure decrease rate and the residual mass in the vessel, so that cvaluation of the fluid steam quality is a very =mportant factor. This fluid quality, however, depends on internal structures in the vessel, such as a steam dryer or steam separator and as there is a gradient of steam void distribution the quality depends also on the height of the nozzle. In addition, the relative ~elocity of water to steam also affetts the fluid quality. Therefore, individual models are to be made for each problem for the enthalpy in out-flowing fluid. At present, it is assumed that water and steam are homogeneously mixed.
0
P,
u2/s)+
,.V2) + ( ~ - =)p,=,(h, + ,,~/2)) + q.
(l~)
It is difficult to determine frictional force Fs. F I. because those include interracial frictional force between liquid and gas. Instead of giving an explicit expression for Fa and F m. these are considered to be unknown variables, and two more equations are added. They are an equation for total frictional force on the pipe wall,
3.2. Flow in a single channel
F,+F,= ~ - Gl~, p ¢,,
(12)
The behavior of fluid in a single channel is expressed by partial differential equations derived from the conservation laws of mass, momentum and energy. For single-phase flow the mass conservation law is
and an equation giving the relation between gas-phase velocity and liquid-phase velocity, for example the slip ratio correlation:
ap ,¥ =
ua=su,.
apu az '
(5)
where spatial coordinate along pipe z and time t are the independent variables. The momentum conservation law is apu
apu 2
aP
&
az
w
Ot
+
og cos 0 - ~ D Oulul.
(6)
(13)
To simplify the expression, eqs. (5) and (1 !) are rewritten as eq. (14), using vector and matrix
A~f+e
=¢.
(14)
And eqs. (12) and (13) are also expressed in the same
K gmmbe / Tkeemal.kydraultc trmui~ts in piping nctww~
way by eq. (15)
443
expressed by the cha~ge in boundary condition:
~-E
(Is)
The main problem is to solve the eqs. (14) and (I 5) for all dmmels simultanmusly with the ordinary differential eqs. (4) for vessels.
a e v~
(19)
(ac,)-(x) aP~ +(r)"
ah°v
4. D e t w e q m t i o m 4.2. Calculation for the whole network
4.1. Calculation for single channel In the finite difference approximation to the differential equations, the implicit method is used to keep the stabifity for the arbitrary size of time step. In this method the differentiation with respect to spatial coordinate is approximated with the difference expression, using the values at each new time step
The state of the whole network is calculated by the following method- The changes in pressure and enthalpy for all vessels are given by finest combinations of d¢~ at the channel inlet and outlet as shown in eq. (20), which is derived from the ordinary differential equations (3),
(4), \
( AP, l
( W)(
Ahv /
Ot
'-' "
(17)
The state in the channel is calculated with the boundary conditions, which are given from vessel pressures and enthalpies at the inlet and outlet. Since the pressures and the enthalpies are unknown at this stage, the increments of these variables are involved in the linearized simultaneous equations, which are derived from eqs. (14) and (15) by substitution of eqs. (16) and (17)
ah~ ,.',
'..,
0
: ,,.8.,'t, I
-(s).
,k
/0 i x /
+
().V
(20)
At '
a~, (¢+a~,)'÷'-(,+a,~) o-7 = 2a:
]
A4, ),~,.,/o.,,.,
.0
",
"t ....
i
(IS)
The increments of variables d¢~ at the channel inlet and outlet are given in eq. (19). That is,
Ap' ( A~,),.,,/.,.,.
(x)
ah~
Apv° +(y). ah o
(21)
where x and y are submatrices of X and Y. The increments of vessel pressures and enthalpies are obtained by solving eqs. (20) and (21) simultaneously. The rank of the coefficient matrix is the order of number of vessels and number of channels, which is far less than the total node number in the channels. After vessel variables dPv, dh v are determined, the channel state is obtained by substitution of these into eq. (19).
0
i dPv
in the simultaneous eqs. (18), the number of unknown variables is larger than that of the equations. After transposing the variables of boundary conditions and using Gaussian elimination, which is easily performed as the coefficients are in a band matrix, eq. (I 8) becomes eq. (19), where the change in channel state is
5. Error estimation As described above, the whole =alculat:.on consists of three parts, namely ordinary differential equations for vessels, partial differential equations for channels, and matrix calculation for the whole tietwork. Among these three calculations, the error for the partial differential equation is of greatest interest and importance. Ignoring the terms of small effects, the partial differential equation (14) can be transformed to three
R. gawabe I T S e r m ~ - ~
444
6, s
independent equations. a~,,, = _ X a~,. Bt " Oz "
(22)
where ~,, is the n-th eigen value of the matrix A - aB, and 4 ' . is the n-th elements of the vector 1/' defined from 4~ as follows 4' = U¢,.
(23)
where U is a matrix formed by the array of eigen vectors of A - tB. The eigen values X,, are pressure wave propagation velocities ( u + c ) , and fluid velocity or void propagation velocity if slip exist in two phase flow. The ~, are pressures P + pcu and fluid specific i~tternal energy i.e. temperature or void fraction. Eq. (22) is rewritten to the expression on the characteristic curve. dz _- X, dt
d~, = 0 , dt
(24)
where n is omitted. F_.q. (24) implies that p e r t u r b a t i o n is
transferred along the characteristic curve without shape change. Expanding the variable 4' at time t in Taylor series with respect to z and substituting into eq. (23), a new variable at advanced time is obtained, where the differentiation is approximated by the implicit method as described in eqs. (16) and (17) 4 ' ( t ) = 4'0 + 4',z + 4'2z ~ + . . .
~(t
+ ~t)=
4'o
(25)
+ ~2~t 4'2 + "'"
+ z'(4', + 3~2At24'.~ - 4 X A z 2 A t 4'4 + " ' " + z,2(4'2 + " ' "
mm~m~ m ~ o m g ~ . ~ t a ~
(26)
r ---~.L L
AA.
= ~'2At
2
~'~.
---~
02~
(27)
a(:,);"
!.0
~.
"-,,.o3s ~,.,,~o,s
................... ~ . . 1 . _ 7 ~ --'.~,..-. -
~
:E ~l~ i
i
s--5
r'"....o1~ "-,.~o2s
According to eq. (24), ~,(t. z) equals ~ ( t + d t , z'). Comparing eq. (25) and (26). the following equation is derived at
~
As an example of the implicit m e t l ~ for a sinsk channel flow, calculaticca for a p ~ wave propagation in two-pha~ flow is ~ in fill. I. A 20 m fan8 straisht pipe is connected to a vessel at the leftend. The friction of the pipe is nesligibly small The system is initially kept at 1 MPa and the void fraction in the pipe is 0.64. The vessel pressure rises in a eosiae function within 0.2 s. Three cases of pressure profile are shown with a different slip ratio, comtant in each calculation. The time step and spatial mesh size are 2 ms and 0.8 m respectively. The remits show the pressure wave propagation with the higher velocity for the larger slip. The calculated results shows slight diffusion of the wave front. Fig. 2 shows the results of other calculations for the same geometry with different time steps. The broken lines show the solution of the diffusion equation of which the diffusion coefficient is determined from eq. (27). The results of the computer program agree well with the broken line. Fig. 3 demonstrates the calculation for a piping network consisting of 3 vessels and 4 pipes with one branching tee, as shown in fig. 3a. Initially the system is kept at 10 MPa and the steam mass fractions in vessel Vt, V2 and V3 are 0.7, 0.1, 005 respectively. Sudden pressure drop is given to the exit. The friction of the piping wall is calculated by using Chisholm's correlation [2] and the slip between liquid and gas phase is ignored. Transients of mass flow rate in the pipe channels are shown in fig. 3b and the transients of pressure in the vessels in fig. 3c. At first, the pressure of vessel V3
where z' = = +
m
=.
S= 7.
t 0 [r... -
. "%. ......
: ~ "~%_ '--":~"~"
........
O. 1.1[_ |
S- 1
This equation signifies that the perturbation diffuses on the characteristic curve with a diffusion coefficient ~2At/2. 0
-
5 Distance
from
10 Vessel
15 ( m )
-
20
Fig. I. Calculation of pressure wave proPalpttion in tw~pha~ flow.
I TAm,~.~c.n.~
R. ~ 1.1 ~
445
The calculation was performed with a time step 0.05 s and the number of total nodes on pipes is 58. The CPU time is 12 S by the compuer IBM 3033.
by ComputerProgram
~
in piplng ~
by Diffusion Eq. •
~
~
-
7. Condmiom
g $
1.0
'st=O'O02i~ . . . . I
i
•
*
I
O
5
i
l
i
I
i
Diitznce (m)
I
I
*
lO
,
15
F~g. 2. Calculation error evaluation for preuare wave in twophase flow (t - 0.4 s, s - !).
follows the exit pressure drop, as V3 is the nearest to the exit. After the discharge flow from V2 established at 0.6 s, as the volume of V3 is small, the V3 pressure is controlled by the branching tee pressure, which is determined from pressure drop in the discharge line. Therefore, the flow from V3 becomes very low and the pressure decreases slower.
An analytical method has been developed for calculating thermal-hydraulic transients of single and twophase flow in arbitrary piping networks. In the method, the piping network is represented by flow channet~ and vessels. The thermal-hydraulic transients in tb:- channel are expressed by partial differential equations, which are solved by the implicit method simultaneously for the whole network with the ordinary differential equations that describe the change of vessel pressures and enthalpies. Numerical calculation error is evaluated in the implicit method for solving the partial differential equations of channel flow. The artificial diffusion of perturbation is superposed on the propagation of the perturbation. The diffusion coefficient is calculated to be Ark2~2, where At is a time step and A denotes the propagation velocity of the perturbation.
Acknowledgement I wish to express my gratitude to Dr. Norihiko Sagawa and Akira Suzuoki for their valuable advice and consultation.
Om
Nomenclature 51~s A-OOlr~ 20m
IOm Pi.p.i Wx
(I) Connectionof Flew Channelsand Vessels
Wx
M
,_ •
Time Is) (b) Transients of Flow hte in Channels
~ .
c D e F G g h
t
o'..,, Time (sl {c) Pressuresof Yessels and BoundaryC~ditton
Fig. 3. Calculation example for piping system.
P q s t
U u W z
velocity of pressure wave hydraulic diameter specific internal energy frictional force mass velocity gravitational acceleration specific enthalpy mass pressure heat transferred to fluid per unit volume and unit time slip ratio time internal energy fluid velocity mass flow rate spatial coordinate
446
a 0 0 0
R. Kawabe / T & ~ n n a l . h ~ c teans~u inpipinSnaworks
void fraction density angle betwe¢'~ flow channel and vertical fine z two-phase fricton multiplier
Superscripts
i ! o
spatial node inlet of flow channel outlet of flow channel
Subscrlpt.~ g gas phase I liquid phase v vessel
Rdalees [I] T.A. ~ b J.H. Murphy and J.A. Rcdfield, Stable num4~ical inlqp'alion of conservation oquations for hydraulic netwocks, Nucl. Sci. Enlp3. 43 (1972) 218-225. 121 g. Akq~wa, Gas-Liquid Two-Phase Flow (Corona Publishing Comp. Ltd.. Tokyo, 1974) p. 82.
AN ANALYTICALMETHOD FOR THERMALHYDRAULXTRANSIEIYTS Nmwcmcs Ryuhei ICAWABE EtuqpRawa&&ikmtwyHitachi,
IN PIPING
Ltd.., 1168 Monyama-cho, Hitachi -shi, Ibamkl- ken, 316 Japan
Raccivod 1 May 1982
To Predict the thermal-hydraulictransients,an analytical method has been developed for single and two-phase flow in arbitrary piping networks. In this method the piping network i8 representedbv vesselsand flow channels. The thermal-hydraulic transients in the channel are described by partial differential quationr derived from mau: momentum and energy cotutervation 18~. The partial differential quotions are solved implicitly, simultaneously for the whole nt?twork, with the ordinary differential equation8 that da&be the change of vessctpressuresand enthalpies. Numerical calculation error is evaluatedin the implicit method for the integration of partial differential equations of channel flow. In the numerical calculation M artificial diffusion appears with a diffusion coefficient Ar X2/2, where Ar is a time step and A denotes the propagation velocity of the perturbation
1. InThe cooling system ?f a nuclear power plant consists of a complex piping ncttwork, in which water, steam or their twephasc mixture transport the heat from the reactor core to iF.e turbine. To estimate thermal shocks on the cooling system components, hydraulic loads on the piping system, and temperature rises of the fuel rods during and after a power change, a reactor scram and other large accidents, computer programs are necessary that can predict thermal-hydraulic transients in the piping network, It is desirable that the computer following two features:
programs have the
(i) arbitrary time step: the time constant of the transients in the piping network can range widely, i.e. from millisecond8 in water hammer phenomena to several minute8 in power changes in the plant, The time step width of the calculation should not be limited by the stability requirement in order to avoid the expensive calculation for plow long term tranoient, and if a very
small time step is adopted, the computer program bhould give the exact solution even for the water hammer phenomena; (ii) arbitrary flow channel network: the modelhag of the system depends on the requirement to the accuracy of the results. For example the effects of parallel chanI
nel should be analyzed in the core heat up calculation, and the core may be rcprescnted by single averaged channel in whole-plant dynamic analysis,
T.A. Porsching, J.H. Murphy and J.rL Redfield (I] show a method that io applicable to an a;rbitrary hydraulic network with an implicit method in numerical integration, i.e. unlimited time step. In their method the flow network is represented by nodes and links, where mass and energy are defined in the node8 and the links represent the flow path between nodes, ThesGfore the dynamics of a system are described by ordinary differential equations. The fluid in pipe is continuous and it8 behavior is described by partial differential quations. When the piping system forms a network, these partial differential quations for each single channel mu8t be solved simultaneously. In thi8 paper, an implicit method is presented that is applicable to arbitrary networks, and estimation of the error in the numerical calculation is dascribcd. 2. Auumptions For the calculation, the hydrsulic network is divided into vessels and ‘flow channels. In each vessel the pres8ure is considered uniform. A one-dimensional expression is applied to the channel flow and thermal quilibrium &St8 between gas and liquid phases.
0029-5493/82/0000-0000/$02.75 Q 1982North-Holland
442
R. Kawab¢ / Thermai.ky&'aulic tramienu be piping ewtworks
The ¢ag,rgy conservation law is
3. Basic equ2tioas
ap(e
3.1. Pressure m a vessel
u'/2)
a~,(e
at
The pressure in a vessel is calculated by an ordinary dJferential equation and is expressed as a function of total mass and internal energy in the vessel P= F(M. U).
(l)
The change rates of mass and energy are related to the flog and the fluid enthalpy in the channels, as follows
.~ = E w,-E w,.
(2)
~' = Eh,W. -- E h , w , .
(3)
a
+
u2/2)
0uP - - - ~ - + q + psu cos O.
(7)
For two phase flow, the pressure, void fraction, and velocities of liquid and gas are the key dependent variables. Hence the mass conservation law is
(op, + (i - , , ) p , ) =
- ~(,,p,u, + ( l--
Q)PlMI}
•
(s)
The momentum conservation law for gas phase is
aap,u,
I
at
where i means the channels the fluid of which flows towards the vessel, and j means the channels where the fluid flow out from the vessel. Differentiation of eq. (I) and substitution of eqs. (2) and (3). give the ordinary differential equation to calculate the vessel pressure.
+
&
=
ae~o=u2= &
aP =~
+
=p,s cos s - F.,
(9)
and that for liquid phase is
a(I-e)p,u~_(l_a)aP
a(i-=)p,u,= at
-
az
az"
+ (! - = ) p , x c o s e - r,.
(io)
The energy conservation law is /
t
/
Steam fraction in out-flowing fluid affects the pres,,ure decrease rate and the residual mass in the vessel, so that cvaluation of the fluid steam quality is a very =mportant factor. This fluid quality, however, depends on internal structures in the vessel, such as a steam dryer or steam separator and as there is a gradient of steam void distribution the quality depends also on the height of the nozzle. In addition, the relative ~elocity of water to steam also affetts the fluid quality. Therefore, individual models are to be made for each problem for the enthalpy in out-flowing fluid. At present, it is assumed that water and steam are homogeneously mixed.
0
P,
u2/s)+
,.V2) + ( ~ - =)p,=,(h, + ,,~/2)) + q.
(l~)
It is difficult to determine frictional force Fs. F I. because those include interracial frictional force between liquid and gas. Instead of giving an explicit expression for Fa and F m. these are considered to be unknown variables, and two more equations are added. They are an equation for total frictional force on the pipe wall,
3.2. Flow in a single channel
F,+F,= ~ - Gl~, p ¢,,
(12)
The behavior of fluid in a single channel is expressed by partial differential equations derived from the conservation laws of mass, momentum and energy. For single-phase flow the mass conservation law is
and an equation giving the relation between gas-phase velocity and liquid-phase velocity, for example the slip ratio correlation:
ap ,¥ =
ua=su,.
apu az '
(5)
where spatial coordinate along pipe z and time t are the independent variables. The momentum conservation law is apu
apu 2
aP
&
az
w
Ot
+
og cos 0 - ~ D Oulul.
(6)
(13)
To simplify the expression, eqs. (5) and (1 !) are rewritten as eq. (14), using vector and matrix
A~f+e
=¢.
(14)
And eqs. (12) and (13) are also expressed in the same
K gmmbe / Tkeemal.kydraultc trmui~ts in piping nctww~
way by eq. (15)
443
expressed by the cha~ge in boundary condition:
~-E
(Is)
The main problem is to solve the eqs. (14) and (I 5) for all dmmels simultanmusly with the ordinary differential eqs. (4) for vessels.
a e v~
(19)
(ac,)-(x) aP~ +(r)"
ah°v
4. D e t w e q m t i o m 4.2. Calculation for the whole network
4.1. Calculation for single channel In the finite difference approximation to the differential equations, the implicit method is used to keep the stabifity for the arbitrary size of time step. In this method the differentiation with respect to spatial coordinate is approximated with the difference expression, using the values at each new time step
The state of the whole network is calculated by the following method- The changes in pressure and enthalpy for all vessels are given by finest combinations of d¢~ at the channel inlet and outlet as shown in eq. (20), which is derived from the ordinary differential equations (3),
(4), \
( AP, l
( W)(
Ahv /
Ot
'-' "
(17)
The state in the channel is calculated with the boundary conditions, which are given from vessel pressures and enthalpies at the inlet and outlet. Since the pressures and the enthalpies are unknown at this stage, the increments of these variables are involved in the linearized simultaneous equations, which are derived from eqs. (14) and (15) by substitution of eqs. (16) and (17)
ah~ ,.',
'..,
0
: ,,.8.,'t, I
-(s).
,k
/0 i x /
+
().V
(20)
At '
a~, (¢+a~,)'÷'-(,+a,~) o-7 = 2a:
]
A4, ),~,.,/o.,,.,
.0
",
"t ....
i
(IS)
The increments of variables d¢~ at the channel inlet and outlet are given in eq. (19). That is,
Ap' ( A~,),.,,/.,.,.
(x)
ah~
Apv° +(y). ah o
(21)
where x and y are submatrices of X and Y. The increments of vessel pressures and enthalpies are obtained by solving eqs. (20) and (21) simultaneously. The rank of the coefficient matrix is the order of number of vessels and number of channels, which is far less than the total node number in the channels. After vessel variables dPv, dh v are determined, the channel state is obtained by substitution of these into eq. (19).
0
i dPv
in the simultaneous eqs. (18), the number of unknown variables is larger than that of the equations. After transposing the variables of boundary conditions and using Gaussian elimination, which is easily performed as the coefficients are in a band matrix, eq. (I 8) becomes eq. (19), where the change in channel state is
5. Error estimation As described above, the whole =alculat:.on consists of three parts, namely ordinary differential equations for vessels, partial differential equations for channels, and matrix calculation for the whole tietwork. Among these three calculations, the error for the partial differential equation is of greatest interest and importance. Ignoring the terms of small effects, the partial differential equation (14) can be transformed to three
R. gawabe I T S e r m ~ - ~
444
6, s
independent equations. a~,,, = _ X a~,. Bt " Oz "
(22)
where ~,, is the n-th eigen value of the matrix A - aB, and 4 ' . is the n-th elements of the vector 1/' defined from 4~ as follows 4' = U¢,.
(23)
where U is a matrix formed by the array of eigen vectors of A - tB. The eigen values X,, are pressure wave propagation velocities ( u + c ) , and fluid velocity or void propagation velocity if slip exist in two phase flow. The ~, are pressures P + pcu and fluid specific i~tternal energy i.e. temperature or void fraction. Eq. (22) is rewritten to the expression on the characteristic curve. dz _- X, dt
d~, = 0 , dt
(24)
where n is omitted. F_.q. (24) implies that p e r t u r b a t i o n is
transferred along the characteristic curve without shape change. Expanding the variable 4' at time t in Taylor series with respect to z and substituting into eq. (23), a new variable at advanced time is obtained, where the differentiation is approximated by the implicit method as described in eqs. (16) and (17) 4 ' ( t ) = 4'0 + 4',z + 4'2z ~ + . . .
~(t
+ ~t)=
4'o
(25)
+ ~2~t 4'2 + "'"
+ z'(4', + 3~2At24'.~ - 4 X A z 2 A t 4'4 + " ' " + z,2(4'2 + " ' "
mm~m~ m ~ o m g ~ . ~ t a ~
(26)
r ---~.L L
AA.
= ~'2At
2
~'~.
---~
02~
(27)
a(:,);"
!.0
~.
"-,,.o3s ~,.,,~o,s
................... ~ . . 1 . _ 7 ~ --'.~,..-. -
~
:E ~l~ i
i
s--5
r'"....o1~ "-,.~o2s
According to eq. (24), ~,(t. z) equals ~ ( t + d t , z'). Comparing eq. (25) and (26). the following equation is derived at
~
As an example of the implicit m e t l ~ for a sinsk channel flow, calculaticca for a p ~ wave propagation in two-pha~ flow is ~ in fill. I. A 20 m fan8 straisht pipe is connected to a vessel at the leftend. The friction of the pipe is nesligibly small The system is initially kept at 1 MPa and the void fraction in the pipe is 0.64. The vessel pressure rises in a eosiae function within 0.2 s. Three cases of pressure profile are shown with a different slip ratio, comtant in each calculation. The time step and spatial mesh size are 2 ms and 0.8 m respectively. The remits show the pressure wave propagation with the higher velocity for the larger slip. The calculated results shows slight diffusion of the wave front. Fig. 2 shows the results of other calculations for the same geometry with different time steps. The broken lines show the solution of the diffusion equation of which the diffusion coefficient is determined from eq. (27). The results of the computer program agree well with the broken line. Fig. 3 demonstrates the calculation for a piping network consisting of 3 vessels and 4 pipes with one branching tee, as shown in fig. 3a. Initially the system is kept at 10 MPa and the steam mass fractions in vessel Vt, V2 and V3 are 0.7, 0.1, 005 respectively. Sudden pressure drop is given to the exit. The friction of the piping wall is calculated by using Chisholm's correlation [2] and the slip between liquid and gas phase is ignored. Transients of mass flow rate in the pipe channels are shown in fig. 3b and the transients of pressure in the vessels in fig. 3c. At first, the pressure of vessel V3
where z' = = +
m
=.
S= 7.
t 0 [r... -
. "%. ......
: ~ "~%_ '--":~"~"
........
O. 1.1[_ |
S- 1
This equation signifies that the perturbation diffuses on the characteristic curve with a diffusion coefficient ~2At/2. 0
-
5 Distance
from
10 Vessel
15 ( m )
-
20
Fig. I. Calculation of pressure wave proPalpttion in tw~pha~ flow.
I TAm,~.~c.n.~
R. ~ 1.1 ~
445
The calculation was performed with a time step 0.05 s and the number of total nodes on pipes is 58. The CPU time is 12 S by the compuer IBM 3033.
by ComputerProgram
~
in piplng ~
by Diffusion Eq. •
~
~
-
7. Condmiom
g $
1.0
'st=O'O02i~ . . . . I
i
•
*
I
O
5
i
l
i
I
i
Diitznce (m)
I
I
*
lO
,
15
F~g. 2. Calculation error evaluation for preuare wave in twophase flow (t - 0.4 s, s - !).
follows the exit pressure drop, as V3 is the nearest to the exit. After the discharge flow from V2 established at 0.6 s, as the volume of V3 is small, the V3 pressure is controlled by the branching tee pressure, which is determined from pressure drop in the discharge line. Therefore, the flow from V3 becomes very low and the pressure decreases slower.
An analytical method has been developed for calculating thermal-hydraulic transients of single and twophase flow in arbitrary piping networks. In the method, the piping network is represented by flow channet~ and vessels. The thermal-hydraulic transients in tb:- channel are expressed by partial differential equations, which are solved by the implicit method simultaneously for the whole network with the ordinary differential equations that describe the change of vessel pressures and enthalpies. Numerical calculation error is evaluated in the implicit method for solving the partial differential equations of channel flow. The artificial diffusion of perturbation is superposed on the propagation of the perturbation. The diffusion coefficient is calculated to be Ark2~2, where At is a time step and A denotes the propagation velocity of the perturbation.
Acknowledgement I wish to express my gratitude to Dr. Norihiko Sagawa and Akira Suzuoki for their valuable advice and consultation.
Om
Nomenclature 51~s A-OOlr~ 20m
IOm Pi.p.i Wx
(I) Connectionof Flew Channelsand Vessels
Wx
M
,_ •
Time Is) (b) Transients of Flow hte in Channels
~ .
c D e F G g h
t
o'..,, Time (sl {c) Pressuresof Yessels and BoundaryC~ditton
Fig. 3. Calculation example for piping system.
P q s t
U u W z
velocity of pressure wave hydraulic diameter specific internal energy frictional force mass velocity gravitational acceleration specific enthalpy mass pressure heat transferred to fluid per unit volume and unit time slip ratio time internal energy fluid velocity mass flow rate spatial coordinate
446
a 0 0 0
R. Kawabe / T & ~ n n a l . h ~ c teans~u inpipinSnaworks
void fraction density angle betwe¢'~ flow channel and vertical fine z two-phase fricton multiplier
Superscripts
i ! o
spatial node inlet of flow channel outlet of flow channel
Subscrlpt.~ g gas phase I liquid phase v vessel
Rdalees [I] T.A. ~ b J.H. Murphy and J.A. Rcdfield, Stable num4~ical inlqp'alion of conservation oquations for hydraulic netwocks, Nucl. Sci. Enlp3. 43 (1972) 218-225. 121 g. Akq~wa, Gas-Liquid Two-Phase Flow (Corona Publishing Comp. Ltd.. Tokyo, 1974) p. 82.