7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Two-fluid model for 1D simulations of water hammer induced by condensation of hot

vapor on the horizontally stratified flow

I. Tiselj* and C. Samuel Martin**

Reactor Engineering Division, Jo2ef Stefan Institute,

Jamova 39, 1000 Ljubljana, Slovenia

iztok.tiselj@ijs.si

**School of Civil and Environmental Engineering, Georgia Institute of Technology,

Atlanta, Georgia 30332, USA

csammartin@ comcast.net

Keywords: condensation induced water hammer, 1D two-fluid model

Abstract

Condensation-induced water hammer has been studied on the experimental device described by Martin et. al. (2007). Transient

flow was initiated by injection of hot ammonia vapor into a 6 m long horizontal pipe with 146.3 mm diameter partially filled with

cold stagnant liquid. Ammonia was used as the working liquid with initial liquid temperatures between -220C and -480C and hot

vapor temperatures around 150C. After the valve opening, hot vapor enters the horizontal test section and (in most cases) a slug

is being created near the open end of the test section. Condensation of the gas bubble captured behind the slug and the closed end

of the test section accelerates the slug and a pressure surge of up to 5 MPa is being observed when the slug hits the closed end of

the pipe. In the present paper, several experimental runs are simulated with 1D computer code WAHA (Tiselj et. al. 2004).

WAHA was developed and tested for the column separation water hammer phenomena. Recent development of the code is being

performed in the field of condensation-induced water hammer in horizontal pipes. WAHA code is based on one-dimensional

six-equation two-fluid model with correlations for heat, mass, momentum transfer in stratified and dispersed flow. Hyperbolic

partial differential equations are solved with a second-order accurate numerical scheme based on high-resolution shock

capturing schemes well known in aerodynamics. Dispersed flow can be bubbly or droplet. Transition from stratified to dispersed

(bubbly or droplet) flow is based on Taitel-Dukler correlation. For the purpose of the present study WAHA was upgraded with

ammonia thermophysical tables and Chato-Dobson correlation for heat and mass transfer in stratified flow. This correlation

includes condensation of the ammonia gas on the cold pipe wall, which is dominant in the first moments of the transient. Several

modifications of the transition criteria from stratified to non-stratified flow are being tested as well as correlations for heat, mass

and momentum transfer near the slug. Calculations are compared with pressure measurements. Results of the simulation

represent a range of the 1D two-fluid model being used in the field of stratified-to-slug flow transition with heat and mass

transfer.

Introduction

Condensation-induced water hammer research at Jo2ef

Stefan Institute has been performed as a part of the EU

research project NURESIM (NUclear REactor SImulations)

and is being continued within the currently running EU

project NURISP (NUclear Reactor Simulation Platform).

Research is related to the simulations of stratified flows in

horizontal pipes. Main attention within NURESIM project

was paid to the CIWH scenario, where cold liquid was slowly

flooding a horizontal pipe filled with hot steam (Strubelj et.

al., 2010). Another type of CIWH scenario assumes injection

of hot steam into horizontal pipe partially filled with cold

liquid: this scenario is the main topic of JSI research within

NURISP project.

1) First type of CIWH can appear when the pipe filed with

hot steam is slowly flooded with cold water. This type of the

CIWH was shown to be a stochastic and thus very

unpredictable phenomena (Bjorge and Griffith, 1984,

Strubelj et. al., 2010). Models used to simulate Type 1 of the

CIWH were 1D two-fluid models and various 3D codes also

based on two-fluid models. Calculated results were

compared to experimental data from Hungarian

KFKI-PMK2 device. It was shown that neither 1D model nor

3D CFD model could accurately predict where and when the

slug will form, or if the slug will form at all.

2) Second type of CIWH can appear when hot steam

enters a pipe that is partially filled with cold liquid. The most

unstable part of the interface is always near the steam inlet

into the pipe, thus, it is more or less known where the slug

will be born. Experimental results for this type of the CIWH

were obtained by C.S. Martin (Georgia Tech, 2007). Since

the position of the slug formation is known, this type of

CIWH is less stochastic and more predictable than the type I

Paper No

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

of the CIWH. In the present paper computer code WAHA is

used as a platform for 1D simulations. The code is solving 6

equation 1D two-fluid model. WAHA code is using a special

type of numerical scheme that aims particularly at the fast

transients (like CIWH). Models that require verification and

validation are criteria that trigger transition from horizontally

stratified into slug regime of two-phase flow. These models

are required in 1D two-fluid models where horizontally

stratified flows are described with different set of

correlations than dispersed flows. The existing WAHA code

physical models were upgraded with correlations for steam

condensation on the cold walls of the horizontal pipes.

Inter-phase heat, mass and momentum exchange models in

horizontally stratified flows are already available in WAHA.

However these models require fine tuning for accurate

description of the CIWH with vapor injection.

Nomenclature

pipe c

pipe c

pipe c

interfa

matrix

matrix

ross-section [m2]

ross-section change due to elasticity [m2]

ross-section covered by phase k [m2]

cial area concentration

- temporal derivatives

- spatial derivatives

A

Ae

Ak

agf

A

B

C,

CD

Cvm

CVMI

D

d

do

E

e

Ffwall

Fk

fk

g

Htk

hk*

kk

Nuk

P

P.

Re

Qk

S

SR

Tk

t

Uk

v

Verit

Vr

X

Xsat

x

Greek l

a vapor volume fraction [m3/m3]

Cld~p

ahV

modified vapor volume fraction [m3/m3]

modified vapor volume fraction [m3/m3]

vector independent variables

p density [kg/m3]

Tg

0

0

Subsripts

f li

vapor source term [kg/(s m )]

inclination of the pipe

relaxation time

quid

g vapor

k liquid or vapor

m mixture

s saturation

i interface

Experimental Facility

An apparatus was designed by Martin et.al. (2007) to

simulate an industrial environment whereby ammonia liquid

is standing in a partially-filled horizontal pipe in thermal

equilibrium with ammonia gas above it. The essential

elements of the test setup consist of a horizontal pipe and a

high pressure tank containing hot ammonia gas as shown in

Fig. 1. The test pipe was a nominal 150 mm diameter, 6 m

long schedule 80 carbon steel pipe, having internal

diameter of 146.3 mm, and wall thickness 11 mm. The

pressure tank contained ammonia gas on top of liquid in

thermal equilibrium at ambient conditions inasmuch as the

entire test facility was outdoors. Between the pressure tank

and test pipe were three valves angle valve, solenoid valve,

and throttle valve -- and a metering orifice. The angle valve

remained fully open, while flow was initiated by a solenoid

valve for a given position of the throttle valve. The flow of

hot gas was controlled by manually positioning the throttle

valve for the existing ambient hot gas pressure. The ammonia

in the insulated test pipe was introduced from an ancillary

system containing a compressor, an auxiliary tank, and

another tank for purging non-condensible gases. For each test,

care was exercised to transfer ammonia liquid to or from the

test pipe to establish the desired depth and equilibrium

temperature. The principal measurements were (1) receiver

gas pressure (ambient temperature), (2) orifice-metered gas

flow (up to 0.45 kg/s), (3) static temperature and pressure of

saturated liquid and gas in test section (225-250 K, -0.45-1.6

bar), (4) dynamic gas pressures and shock pressures in test

section that reached up to 50 bar.

Fig. 2 shows schematic description of the

condensation-induced water hammer phenomena. Top of the

Fig. 2 shows initial conditions in the system cooled to

225-250 K. Initial depth of the liquid and initial temperature

and corresponding saturation pressure were modified in

different experimental runs. When the valve on the hot gas

inlet pipe is opened, hot gas enters the test section and can

induce slugging if the relative velocity between phases is

high enough. Gas velocity above the liquid interface depends

on inlet mass flow rate and on the condensation rate of the

gas on the cold walls of the pipe and cold surface of the liquid.

Once the slug is formed, very efficient heat transfer at the

head of the slug causes pressure difference between the tail

and the head of the slug and pushes the slug towards the

closed end. During the slug propagation the mass of the

liquid inside the slug grows and the water hammer pressure

Paper No

inter-phase drag coefficient [kg/m4]

dimensionless interfacial friction coefficient

virtual mass coefficient [kg/m3]

virtual mass term [N/m3]

pipe diameter [m]

pipe wall thickness [m]

average droplet diameter [m]

pipe elasticity modulus of the material [N/m2]

specific total energy [J/kg]

wall friction forces [N/m3]

drag force on phase k [N/m3]

dimensionless friction factor

gravity [m/s2]

volumetric heat transfer coefficient [W/(m3 K)]

specific internal enthalpy [J/kg]

thermal conductivity [W/(m K))

Nusselt number

pressure [Pa]

interfacial pressure [Pa]

Reynolds number

volumetric heat flux to phase k [W/m3]

stratification factor

source term vector

non-relaxation source term vector

relaxation source term vector

temperature [K]

time [s]

specific internal energy of phase k [J/kg]

velocity [m/s]

critical velocity [m/s]

relative velocity [m/s]

vapor quality

vapor quality at saturation

spatial coordinate [m]

letters

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

surge is registered when the slug hits the closed pipe end.

TRANSDUCERS

Static Pressure: PT(Gage); Pl; P2; PACE1; PACE2; PACE3; PACE4; PACE5

Dynamic Pressure: PCB1; PCB2; PCB3; PCB4

Static Temperature: RTD1; RTD2; RTD3; RTD4; RTD5

Dynamic Temperature: Thermocouples TC1; TC2; TC3; TC4; TC5

Fig. 1: Schematic of test pipe, orifice and pressure tank

I hot steam inlet

closed end

slug head

Fig. 2: Typical transient: red hot gas, blue cold liquid, top -

initial state, bottom slug.

Two-Phase Flow Model Of The WAHA Code

Mathematical model of the WAHA code is ID six-equation

two-fluid model similar to the models of RELAP5 (Carlson

et al, 1990) or CATHARE (Bestion, 1990) computer codes.

The basic equations are mass, momentum and energy

balances for vapor and liquid, with terms for pipe elasticity

and without terms for wall-to-fluid heat transfer. Continuity

equations for liquid and vapor (gas) phase are:

(1-a) + -p, ap +( I(1-a) p,v,

8a(t i-a) p + x

8t 8t 8x

OA(1-a) p v

St

8A(1-a)p pv/ ap

+A(1-a)--

8 x 8x

aa

A*CVM Ap, =AC, v, v,

8x

-AF v, +A(1-a)p, gcos0 -AF,

SAa pv, Aa pggv +Aap +

at 8x ax

aa

+A* CVM+Ap,- =-AC, Iv,v, +

ax

+AFgv, +Aap,gcosO-AFg,,

Internal energy balance equations for both phases are:

u 8 u f uf 9a p

ot ox 8t 8t

+(1 a) s--t I a) Pf Vx p-+p(l -)K +

a(l-a)vf ap

+P- +p(l-a)vfK K=

8x 8x

1 dA(x)

Qf -r,(u -u)+v,F,w-(1 -a)vpA d(x

A(x) dx

8ug 8ug 8a 8p aavg

a P a --+ap v--+p--+paK +p- a+

S8t 8x 8t 8t 8x

+pavgK d=

dA(x)

g g g g+vF,^ -avp dx

A(x) dx

+(1-a) pfvK f

ax

1 dA(x)

-Fg -(1-a) Pv Ax) dx

A(x) dx

(1) Specific total energy of liquid or gas is:

aapg a p +aapgvg ap

--ps +ap K +--p-v+ap, v,K 1P

at a 8t 8x x sx

g 1 dA(x)

S A(x) dx (2)

where the temporal changes of cross-section A(x,t) are

neglected in the denominators of the last term of equations

(1) and (2) (and also in Eqs. (5) and (6)). Momentum balance

equations for both phases are:

e=u+v2/2 (7)

Differential terms are collected on the left-hand side of the

equations and the non-differential terms are collected on the

right-hand side.

Terms that include constant K in Eqs. from (1) to (6) are

due to the elasticity of the pipe walls. According to Wylie and

Streeter (1978), speed c of a small pressure wave in 1D

elastic pipe filled with single-phase fluid should be reduced.

The following modification is thus introduced. Pipe

Paper No

Solenoid HEV

Angle Volve

Paper No

cross-section A(x,t) can vary along the coordinate

function of initial pipe geometry A(x) and due to the

change A,(p(x,t)).

A(x,t) = A(x) + A (p(x,t))

The pressure pulse changes the pipe cross-sec

accordance with linear relation:

dAe D dp

=K-dp

A(x,t) d E

where D is diameter, d is thickness and E is Young's n

of the elasticity of the pipe.

Closure relations

The WAHA code uses several different;

non-differential closure relations. Closure relati

two-phase flow are used to describe interfacial hea

and momentum transfer, wall friction, interfacial p

virtual mass term, equations of state etc. The equa

state for each phase k, where k is f for fluid and g fo

are:

dp,-(a P

d =- --

aP d,

Derivatives on the right-hand side of the Eq. (

determined by the ammonia property subroutines de

for the WAHA code using pressure and tempera

specific internal energy as input. Ammonia propel

pre-tabulated with subroutines developed on the I

IAPWS recommendations (Bukes, Dooley, 2001), an

at approximately 300 pressures (0 bar 2500 bar)

temperatures (195.5 K 1714 K).

The virtual mass term CVM in Eqs. (3) and (4) is

obtain hyperbolicity of the system:

av, t ax at gv

Value of the coefficient Cv, was tuned to ens

hyperbolicity of the two-fluid model equations.

virtual mass term does not ensure uncon

hyperbolicity of the equations. For very large

velocities (comparable to sonic velocity) c

eigenvalues may appear, however these velocities

relevant in realistic two-phase flows.

In dispersed flow it is assumed that the pressure

phases is the same. The water surface in hori

stratified flow can be wavy, therefore the interfacial

termp, is applied to describe pressure gradients:

x as a

pressure

(8)

tion in

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

interface heat transfer terms (Q,, Q,).

Terms due to the variable pipe cross-section.

Terms with wall friction (Ffwai, Fg,wai).

Term with volumetric forces (g cosO).

Sources from the points 1. and 2. are the so-called

relaxation source terms. They play crucial role in the case of

condensation-induced water hammer and therefore their

detailed description is given in the next section. Other,

non-relaxation source terms are terms with wall friction and

volumetric forces.

(9) Relaxation source terms

Relaxation source terms are inter-phase mass, momentum,

nodulus and energy exchange terms, which tend to establish thermal

and mechanical equilibrium between the phases.

Characteristic time scale of relaxation source terms can be

much shorter than the characteristic time scale of the acoustic

al and waves (terms are stiff and need special numerical treatment).

ons in Relaxation source terms are flow regime dependent.

it, mass Avery crude flow regime map (described below) has been

pressure, applied in the WAHA code, which is actually nothing more

tions of than search for the best fit of the macroscopic data, and is

ir steam open for improvement with the further comparison with the

experimental data. More detailed flow regime maps were

abandoned as they are developed for the steady-state flow

regimes. The accuracy of the existing more detailed flow

(10) regime maps in the area of fast transients comparing to our

crude flow regime map, is in our opinion not significantly

10) are higher, and does not justify their use in the WAHA code. The

veloped main goal of the WAHA correlations is to have correct

iture or correlations in the limit of high and low vapor volume

ties are fractions with their smooth transition into the single-phase

basis of flow, with possibility of their further tuning on the basis of

d saved the experiments. It is important to note that even the

and 350 "standard" single-phase wall friction correlations (RELAP,

CATHARE) turn out to be insufficient in the area of the fast

used to transients. The WAHA code offers an option of the "unsteady

wall friction model" that takes into account additional wall

friction due to the unsteadiness of the flow.

9 The WAHA code distinguishes two flow regimes (Fig. 3):

) 1 dispersed flow with stratification factor S= 0 and horizontally

11 stratified flow with S=1. There is also transition area

between both regimes with 0

~~
ure the divided into bubbly flow (a< 0.5), droplet flow (a > 0.95) and~~

Applied transitional bubbly-to-droplet flow.

utiuInal

relative

complexx

are not

of both

zontally

pressure

p, = Sa(1 a)(pf p,)gD (12)

where S presents the stratification factor.

Terms that don't include derivatives are source terms and

they are flow regime dependent. Source terms in Eqs. (1)-(6)

are:

Terms with inter-phase drag (C,).

Terms with inter-phase exchange of mass and energy:

vapor generation rate (Fg),

Dispersed flow

Horizontally stratified Transitional area S=

Transitional area

flow

ow 1 > S > > 0 95 Droplet flow

S = 1 0 95 > a > 05 Transitional flow

a< 0 5 Bubbly flow

0.5 Vcnt Vent

0 Vr

Fig. 3: WAHA flow regime map.

Flow is dispersed with S= 0 if at relative velocity v, is larger

than the critical velocity:

l"vr "c rt

J ia (1-a)1

vcrit= gD(pf pg) +

S gP 9 ) (13)

This expression is an approximation based on the

Kelvin-Helmholtz instability. This critical velocity is at the

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

same time maximum relative velocity, where two-fluid

model with applied interfacial pressure term and without

virtual mass term, is still hyperbolic. Flow is horizontally

stratified with stratification factor S = for vr < vcr,,rc/2 .

Flow is transitional between dispersed and horizontally

stratified, if v ..I ,/2< Ivv, Stratification factor S is

linearly interpolated between 0 and 1 in such case. This

approach is similar to the well-known Taitel-Dukler

correlation for transition from stratified to slug flow (Taitel,

Dukler 1976). From the standpoint of condensation-induced

water hammer modelling, this model represents a possible

area of the future work, however it was not changed in the

present study.

The most important set of correlations for the present

research are stratified flow correlations that are crucial for

accurate description of the initial phase of the transient

including the formation of the slug. Slug formation means

transition in different flow regime that requires different set

of correlations. Despite the fact that this work considers

modelling of 1D slug flow with two-fluid models, existing

slug flow correlations (see Lin, Hanratty 1986, or Issa,

Kempf, 2003 for example) are not directly applicable as they

give an averaged heat, mass and momentum transfer

correlations instead of instantaneous local values needed for

the present study. A single slug is being explicitly followed in

the present study with 1D two-fluid model. Such tasks are

expected to be performed with multidimensional CFD

analysis using some of the free surface tracking algorithms

and not in 1D two-fluid models. However, while the current

CFD codes might be able to describe the stratified flow with

condensation and also slug formation and development, one

can certainly expected problems with modelling of the thin

condensate film on the walls and water hammer shock waves

followed by the flashing of the liquid. And while we do

intend to test CFD models for condensation-induced water

hammer in the future, our present goal represents

development of the suitable 1D two-fluid model.

Inter-phase momentum transfer

The interfacial friction coefficient C, in momentum

equations is calculated from correlations, which are valid

for two-phase flow water-vapor and for two-component

flow water-ideal gas (similar to RELAP5 model). The

original WAHA correlations remain unchanged for the

present simulations and analyses have shown rather low

sensitivity of the results to the inter-phase friction

coefficients in stratified, dispersed and transitional flow.

Horizontally stratified flow interfacial friction coefficient

is calculated from the equation, which states that magnitude

of the drag force of the gas on the liquid is equal to the drag

force of the liquid on gas:

F,= F = C, (v -v ) (14)

interfacial friction coefficient is then calculated as:

1 (, v)2

( Pf, V)2

k =g,f

the vapor volume fraction a:

1. a < 0.5 (Bubbly flow):

C, pCaf

8

with drag coefficient of the slug:

C = 24(1 + 0.1Re 75)/Re

and interfacial area concentration:

ad = 3.6ab,, /do

where abb is modified vapor volume fraction, do is average

slug diameter and Re is Reynolds number.

2. a > 0.95 (Droplet flow):

C, =maxlp,gCag, 0. 1 (19)

(8

with the drag coefficient of the droplet:

C, =min(24(l+0.1Re 75)/Re, 0.5)

and the interfacial area concentration:

ag. = 3 .'I.,i. 10 4)/d0

where adp is modified liquid volume fraction and do is

average droplet diameter.

3. 0.5
inter-phase friction coefficient is calculated with

interpolation:

=(c- bubbly" (c droplet)( q) (22)

with exponent q:

( 0.95-a

q 0.95-0.5)

that was chosen to ensure smooth transition between

correlations in Eq. (16). and Eq. (19).

Dispersed-to-horizontally stratified interfacial friction

coefficient is calculated with interpolation:

C, =S(C,sied) + (1-S)(Cdpe) (24)

Inter-phase heat and mass transfer

Calculation of inter-phase heat and mass transfer were

significantly modified for the present condensation-induced

water hammer research. Original WAHA correlations do not

take into account wall heat transfer, which is an important

mechanism for the present work. are valid only for

water-vapor two-phase flow. The "standard" inter-phase

mass transfer (vapor generation rate F ) is calculated as:

where fk are friction factors, v, is interface velocity and

a g is interfacial area.

Dispersed flow coefficients are further divided according

F -f + ,g

g h -hf

9 f

where hk are specific enthalpies and Q,k are

Paper No

Paper No

liquid-to-interface and gas-to-interface heat fluxes. The

volumetric heat fluxes are calculated as:

Q,, = H,,(T -T, )

k fg

The heat transfer coefficients Hk depend on flow regime.

Beside the interphase heat and mass transfer,

condensation on the wall is taken into account as:

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

if slug head identified:

H,f = -Ca(1- a) Va f(T, Tg) v, c3

else

H, -0

enf =

endif

IfT Tf T I C4

f(Tf,ITg)= f T- "4

Tgr iS

Fg wal H,g-wao (Twall- T,)

F g wall : *

h* -h

9

T
;T > T

using Chato-Dobson correlation (Dobson, Chato, 1998) for

condensation rate in the horizontal pipe:

Nu = NU +( -(l f /rT)Nfod =

0.23Re'o2 Ga Pr 025

1+1.11x28 Ja

+(1- ,/;7T)0.0195Re8 Pr4qI4 fA(X,)

4k

H,g wall = Nulm D4 -

4k (28)

4k

H, =(1 -O /,r),Nfo~Cd D2

The first term in Chato-Dobson correlation represents part

of the coefficient, which represents contribution of the film

condensation on the wall and the second term represents

contribution of the forced-convective condensation on the

liquid-gas interface. The following variables are used in Eq.

(28): Rego Gas only Reynolds number, Ga Galilei number,

Ja Jacob number, <(/Xd) function of turbulent-turbulent

Lockhart-Martinelli parameter given in Dobson, Chato,

1998. Each term of Chato-Dobson correlation is used

separately: film condensation term in Eq. (27) and

forced-convection inter-phase exchange coefficient in Eq.

(26).

The vapor heat transfer coefficient H,g is calculated as

(similar in the RELAP5 code):

k

H,g = a, g 0.023Re8 (29)

D

Dispersed flow heat transfer is actually not present during

the transient the closest "approximation" of dispersed flow

can be seen at the head of the slug, where wave breaking

appears and causes much more efficient inter-facial heat

transfer than predicted by Chato-Dobson correlation. This

phenomena can be seen from the air-water experiments and

simulations of Bartosiewitz (2008) and experiment of

Vall6e et. al. (2010). Thus a very crude inter-phase heat and

mass transfer is used at the location of the slug head:

Values C1=1.5, C2=1, C3=0, C4=1, C5=1 were used in the

calculations collected in the Table 1 below. Head of the slug

is located with the gradient of the liquid superficial velocity:

if V((1- a)vv) < -0.05, slug head identified (31)

This is rather unusual approach for 1D two-fluid model and

actually represents some kind of inter-phase tracking within

the 1D two-fluid model. However, according to our

experience, this is the best way to perform

condensation-induced water hammer simulations with 1D

two-fluid model. If the slug, and especially the head of the

slug (where stratified heat and mass transfer correlations are

not applicable) is not successfully recognized and

condensation rate increased in that area, the simulations

exhibit rather poor results.

The vapor heat transfer coefficient H,g in dispersed

flow is calculated similar in the RELAP5 code, with a single

goal to almost instantaneously bring vapor to equilibrium:

H = (1+ 7.(100+ 25. 7)), (32)

where 77= max(-2,T, -7T) .

Unlike the standard WAHA code where transition

between horizontally stratified and dispersed flow

(stratification factor 0 < S <1) inter-phase heat transfer

coefficients are calculated with interpolation:

Hk =S(Hk-stried) +(1-S)(Hk-diped) (33)

the modified heat transfer coefficients are calculated as:

H, = max(H,k-satfied Hk-dspersed)

Inter-phase exchange correlations described with Eqs. (25)-

(34) are applied in WAHA code with minor correction

factors, that act at very low vapor or liquid volume fractions

and prevent negative values or extremely large values of

heat transfer coefficients.

Numerical Method

The system of six-equations model (Eqs. (1)

written in vectorial form:

- t &x

(6)) can be

where represents the non-conservative vector of the

independent variables:

vy = (p,a,f,vg,uf,u ) (36)

Paper No

further A and B are matrices of the system, and S is a

vector with non-differential terms in the equations.

Non-conservative schemes are known to converge to the

wrong solutions when shocks are present in the flow field,

however according to our experience (Tiselj, Petelin, 1997,

Tiselj et. al., 2004), non-conservative scheme does not seem

to be a big deficiency for short transients like water hammer

events.

The numerical scheme of the WAHA code is based on the

two step operator splitting and characteristic upwind method;

i.e. convection with non-relaxation source terms and

relaxation sources in Eq. (35) are treated separately:

A + B SN-R

t -x

A d -t

Each step of the operator splitting is performed with second

order accuracy. The applied operator splitting method is

formally first-order accurate. However, the numerical tests

showed, that despite the formally lower order of accuracy, the

results are practically the same as with the second-order

Strang operator splitting method.

Results and Discussion

Result of two simulations that were performed only with

stratified flow correlations for heat and mass transfer in Figs.

4. and 5. show a reasonable agreement between the pressure

measured in the experiments and simulation. First

experiment shown in Fig. 4 is case 03.01.01-13 with strong

pressure peak of 51 bar (see Table 1 for details of all

experimental runs and main results of the simulations). The

case 11.14.00-13 (Fig. 5) is experiment without recorded

pressure peak. Result in Fig. 4 shows that Chato-Dobson

correlation predicts reasonable condensation rate in the

horizontal test section, which results in a reasonable pressure

prediction. According to the calculations, slug is formed in

both cases due to the sufficiently large relative velocity

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

between the phases. However, since only stratified flow

correlations are used, the condensation rates ahead and

behind the slug are more or less the same and the slug is not

accelerated towards the closed end, but approaches the end

very slowly. Experiment 11.14.00-13 (Fig. 5) shows rather

smooth pressure growth, which means that liquid slug did not

develop in that experiment and that slug predicted by the 1D

two-fluid model does not exist. As is also shown later, the

existing 1D two-fluid model is not very accurate for the

experimental cases with small initial amount of the liquid.

It should be noted that static pressure measurement shown

in Fig. 4 does not include pressure peaks of the water hammer,

only some pressure fluctuations right behind and after the

pressure surge of the case 03.01.01-13 pressure are visible. It

is interesting to note that in the case 03.01.01-13 pressure

growth is stopped at time -1.4 s when the slug is formed and

stronger condensation at the head of the slug starts. However,

after the bubble behind the slug rapidly condenses and water

hammer peak is over (-1.9 s), the condensation rate remains

similar as in the Chato-Dobson based calculation.

2UUUUU

0 1 2 3 4 5 6

Fig. 4: Pressure (Pa) vs. time (s) at PACE3 static gas

pressure sensor: measurement (dashed) and calculation (solid

line) for test case 03.01.01.13.

Table 1: Overview of the selected experimental cases and corresponding simulations. Times of pressure peaks are given with

respect to the start of the gas injection. Times in Figures are given with respect to the starting time of the measurements.

experiment initial initial initial approximate hot gas measured time of calculated time of

pressure temperature vapor hot gas temperature pressure measured pressure calculated

in the in the test volume mass flow (K) peak pressure peak (bar) pressure

test section (K) fraction rate (kg/s) (bar) peak (s) peak (s)

section

(bar)

03.01.01-13 0.51 226.8 0.4736 0.31 293 51 0.58 56 0.58

03.01.01-14 0.55 228.6 0.4736 0.11 293 20 2.39 29 1.01

03.01.01-16 1.36 244.8 0.4736 0.32 295 27 0.63 17 0.92

05.21.01-12 0.45 224.8 0.4736 0.081 296 22 1.94 27 1.12

05.21.01-14 0.50 225.1 0.4736 0.15 296 44 0.81 43 0.83

11.13.00-31 0.52 227.8 0.4736 0.14 288 41 0.88 39 0.85

11.14.00-11 0.53 231.7 0.793 0.34 286 no shock / 47 0.76

11.14.00-12 0.57 232.1 0.793 0.42 287 25 0.98 48 0.66

11.14.00-13 0.65 233.9 0.883 0.43 287 no shock / 21 0.76

Paper No

350000

300000

250000

200000

150000

100000

50000 1-r---- ----

0 05 1 15 2 25 3 35 4 45 5

Fig. 5: Pressure (Pa) vs. time (s) at PACE3 static gas

pressure sensor: measurement (dashed) and calculation

(solid line) for test case 11.14.00-13.

Another effect is responsible for the differences in the

condensation rate at later times seen in Figs. 4. and 5.: pipe

wall, which is assumed to be at constant initial temperature

in the simulations, is in fact slowly warming up, reducing

the efficiency of the condensation. A conjugate heat transfer

model could be added in the two-fluid model to avoid that

effect, however, all the pressure peaks appear before the

times when these differences become relevant.

Similar agreement as seen in Fig. 4 and Fig. 5 exists for

other test cases collected in Table 1 and calculated with only

stratified flow heat and mass transfer correlations.

Water hammer modelling

An open issue remains upgrade of the two-fluid model

with a procedure that can recognize the head of the slug,

where inter-phase heat and mass transfer is much more

efficient than predicted by the Chato-Dobson correlation.

Slug head identification model from Eq. (31) was found to

predict the area of the slug head. Increase of the heat transfer

coefficients in the slug head region is performed with a

general model of Eq. (30). These two models with the

current values of the coefficients are being developed by

fitting of the calculations with the measurements and remain

open for further improvements.

Capabilities of the current form of both models gives

predictions of the condensation-induced water hammer

pressure peaks with accuracy shown in the Table 1, where

last 4 columns show measured and calculated magnitude

and time of the pressure peak. Good agreement of pressure

peak and timing is seen for cases 03.01.01-13, 05.21.01-14

and 11.13.00-31. Fig. 6 shows measured and calculated

pressure for the case 03.01.01-13, which can be considered

as a successful simulation. Secondary shock waves are seen

in computation and experiment. They are caused by a

classical "water column separation" mechanism, where

WAHA code is well tested and accurate. As shown in Fig. 7

that presents the same case, pressure peak causes only minor

changes in the integral condensation rate. Comparison of

Fig. 7 and Fig. 4. shows very similar calculated pressure

histories, despite the absence of the pressure peaks in Fig. 4.

Case 03.01.01-14 performed with slightly lower hot gas

mass flow rate and similar initial pressure and volume

fraction than cases 05.21.01-14 and 11.13.00-31 gives much

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

earlier time of water hammer than measured. Surprisingly,

comparing to the 03.01.01-14 slightly higher pressure peak

at earlier time was measured in the case 05.21.01-12, despite

even lower hot gas mass flow rate. This leads us to

conclusion, that tuning of the models cannot be performed

on one single experimental case due to the rather stochastic

nature of the whole phenomena. Thus, all the changes in the

models are continuously tested for all 9 test cases of the

Table 1.

Figs. 8 an 9 show an example of a simulation of modest

accuracy 05.21.01-12.

'11 14 00-testl 1strat out' u 1 5

/' FEB/EXPERIMENT/11 14 00/pace-pcb4-t11 exp u 1 4

/

Fig. 6: Measured (dashed) and calculated (solid) pressure

(Pa) vs. time (s) at the closed end for test 03.01.01-13.

500000

1 15 2 25 3 35 4 45 5

Fig. 7: Measured and calculated pressure(Pa) vs. time (s) in

the middle of the pipe (PACE3) for test 03.01.01-13.

The case 03.01.01-16 was performed at higher initial

temperature than other tests. Similar pressure peak is

measured and calculated, however, calculated pressure peak

occurs too late.

The worse results are obtained for high initial vapor

volume fractions. An example of poor simulation is given in

Figs. 10. and 11. for the case 11.14.00-11. Despite low

amount of liquid, 1D two-fluid model predicts formation of

the slug in all 3 cases 11.14.00-11, 12, 13. Slug formation is

followed by a strong water hammer, which is not seen in the

measurements at all, except a pressure peak of medium

magnitude in the case 11.14.00-12. As seen in Fig. 5 the

problem might not stem from the inter-phase heat and mass

transfer correlations but from the basic two-fluid equations

5e+06

4e+06

3e+06

2e+06

le+06

'03 01 01-test3 out' u 1 2

/ /FEB/ XPERIMENT/03 01 01/pace-pcb4-tl3dp exp' u 1 5 -----

1v

Paper No

and their capabilities to model stratified flows. Non-existing

slug in the case of Fig. 5 simulation is predicted even with

correlations for stratified flow. Pressure interface term that

makes the two-fluid model to behave like a shallow water

equations when stratification is assumed, is more accurate at

vapor volume fractions around 0.5, where circular pipe

behaves like a rectangular channel. At low (less than 0.2) or

high (higher than 0.8) vapor volume fractions the pressure

interface term might have a different form, which would

influence the dynamics of the large interfacial waves (slugs).

25e+06

2e+06

1 5e+06

le+06 F

500000

1 15 2 25 3 35 4 45

Fig. 8: Measured (dashed) and calculated (solid) pressure

(Pa) vs. time (s) at the closed end for test 05.21.01-12.

'05 21 01-test12NOV out u 1 5

140000 /FEB/EXPERIMENT/05 21 01/pace-pcb4-tl2dp exp'u 1 2---

120000

100000

80000

60000

40000

20000

0

1 15 2 25 3 35 4 45

Fig. 9: Measured and calculated pressure (Pa) vs. time (s) in

the middle of the pipe (PACE3) for test 05.21.01-12.

5e+06 ----------------------- ----

5e+ 6 11 14 00-test11 out u 1 2

/ /FEB/EXPERIMENT/11 14 00/pace-pcb4-t11 exp u 1 5 ----

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

400000

S11 14 00-testll out u 15 ---

/ I FEB/EXPERIMENT/11 14 00/pace-pcb4-t11 exp' u 1.--'-

350000 -

300000 -

250000

200000

150000

100000 -

50000 --- l -----------

05 1 15 2 25 3 35 4 45

Fig. 11: Measured and calculated pressure (Pa) vs. time (s) in

the middle of the pipe (PACE3) for test 11.14.00-11.

All computations were performed with an input model

consisting of -170 volumes with the horizontal test section

discretized in 60 volumes. Grid refinement was performed

for all cases in the Table 1 with test section discretized into

120 volumes. Pressure peaks and times of the peaks obtained

on the refined grid were typically up to 5% different.

Conclusions

Condensation-induced water hammer has been studied on

the experimental device described by Martin et. al. (2007)

and simulated with 1D two-fluid model of the computer

code WAHA (Tiselj et. al. 2004). For the purpose of the

present study WAHA was upgraded with ammonia

thermo-physical tables and Chato-Dobson correlation for

heat and mass transfer in stratified flow. This correlation

includes condensation of the ammonia gas on the cold pipe

wall, which is dominant in the first moments of the transient.

Several modifications of the transition criteria from

stratified to non-stratified flow are still being tested as well

as correlations for heat, mass and momentum transfer near

the head of the slug. Current models for condensation in the

slug head are being developed as a best fit with various

experimental runs. Results of the simulation show that 1D

two-fluid model can capture the main phenomena of the

condensation-induced water hammer, however, reliable

prediction of the condensation-induced water hammer in the

current configuration is still not possible. Behaviour at

various initial temperatures, pressures and hot gas flow rates

are well described for initial filling of the pipe of around

50%, while non-existent condensation-induced water

hammers are predicted at low liquid fillings (10-20%).

Acknowledgements

This research was financially supported by the Ministry of

Higher Education, Science and Technology, Republic of

Slovenia, project no. J2-1134 and research project of the EU

7th FP NURISP

05 1 15 2 25 3 35 4 45

Fig. 10: Measured (dashed) and calculated (solid) pressure

(Pa) vs. time (s) at the closed end for test 11.14.00-11.

"05 21 01-testl2NOV out u'l 2

/ /FEB/EXPERIMENT/05 21 01/pace-pcb4-tl2dp exp u 1 5 ----

S --.- ._- ,-

Paper No 7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

References 1978.

Bartosiewicz Y., Seynhaeve J.-M., Vall6e C., H6hne T., Vall6e C., Lucas D., Beyer M., Pietruske H., Schiitz P., Carl

Lavi6ville J.-M., Modeling free surface flows relevant to a H., Experimental CFD grade data for stratified two-phase

PTS scenario: comparison between experimental data and flows, Nuclear Engineering and Design, in press (2010)

three RANS based CFD-codes. Comments on the (doi:10.1016/j.nucengdes.2009.11.011)

CFD-experiment integration and best practice guideline,

proceedings of XCFD4NRS, Grenoble, France, 2008. Dobson M.K, Chato J.C, Condensation in smooth horizontal

tubes, Journal of Heat Transfer ASME 120, 193-213, 1998.

Bestion D., The Physical closure laws in the CATHARE

code, Nuclear Engineering and Design, 124, pp. 229-245,

1990.

Bjorge R.W., Griffith P, Initiation of Waterhammer in

horizontal and nearly horizontal pipes containing steam and

subcooled water, ASME Journal of Heat Transfer 106,

pp8350840, 1984.

Carlson K.E., Riemke R.A., Rouhani S.Z., Shumway R.W.,

Weaver W.L., RELAP5/MOD3 Code Manual, Vol. 1-7,

NUREG/CR-5535, EG&G Idaho, Idaho Falls, 1990.

Issa R.I., Kempf M.H.W., Simulation of slug flow in

horizontal and nearly horizontal pipes with the two-fluid

model, International Journal of Multiphase Flow 29, pp

69-96, 2003.

Lin PY., Hanratty T.J., Prediction of the initiation of slugs

with linear stability theory, International Journal of

Multiphase Flow 12(1), pp 79-98, 1986,

Martin C. S., R. Brown R., J. Brown J.,

Condensation-Induced Hydraulic Shock Laboratory Study,

ASHRAE 970-RP, June 2007

Martin C.S. Condensation-induced water hammer in a

horizontal refrigerant pipe, Proceedings of BHR Pressure

Surges,, 581-589, 2004.

Rukes B., Dooley R.B., Guideline on the IAPWS

Formulation 2001 for the Thermodynamic Properties of

Ammonia-Water Mixtures,

http://www.iapws.org/relguide/nh3h2o.pdf, (2001)

Strubelj L., Ezsol G, Tiselj I., Direct contact condensation

induced transition from stratified to slug flow. Nucl. Eng.

Des.. vol. 240, no. 2, str. 266-274, doi:

10.1016/j.nucengdes.2008.12.004, 2010

Taitel Y., Dukler A.E., A model for predicting flow regime

transitions in horizontal and near horizontal gas-liquid flows,

AICHE Journal 22, pp. 47-55, 1976.

Tiselj I., Horvat A. Cerne G, Gale J., Parzer I., Mavko B.,

Giot M., Seynhaeve J. M., Kucienska B., Lemonnier H.,

WAHA3 Code Manual, Jo2ef Stefan Institute Report,

IJS-DP-8841, March 2004

Tiselj I., Petelin S., Modelling of Two-Phase Flow with

Second-Order Accurate Scheme. Journal of Computational

Physics, vol. 136, pp. 503-521, 1997

Wylie E.B., Streeter VL., Fluid Transients, McGraw-Hill,