1.1 Fuel and core coolant channel temperatures

Before starting the experiments the thermal power released by the core was calibrated, according with the methodology developed by Mesquita et al. (2007). The calibration method used consisted of the steady-state energy balance of the primary cooling loop. For this balance, the inlet and outlet temperatures and the water flow in this primary cooling loop were measured. The heat transferred through the primary loop was added to the heat leakage from the reactor pool. The temperature measurements lines were calibrated as a whole, including sensors, cables, data acquisition cards and computer. The uncertainties for the temperature measurement circuit were ±0.4 oC for resistance temperature detectors, and ±1.0 oC for thermocouples circuits. The adjusted equations were added to the program of the data acquisition system (DAS). The sensor signs were sent to an amplifier and multiplexing board of the DAS, which also makes the temperature compensation for the thermocouples. The temperatures were monitored in real time on the DAS computer screen. All data were obtained as the average of 120 readings and were recorded together with their standard deviations. The system was developed to monitor and to register the operational parameters once a second in a historical database (Mesquita & Souza, 2010).

The original fuel element at the reactor core position B1 was removed and replaced by an instrumented fuel element. Position B1 is the hottest location in the core (largest thermal power production), according to the neutronic calculation (Dalle et al., 2002). The instrumented fuel element is in all aspects identical to standard fuel elements, except that it is equipped with three chromel-alumel thermocouples (K type), embedded in the fuel meat. The sensitive tips of the thermocouples are located along the fuel centreline. Their axial position is one at the half-height of the fuel meat and the other two 2.54 mm above and 2.54 mm below. Figure 6 shows the diagram of the instrumented fuel element and the Table I presents its main characteristics (Gulf General Atomic, 1972). Figure 7 shows the instrumented fuel element and one thermocouple inside a core channel.

979.0 mm

STAINLESS STEEL

TOP F TT NG

381.0 mm

LEAD-OUT TUB NG

REFLECTOR

SPACER

THERMOCOUPLES

THERMOCOUPLE

STANLESS STEEL

BOTTOM F TT NG

Fig. 6. Diagram of the instrumented fuel element

Parameter	Value
Heated length	38.1 cm
Outside diameter	3.76 cm
Active outside area	450.05 cm2
Fuel outside area (U-ZrHi.6)	434.49 cm2
Fuel element active volume	423.05 cm3
Fuel volume (U-ZrHr^)	394.30 cm3
Power (total of the core = 265 kW)	4.518 kW

Table 1. Instrumented fuel element features

Fig. 7. IPR-R1 core top view with the instrumented fuel element in ring B and one probe with thermocouple inside the core.

The instrumented fuel element was replaced to new positions and measures the fuel temperature in each one of the core fuel rings (from B to F). At the same way, two thermocouples were replaced to channels close to the instrumented fuel element to measure the coolant channel temperature. Experiments were carried out with the power changing from about 50 kW to 250 kW in 50 kW steps for each position of the instrumented fuel element. The fuel and coolant temperatures were monitored as function of the thermal power and position in the core.

In the TRIGA type reactors the buoyancy force induced by the density differential across the core maintains the water circulation through the core. Countering this buoyancy force are the pressure losses due to the contraction and expansion at the entrance and exit of the core as well as the acceleration and friction pressure losses in the flow channels. Direct measurement of the flow rate in a coolant channel is difficult because of the bulky size and low accuracy of flow meters. The flow rate through the channel may be determined indirectly from the heat balance across the channel using measurements of the water inlet and outlet temperatures. Two type K (chromel-alumel) thermocouples fixed in two rigid aluminum probes (7.9 mm of diameter), were inserted into the core in two channels close to position B1 (Channel 1 and 1′ in Fig. 4 and Fig. 5) and measured the inlet and outlet coolant channel temperatures. The probes penetrated axially the channels through small holes in the core upper grid plate. The probes were positioned in diametrically opposite channels, so that when a probe measured the channel entrance temperature, the other one registered the channel exit temperature. In a subsequent run, the probe positions were inverted. This procedure was used also for the Channels 1′, 2′, 3′, 4′ and 5′ (Fig. 5). There is no hole in the top grid plate in the direction of the Channel 0; so it was not possible to measure its temperature. The inlet and outlet temperatures in Channel 0 were considered as being the same of Channel 1. For the other channels there are holes in the top grid plate where it was possible to insert the temperature probes. To found the bulk coolant temperature axial profile at hot channel, with the reactor operating at 250 kW, the probe that measures the channel inlet temperature was raised in steps of 10 cm and the temperature was monitored. The same procedure was done with the reactor operating at 100 kW, but the probe was raised in steps of 5 cm.

In 1934 Fermi bombarded an uranium target with neutrons slowed down in paraffin in an attempt to produce transuranic elements. The first impression after the experiment was that uranium did undergo neutron capture and the reaction product was beta radioactive. Subsequent investigation of this reaction showed that the final activity produced included a range of different half-lives. This puzzle triggered intensive research from 1935 to 1939.

The identification of one of the activities produced as the rare-earth lantanum, first by Curie and Savitch in 1938 and then by Hahn and Strassmann in 1939, started to shed light on the puzzle. Indeed it was this fact that lead Hahn and Strassman to interpret the experimental activities as barium, lanthanum and cerium instead of radium, actinium and thorium. Shortly afterwards Meitner and Frisch (1939) suggested that the uranium nucleus, after the absorption of a neutron, splits itself into two nuclei of roughly equal size. Because the resemblance with the biological process in a living cell, the process was called fission. A typical example of the splitting is represented in Equation 1. Later measurements established the asymetric character of the process, the large energy release (~ 200 MeV) and the emission of prompt neutrons, which could trigger new fission processes and produce a chain reaction.

The first self-sustained chain reaction was achieved by Fermi in 1942 at the University of Chicago, which marked the begining of the nuclear age.

Since then many types of reactor have been developed for research, military and civil applications. Some examples include the Gas-Cooled reactor (GCR), Light Water Reactor (LWR), Heavy Water Reactor (HWR), Boiling Water Reactor (BWR), Liquid Metal Cooled Fast Breeder Reactor (LMCFBR) and so on. Independently of the kind of reactor one is considering there are important design and operating criteria which require a knowledge of the energy released in the decay of the fission products. In Table 1 an approximate distribution of the energy released in the fission of 235 U is presented. This table shows that from the total energy released inside the reactor (190 MeV), approximately 7 % is due to the beta decay of the fission products in the form of gamma and beta radiation. This source of energy is commonly called decay heat and depends on the fuel used in the reactor.

Distribution	MeV
Kinetic energy of light fission fragments	100
Kinetic energy of heavy fission fragments	67
Energy of prompt neutrons	5
Energy of prompt gamma rays	5
Beta energy of fission products	7
Gamma energy of fission products	6
Subtotal	190
Energy taken by the neutrinos	11
Total	201

Table 1. Approximate distribution of the energy released in the fission process of 235 U.

Experiments to understand the behavior of the nuclear reactors operational parameters allow improve model predictions, contributing to their safety. Developments and innovations used for research reactors can be later applied to larger power reactors. Their relatively low cost allows research reactors to provide an excellent testing ground for the reactors of tomorrow.

The experiments described here confirm the efficiency of natural convection in removing the heat produced in the reactor core by nuclear fission. The data taken during the experiments provides an excellent picture of the thermal performance of the IPR-R1 reactor core. The IPR-R1 TRIGA core design accommodates sufficient natural convective flow to maintain continuous flow of water throughout the core, which thereby avoids significant bubbles formation and restricts possible steam bubbles to the vicinity of the fuel element surface. The spacing between adjoining fuel elements was selected not only from neutronic considerations but also from thermohydrodynamic considerations. The experimental data also provides information, which allows the computation of other parameters, such as the fuel cladding heat transfer coefficient (Mesquita & Rezende, 2007). The theoretical temperatures and mass flux were determined under ideal conditions. There is a considerable coolant crossflow throughout the channels. Note that the natural convection flow is turbulent in all channels near the centre. The temperature measurements above the IPR-R1 core showed that water mixing occurs within the first few centimeters above the top of the core, resulting in an almost uniform water temperature. The temperature at the primary loop suction point at the pool bottom, as shown in Fig. 3, has been found as the lowest temperature in the reactor pool. Pool temperature depends on reactor power, as well as on the external temperature because it affects the heat dissipation rate in the cooling tower. The results can be considered as typical of pool-type research reactor. Further research could be done in the area of boiling heat transfer by using a simulated fuel element heated by electrical current (mock-up). The mock fuel element would eliminate the radiation hazard and allow further thermocouple instrumentation. By using a thermocouple near the fuel element surface, the surface temperature could be measured as a function of the heat flux.

It is suggested to repeat the experiments reported here, by placing a hollow cylinder over the core, with the same diameter of it, to verify the improvement of the mass flow rate by the chimney effect. These experiments can help the designers of the Brazilian research Multipurpose Reactor (RBM), which will be a pool reactor equipped with a chimney to improve the heat removal from the core (CDTN/CNEN, 2009).

Interfacial area concentration is defined as interfacial area per unit volume of two-phase flow. Therefore, the term "interfacial area concentration" is usually used in the meaning of

volume averaged value and denoted by ai. For example, one considers the interfacial area concentration in bubbly flow as shown in Fig.1. In this figure, Ai is instantaneous interfacial area included in volume, V. The volume averaged interfacial area concentration is given by

For simplicity, bubbles are sphere of which diameter is db, interfacial area concentration is given by

-V :rcNdb2 6a

i = V = db

Here, N is number of bubbles in volume V, and a is void fraction (volumetric fraction of bubbles in volume V).

Total Surface Area of Bubbles, Ai

Total Volume of Two-Phase Flow, V

Fig. 1. Interfacial area in bubbly flow

Similarly, time averaged interfacial area concentration, a i and statistical averaged interfacial area concentration, ai can be defined. The transport equation of interfacial area concentration is usually given in time averaged form in terms of time averaged interfacial area concentration, ai. However, for the derivation of the transport equation, it is desirable to formulate interfacial area concentration and its transport equation in local instant form.

Kataoka et al. (1986), Kataoka (1986) and Morel (2007) derived the local instant formulation of interfacial area concentration as follows.

One considers the one dimensional case where only one plane interface exists at the position of x=x0, as shown in Fig.2. In the control volume which encloses the interface in the width of Ax, as shown in Fig.2, average interfacial area concentration is given by

When one takes the limit of Ax ^ 0, local interfacial area concentration, ai, is obtained. It takes the value of zero at x Ф x0 and infinity at x=x0 . This local interfacial area concentration is given in term of delta function by

(4)

This formulation can be easily extended to three dimensional case. As Shown in Fig.3, three dimensional interface of gas and liquid is mathematically given by

f(x, y,z, t)=0

f(x, y,z, t)>0 (gas phase), f(x, y,z, t)<0 (liquid phase)

Fig. 3. Mathematical representation of three-dimensional interface

As shown in Fig.4, one considers the control volume which encloses the interface by following two surfaces.

f(x, y,z, t)=Af/2 (7)

f(x, y,z, t)=-Af/2 (8)

By the differentia! geometry, the width of the control volume is given by

Af/|grad f(x, y,z,,t)|

Then, average interfacial area concentration in this control volume is given by

ai = |grad f(x, y,z, t)| / Af

When one takes the limit of Af ^ 0, local interfacial area concentration, ai, is obtained by

ai=|grad f(x, y,z, t)| S(f(x, y,z, t)) (10)

where S(w) is the delta function which is defined by

where g(w) is an arbitrary continuous function..

In relation to local instant interfacial area concentration, characteristic function of each phase (denoted by (k) is defined by

(G =h(f(x, y,z, t)) (gas phase) (12a)

фL =1-h(f(x, y,z, t)) (liquid phase) (12b)

where suffixes G and L denote gas and liquid phase respectively. фк is the local instant void fraction of each phase and takes the value of unity when phase k exists and takes the value of zero when phase к doesn’t exist. Here, h(w) is Heaviside function which is defined by

h(w) =1 (w>0)

=0 (w<0) (13)

Heaviside function and the delta function are related by

Using above equations, the derivatives of characteristic function are related to interfacial area concentration as follows.

Eq.(15), local instant interfacial concentration is related to the derivative of characteristic function of each phase. Here, directional differentiation of characteristic function is considered. Spatial coordinate (x, y,z) is denoted by vector x and displacement vector of r is defined.

Fig. 5. Unit normal outward vector of phase k

x=(x, y,z)	(17)
r=(rx, ry, rz)	(18)

At position x, directional differentiation of characteristic function фі<(х) in r direction

d

(denoted by — ) is defined by dr

d

— <Mx) = nr • grad<Mx)

dr

= -nr • nkiai (19)

= — cos 9 ai

Here, nr is unit vector of r direction and 9 is the angle between nr and nki as shown in Fig. 6.

фк(х+г)=0 Фк(х+г)=1

In view of Fig.6, the product of фк(х+г) and Eq(19) is given by

(0<9<я /2) = — cos9 ai (я /2 < 9 < я)

Equation (19) is rewritten by

d 1

фк(х+r)— Фk(x) =—(-cos9 + |cos9 )ai dr 2

From Eqs.(19) and (21), one obtains

d d

— фk(x)-2фk(x + rb — фk(x) =-|cos9| ai dr dr

Integrating Eq.(22) for all r directions, one obtains

/Г J0{‘dr фk(x) — 2<^(x+r) dr фk(x)}sin 9d9dh =JO" jo-icos9 aisin mm^-2^

Rearranging Eq(23), one obtains

ai = — j^d” <k(x) — 2фк(х + r)d“ <k(x)}sin 9d9dh

d

As stated above, — is directional differentiation of characteristic function фк(х) dr

direction.

When one approximates the directional differentiation of characteristic function in Eq.(24) in the interval of | r |, one obtains,

— фк(х) *Mx±£bMx) dr r

<k(x+r) _d <k(x) .. фк(х+r)<k(x+r)-фк(х+r)<k(x) = фк(х+r) — фк(х+r)<k(x) dr |r| |r|

Then, the integrated function in Eq.(24) can be given by

— фк(х) — 2фк(х + r) — фк(х) ,фк(х){фк(х + r) -1} +фк(х + r){<Mx) -1} dr dr r

(„ |cos9| ,.фк(x){1 — фк(x + r)l + фк(x + r){1 — фк(x))^

V

Using this relation, Eq.(24) can be rewritten by

-=,1, f 2яг iim{ фк (х){! — фк(х+r)}+фк(x+r)(! — фк(x)} 1 2rcfo fo |r|—>o{ |r|

In the right hand side of Eq.(29), the term

фк(х){1 — фк(х + r)} + фк(х + r){1 — фк(х)}

represents the probability where gas-liquid interface exists between х and х+r as shown in Fig.7.

2. Basic transport equations of interfacial area concentration

Based on the rigorous formulation of interfacial area concentration, one can derive transport equation of interfacial area concentration. The transport equations of interfacial area concentration consist of two equations. One is the conservation equation of interfacial area concentration and the other is the conservation equation of interfacial velocity (velocity of interface), Vi.

Kataoba (2oo8) derived the local instant conservation equation of interfacial area concentration based on the formulation given by Eq.(24). In order to obtain the local instant conservation equation of interfacial area concentration, characteristic function of each phase

(denoted by фк) given by Eqs.(11) or (12) is needed. Kataoka(1986) also derived local instant formulation of two-phase flow which gives the local instant conservation equations of mass momentum and energy in each phase. The conservation equation of characteristic function of each phase (denoted by ф^ given by

Я

—(<kPk) + ^MkPk^ = — Pki(vki — Vi)• nkiai (k = G, L) (30)

ot

Here, pk, vk are density, velocity of each phase. Suffix ki is value of phase k at interface. Using Eqs.(24) and (30), Local instant conservation equation of interfacial area concentration is given by

Я-(ai) + grad(aiVi) = П —vkgrad(9k) — 2фk(x + r)0vkgrad(9k)}sinededq (31)

ot 2-J0 J0 or or

Averaging Eq.(31), one obtains the conservation equation of averaged interfacial area concentration by

—t(ai) + grad(aiVi) = 2- Ґ П——Гkgrad(фk) — 2фk(x + r)°ipgraii^J}sinededq (32)

ot 2—10 J0 or or

where averaged interfacial velocity Vi is defined by

V = (aiV)/a“ (33)

The right hand side of Eqs.(31) and(32) represent the source term of interfacial area concentration due to the deformation of interface. In the dispersed flows such as bubbly flow and droplet flow, this term correspond breakup or coalescence of bubbles and droplets.

Morel (2007) derived the conservation equation of averaged interfacial area concentration based on the detailed geometrical consideration of interface.

+ V^aiVi =ai(Vi • nGi )V • nGi (34)

Here, denotes time averaging and Vi is the time averaged velocity of interface which is given by

Vi = ai(Vi • nGi )nGi / ai (35)

The research group directed by Prof. Ishii in Purdue university derived the transport equation of interfacial area concentration of time averaged interfacial area concentration based on the transport equation of number density function of bubbles (Kocamustafaogullari and Ishii (1995), Hibiki and Ishii (2000a)) . It is given by

o“ 4

-a-+V • ai Vi = ^^

ot j = 1

Here, the first term in right hand side of Eq.(36) represent the source and sink terms due to bubble coalescence and break up. Interfacial area decreases when bubbles coalescence and increases when bubbles break up. This term is quite important in interfacial area transport. Therefore, the constitutive equations of this term are given by Hibiki and Ishii (2000a, 2000b) and Ishii and Kim (2004) based on detailed mechanistic modeling. The second term in right hand side of Eq.(36) represent the source and sink terms due to phase change. Equation (36) is practical transport equation of interfacial area concentration.

As for the conservation equation of interfacial velocity, Kataoka et al. (2010,2011a) have derived rigorous formulation based on the local instant formulation of interfacial area concentration and interfacial velocity, which is shown below. Since interface has no mass, momentum equation of interface cannot be formulated. Therefore, the conservation equation of interfacial velocity (or governing equation of interfacial velocity) has to be derived in collaboration with the momentum equation of each phase. Since interfacial velocity is only defined at interface, local instant formulation of interfacial velocity must be expressed in the form of

aiVi

Using Eq.(22), interfacial velocity is expressed by

aiVi|cos0| = — Viф^) + 2^фк(х + r)ck(x) (37)

or or

Considering Fig.7, interfacial velocity is approximated by velocity of each phase without phase change.

Vi ~<k(x)U-<k(x + r)}vk(x) + (k(x + r){l-cMx)}vk(x + r) (38)

From Eqs.(27) and (38) with some rearrangements, one obtains following approximate expression.

When one takes the limit of |r| ^ 0, one obtains

a V lcos0i lim ck(x){1 — ck(x + r)}vk(x) + ck(x + r){l -ck(x)}vk(x + r) (40)

i ^ |r| ^0 |r|

Integrating Eq.(40) in all direction, one finally obtains

aiVi = — l f2" f"limCk(x){1 -<Mx + r)}vk(x) + {1 — Mx)}<Mx + r)vk(x + r)si^edQdC (41)

i i 20 Jo|r|^0 r

Averaging Eq.(41), averaged interfacial velocity is given by

Viai = .! f2" f" lim (x){l — ck(x + r)}vk(x) + {1 — ck (x)}ck(x + r)vk(x + r) sinQdedC (42)

i i 20 Jo|r|^0 r T y ’

On the other hand, using Eq.(29), following relation can be obtained for averaged velocity of each phase, vk by

where vk is defined by

vk = №kvkV (44)

Using, Eqs.(42) and (43), the difference between time averaged interfacial velocity, Vi and time average velocity of phase k, vk is given by

(Vi — vk)ai

= 1 f2я Гlim ^k(x){1 -<Mx + r)lK(x)- vk(x)l +l1 — Фk(x)}Фk(x + r)K(x + r)- vk(x)l sin0ded, (45)

2„fo fo |r|—о |r| Ф

Rearranging the term in integration in the right hand side of Eq.(45) one obtains

lim Фk(x){1 — <Mx + r)}{vk(x) — vk(x)} +{1 — Фk(x)} <Mx + r){vk(x + r) — vk(x)} |r|—о |r|

} } } }

= — lim фk (x)фk(x + r)vk (x + r) + фk(x)фk (x + r)vk (x)

|r| —0 Irl

Here, фk and vk are fluctuating terms of local instant volume fraction and velocity of phase k which are given by

(47)

=ф — ф< (48)

Equations (45) and (46) indicate that the difference between time averaged interfacial velocity, Vi and time averaged velocity of phase k, vk is given in terms of correlations between fluctuating terms of local instant volume fraction and velocity of phase k which are related to turbulence terms of phase k.

Then, it is important to derive the governing equation of the correlation term given by Eq.(46). In what follows, one derives the governing equation based on the local instant basic equations of mass conservation and momentum conservation of phase k which are given below (Kataoka (1986)). In these conservation equations, tensor representation is used. Einstein abbreviation rule is also applied. When the same suffix appear, summation for that suffix is carried out except for the suffix k denoting gas and liquid phases.

(Mass conservation)

Svkp

*^-0 <49>

(Momentum conservation)

Averaging Eqs.(49) and (50) , one obtains time averaged conservation equation of mass and momentum conservation of phase k.

(Time averaged mass conservation)

(Time averaged momentum conservation)

Subtracting Eqs(51) and (52) from Eqs.(49) and (50), the conservation equations of fluctuating terms are obtained.

(Conservation equation of mass fluctuation)

я “7= vkpinkpiai

Sxp Фk

(Conservation equation of momentum fluctuation)

Using Eqs(53) and (54), one can derive conservation equation of

r r ft

(x)(Mx+r)vk(x+r)+<Mx)(k(x+r)vk(x)

Then, conservation equation of the difference between interfacial velocity and averaged velocity of each phase is derived. The result is given by

1 1 i*2" a" 1 ,,, , JPkr і JPk ■ «.пн = — "2 Jn Jn І1|П0іТ(^кФкг ^JkL + ФкФкг J-L)smЄаЄаФ pk 2%J0 J° |r|^o|L J-a J—a

"^Jo2* J0 йІП0|^(фкфкгр’kar +фкфкгР’ka)sin0аеаф 1 A" f2p

+vkaфk (фкг VkRг + фкг vkR )} sin еаеаф

J—R

As shown above, the formulation of governing equation of interfacial velocity is derived. Then, the most strict formulation of transport equations of interfacial area concentration is given by consenvation equation of interfacial area concentration (Eq.(32) , Eq(34), o! Eq.(36)) and TO^e^ation equation of interfacia1 velocity (Eq.(55)). As shown in Eq.(55), the conservation equation of interfacial velocity consists of various congelation terms of fluctuating terms of velocity and local instant volume fractions. These congelation temis represent the tuAulent transport of interfacial area, which reflects the interactions between gas liquid interface and tuAulence of gas and liquid phases. Equation (55) represents such tuAulence transport temis of interfacial area concentration. Accurate predictions of interfacial area transport can be possible by solving the transport equations derived here. Howeveq Eq.(55)
consists of complicated correlation terms of fluctuating terms of local instant volume fraction, velocity, pressure and shear stress. The detailed knowledge of these correlation terms is not available. Therefore, solving Eq.(55) together with basic equations of two-fluid model is difficult at present. More detailed analytical and experimental works on turbulence transport terms of interfacial area concentration are necessary for solving practically Eq.(55).

1.2 Experimental apparatus and procedure

The optical system of the Mach-Zehnder interferometer, MZC-60S is shown in figure 7 and figure 8 to visualize the exchange flow. After being rejoined behind the splitter, the test and reference laser beams interfere, and the pattern of interference fringes appears on the screen.

If the density of the test section is homogeneous, the interference fringes are parallel and equidistant (Keulegan, 1958). If it is not homogeneous, the interference fringes are curved. An inhomogeneity in the test section produces a certain disturbance of the non-flow fringe pattern. The digital camera and high-speed camera using D-file is able to attach to the interferometer.

The mass flow rate through the core coolant channels was determined indirectly from the heat balance across each channel using measurements of the water entrance and exit temperatures. Although the channels are laterally open, in this work cross flow or mass transfer between adjacent channels was ignored. Inlet and outlet coolant temperatures in channels were measured with two rigid aluminum probes with thermocouples. They were inserted in the upper grid plate holes (Fig. 5). Figure 8 illustrates schematically the general natural convection process established by the fuel elements bounding one flow channel in the core. The core coolant channels extend from the bottom grid plate to the top grid plate. The cooling water flows through the holes in the bottom grid plate, passes through the lower unheated region of the element, flows upwards through the active region, passes through the upper unheated region, and finally leaving the channel through the differential area between a triangular spacer block on the top of the fuel element and a round hole in the grid. As mentioned, in natural convection the driving force is supplied by the buoyancy of the heated water in the core channels.

In a typical TRIGA flow channel entire fuel element is cooled by single phase convection as long as the maximum wall temperature is kept below that required to initiate boiling. However, at higher power levels the inlet and outlet regions of the core, where the heat fluxes are the lowest, the channels are cooled by single phase convection. In the central region, where the axial heat flux is highest, the mode of heat transfer is predominantly subcooled boiling (Rao et al., 1988 and Mesquita et al. 2011).

The channel heating process is the result of the thermal fraction contributions of the perimeter of each fuel around the channel. So there was an average power of 4.518 kW dissipated in each stainless steel cladding fuel element and 4.176 kW dissipated in each aluminum cladding fuel element at 265 kW core total power. The values are multiplied by the fuel element axial power distribution and core radial power distribution factors as shown in profiles of Fig. 9.

Graphite

Lower Grid Plate _

(adapted from Veloso, 2005)

Fig. 8. A scheme of one flow channel in the TRIGA core

The power axial distribution factor in the fuel is 1.25, according with Marcum (2008). Figure 10 shows in detail the coolant channels geometry. The core radial power distribution factors, shown in Fig. 10, were calculated by Dalle et al. (2002) using WIMS-D4 and CITATION
codes. The products are multiplied by the fractions of the perimeters of each fuel in contact with the coolant in each channel. The two hottest channels in the core are Channel 0 and Channel 1′. Channel 0 is located closer to the core centre, where density of neutron flux is larger, but there is no hole in the top grid plate in the direction of this channel. Table 1 gives the geometric data of the coolant channels and the percentage of contribution relative to each fuel element to the channels power (Veloso, 2005 and Mesquita, 2005).

Fig. 10. Core coolant channels geometry and radial power distribution

Channel Number	Area [cm2]	Wetted Perimeter [cm]	Heated Perimeter [cm]	Hydraulic Diameter [cm]	Channel Power [%]
0	1.5740	5.9010	3.9060	1.0669	1.00
1′	8.2139	17.6427	15.1556	1.8623	3.70
2′	5.7786	11.7456	11.7456	1.9679	2.15
3′	5.7354	11.7181	11.7181	1.9578	1.83
4′	5.6938	11.7181	8.6005	1.9436	1.13
5′	3.9693	10.8678	3.1248	1.4609	0.35

Table 2. Channel geometry and hydraulic parameters (Veloso, 2005; Mesquita, 2005)

The mass flow rate in the hydraulic channel (m) in [kg/ s] is given indirectly from the thermal balance along the channel using measurements of the water inlet and outlet temperatures:

Where qc is the power supplied to the channel [kW], Cp is the isobaric specific heat of the water [J/kgK] and AT is the temperature difference along the channel [oC]. The mass flux G is given by: G = m / channel area. The velocity u is given by u = G / p, where p is the water density (995 kg/ m3). The values of the water thermodynamic properties were obtained as function of the bulk water temperature at the channel for the pressure 1.5 bar (Wagner & Kruse, 1988) Reynolds number (Re), used to characterize the flow regime, is given by:

Re = GDw (2)

Where G is the mass flux in [kg/m2s], Dw is the hydraulic diameter in [m] and ц is the dynamic viscosity [kg/ms].

Once the reactor is shutdown, the energy released in radioactive decay provides the main source of heating. Hence, the coolant needs to be maintained after the termination of the neutron-induced fission process in a reactor. The decay heat varies as a function of time after shutdown and can be determined theoretically from known nuclear data. Such computations are presently based on the inventory of nuclei created during the fission process and after reactor shutdown, and their radioactive decay characteristics:

f (t) = E(% + E7i + E«i AN (t) (2)

і

where f (t) is the power function, E; is the mean decay energy of the ith nuclide (в, Y and a components), A; is the decay constant of the ith nuclide (A; = ln(2)/T1/2;), and N;(t) is the number of nuclide i at cooling time t. In the summation calculations 2, the first step is the determination of the inventory of nuclei N;(t), which can be obtained by solving a linear system of coupled first order differential equations that describe the build up and decay of fission products:

(A; + фЩ + E fj^i AjNj + E ¥k^i ф Nk + yiF
j k

where Ni represents the number of nuclides i, Xi stands for the decay constant of nuclide i, Hi is the average capture cross section of nuclide i, ф is the neutron flux, fj^i is the branching ratio of the decay from nuclide j to i, Цк^і is the production rate of nuclide i per one neutron capture of nuclide k, yi is the independent fission yield of nuclide i and F is the fission rate. These calculations require extensive libraries of cross sections, fission yields and decay data.

As mentioned earlier, an accurate assessment of the decay heat is highly relevant to the design of nuclear facilities. Consider for example the safety analysis of a hypothetical loss-of-coolant accident (LOCA) in a light water reactor or the cooling needs of a spent-fuel pool. Calculations of the decay heat are also important for the design of the shielding of discharged fuel, the design and transport of fuel-storage flasks and the management of the resulting radioactive waste. This assessment is obviously not only relevant to safety, it also has economic and legislative consequences. For example, the accuracy of the presently available decay data is still not high enough and this situation translates into higher safety margins implying greater economic costs. Reducing the uncertainty in available decay data is one of the main objectives of the work devoted to decay heat in present-day research.

Nowadays the most extended way to calculate the decay heat in reactors is the summation calculation method based on equations 2 and 3. As can be seen from 2 this method is simply the sum of the activities of the fission products produced during the fission process and after the reactor shutdown weighted by the mean decay energies. As a consequence, the summation calculation method relies heavily on the available nuclear data and this is the reason why in the early days of the nuclear age this method was not satisfactory. A different method, the so-called statistical method was the prefered way of evaluation of the decay heat at that time.

The statistical method was based on the work of Way and Wigner (Way & Wigner, 1948). They considered fission products as a sort of statistical assembly and relying on mean nuclear properties, they deduced empirical relations for the radioactive half-lives and atomic masses of fission products. With these relations the gamma (Py MeV/fission-s) and beta plus gamma (Pt MeV/fission-s) decay heat power functions in 235 U were determined as follows:

PY (t) = 1.26t—12	(4)
Pt (t) = 2.66t-12	(5)

These relations were considered valid for decay times t in seconds in the range of 10 seconds to 100 days. For many years this was the only available method for calculations of the decay heat. Other fissile materials were supposed to behave similarly to 235 U. This method is clearly simpler than the summation calculations, however it is incomplete by nature and less accurate for longer cooling times. In a natural way, with an increasing volume of nuclear data, the statistical method was gradually superseded by the summation calculations. Even though it is not used anymore, the work of Way and Wigner (Way & Wigner, 1948) was a seminal article, and the interested reader is encouraged to read it. Their predictions are approximatelly off by 60 % from todays accepted values, but the form of the parametrization of the decay heat as function a x t-b is still used for benchmarks determination (benchmarks are obtained as a sum of functions of the form a x t-b that cover the full cooling period, which are adjusted to experimental data).

It should be noted that in the beta part of equations 4 , 5 as well as in 2 the neutrino energy is not included.

Salah El-Din El-Morshedy

Reactors Department, Nuclear Research Center, Atomic Energy Authority Egypt;

1. Introduction

Research reactors with power between 1 MW and 50 MW especially materials testing reactors (MTR), cooled and moderated by water at low pressures, are limited, from the thermal point of view, by the onset of flow instability phenomenon. The flow instability is characterized by a flow excursion, when the flow rate and the heat flux are relatively high; a small increase in heat flux in some cases causes a sudden large decrease in flow rate. The decrease in flow rate occurs in a non-recurrent manner leading to a burnout. The burnout heat flux occurring under unstable flow conditions is well below the burnout heat flux for the same channel under stable flow conditions. Therefore, for plate type fuel design purposes, the critical heat flux leads to the onset of the flow instability (OFI) may be more limiting than that of stable burnout. Besides, the phenomenon of two-phase flow instability is of interest in the design and operation of many industrial systems and equipments, such as steam generators, therefore, heat exchangers, thermo-siphons, boilers, refrigeration plants and some chemical processing systems. In particular, the investigation of flow instability is an important consideration in the design of nuclear reactors due to the possibility of flow excursion during postulated accident. OFI occurs when the slope of the channel demand pressure drop-flow rate curve becomes algebraically smaller than or equal to the slope of the loop supply pressure drop-flow rate curve. The typical demand pressure drop-flow rate curves for subcooled boiling of water are shown in Fig. 1 (IAEA-TECDOC-233, 1980). With channel power input S2, operation at point d is stable, while operation at point b is unstable since a slight decrease in flow rate will cause a spontaneous shift to point a. For a given system, there is a channel power input Sc (Fig. 1) such that the demand curve is tangent to the supply curve. The conditions at the tangent point c correspond to the threshold conditions for the flow excursive instability. At this point any slight increase in power input or decrease in flow rate will cause the operating point to spontaneously shift from point c to point a, and the flow rate drops abruptly from M to Mc. For MTR reactors using plate-type fuel, each channel is surrounded by many channels in parallel. The supply characteristic with respect to flow perturbations in a channel (say, the peak power channel) is essentially horizontal, and independent of the pump characteristics. Thus, the criterion of zero slope of the channel demand pressure drop-flow curve is a good approximation for assessing OFI, i. e.

Channel flowrate

Fig. 1. Typical S-curves to illustrate OFI, (IAEA-TECDOC-233, 1980)

Functionally, the channel pressure drop-flow curve depends on the channel geometry, inlet and exit resistances, flow direction, subcooled vapor void fraction, and heat flux distribution along the channel.

2. Background

There is a lot of research work in the literature related to flow instability phenomenon in two-phase flow systems. (Ledinegg, M., 1938) was the first successfully described the thermal-hydraulic instability phenomenon later named Ledinegg instability. It is the most common type of static oscillations and is associated with a sudden change in flow rate. (Whittle & Forgan, 1967) and (Dougherty et al., 1991) were performed an experimental investigations to obtain OFI data in a systematic methodology for various combination of operating conditions and geometrical considerations under subcooled flow boiling. (Saha et
al., 1976) and (Saha & Zuber, 1976) carried out an experimental and analytical analysis on the onset of thermally induced two-phase flow oscillations in uniformly heated boiling channels. (Mishima & Nishihara, 1985) performed an experiment with water flowing in round tube at atmospheric pressure to study the critical heat flux, CHF due to flow instability, they found that, unstable-flow CHF was remarkably lower than stable-flow CHF and the lower boundary of unstable-flow CHF corresponds to the annular-flow boundary or flooding CHF. (Chatoorgoon, 1986) developed a simple code, called SPORTS for two-phase stability studies in which a novel method of solution of the finite difference equations was devised and incorporated. (Duffey & Hughes, 1990) developed a theoretical model for predicting OFI in vertical up flow and down flow of a boiling fluid under constant pressure drop, their model was based on momentum and energy balance equations with an algebraic modeling of two-phase velocity-slip effects. (Lee & Bankoff, 1993) developed a mechanistic model to predict the OFI in transient sub-cooled flow boiling. The model is based upon the influence on vapor bubble departure of the single-phase temperature. The model was then employed in a transient analysis of OFI for vertical down-wards turbulent flow to predict whether onset of flow instability takes place. (Chang & Chapman, 1996) performed flow experiments and analysis to determine the flow instability condition in a single thin vertical rectangular flow channel which represents one of the Advanced Test Reactor’s (ATR) inner coolant channels between fuel plates. (Nair et al., 1996) carried out a stability analysis of a flow boiling two-phase low pressure and down flow relative to the occurrence of CHF, their results of analysis were useful in determining the region of stable operation for down flow in the Westinghouse Savannah River Site reactor and in avoiding the OFI and density wave oscillations. (Chang et al., 1996) derived a mechanistic CHF model and correlation for water based on flow excursion criterion and the simplified two-phase homogenous model. (Stelling et al., 1996) developed and evaluated a simple analytical model to predict OFI in vertical channels under down flow conditions, they found a parameter, the ratio between the surface heat flux and the heat flux required to achieve saturation at the channel exit for a given flow rate, is to be very accurate indicator of the minimum point velocity. (Kennedy et al., 2000) investigated experimentally OFI in uniformly heated micro channels with subcooled water flow using 22 cm tubular test sections, they generated demand curves and utilized for the specification of OFI points. (Babelli & Ishii, 2001) presented a procedure for predicting the OFI in down ward flows at low-pressure and low-flow conditions. (Hainoun & Schaffrath, 2001) developed a model permitting a description of the steam formation in the subcooled boiling regime and implemented it in ATHLET code to extend the code’s range of application to simulate the subcooled flow instability in research reactors. (Li et al., 2004) presented a three dimensional two-fluid model to investigate the static flow instability in subcooled boiling flow at low-pressure. (Dilla et al., 2006) incorporated a model for low — pressure subcooled boiling flow into the safety reactor code RELAP5/Mod 3.2 to enhance the performance of the reactor code to predict the occurrence of the Ledinegg instability in two-phase flows. (Khater et al., 2007a, 2007b) developed a predictive model for OFI in MTR reactors and applied the model on ETRR-2 for both steady and transient states. (Hamidouche et al., 2009) developed a simple model based on steady-state equations adjusted with drift-flux correlations to determine OFI in research reactor conditions; they used RELAP/Mod 3 to draw the pressure drop characteristic curves and to establish the conditions of Ledinegg instability in a uniformly heated channel subject to constant outlet
pressure. From the thermal-hydraulic point of view, the onset of significant void (OSV) leads to OFI phenomena and experimental evidence shows also that OSV is very close to OFI (Lee & Bankoff, 1993; Gehrke & Bankoff, 1993). Therefore, the prediction of OFI becomes the problem of predicting OSV. The first study that addressed the OSV issue was performed by (Griffith et al., 1958), they were the first to propose the idea that boiling in the channel could be divided into two distinct regions: a highly subcooled boiling region followed by a slightly subcooled region, they defined the OSV point as the location where the heat transfer coefficient was five times the single-phase heat transfer coefficient. A few years later, (Bowring, 1962) introduced the idea that OSV was related to the detachment of the bubbles from the heated surface and the beginning of the slightly subcooled region was fixed at the OSV point. (Saha & Zuber, 1974) developed an empirical model based on the argument that OSV occurs only when both thermal and hydrodynamic constraints are satisfied, where a general correlation is developed to determine OSV based on the Peclet and Stanton numbers. (Staub, 1968) postulated that OSV occurs when steam bubbles detach from the wall and assumed a simple force balance on a single bubble with buoyancy and wall shear stress acting on detach the bubble with surface tension force tending to hold it on the wall. He also postulated that the bubble could grow and detach only if the liquid temperature at the bubble tip was at least equal to the saturation temperature. (Unal, 1977) carried out a semi-empirical approach to determine and obtain a correlation of OSV point for subcooled water flow boiling. (Rogers et al., 1986; Chatoorgoon et al., 1992) developed a predictive model which relates the OSV to the location where the bubble first detaches assuming that bubble grow and collapse on the wall in the highly sub-cooled region. (Zeiton & Shoukri, 1996, 1997) used a high-speed video system to visualize the sub-cooled flow­boiling phenomenon to obtain a correlation for the mean bubble diameter as a function of the local subcooling, heat flux, and mass flux. (Qi Sun et al., 2003) performed a predictive model of the OSV for low flow sub-cooled boiling. The OSV established in their model meets both thermodynamic and hydrodynamic conditions. Several coefficients involved in the model were identified by Freon-12 experimental data.

It is clear that, there are several predictive models for OSV and OFI have been derived from theoretical and experimental analysis in the literature. However, their predictions in vertical thin rectangular channels still have relatively high deviation from the experimental data. Therefore, the objective of the present work is to develop a new empirical correlation with lower deviation from the experimental data in order to predict more accurately the OFI phenomenon as well as void fraction and pressure drop in MTR reactors under both steady and transient states.

As shown in the previous section, the rigorous formulation of transport equation of interfacial area concentration are given by conservation equation of interfacial area concentration (Eq.(32), Eq.(34) or Eq(36)) and conservation equation of interfacial velocity (Eq.(55)). However, Eq.(55) consists of complicated correlation terms of fluctuating terms of local instant volume fraction, velocity, pressure and shear stress. The detailed knowledge of these correlation terms is not available. Therefore, solving Eq.(55) together with basic equations of two-fluid model is difficult at present. More detailed analytical and experimental works on turbulence transport terms of interfacial area concentration are necessary for solving practically Eq.(55). From Eqs.(45) and (46), interfacial velocity is related to averaged velocity of phase k (gas phase or liquid phase)by following equation.

Viai = vkai -^-ї271 Ґ1 Um^cKkTkrVkr + <k<krv’k)sineded< (56)

When one considers bubbly flow and phase k is gas phase, Eq.(56) can be rewritten by

Viai = vGai -^-{0П |0 hm A(<G<Grv’ Gr + <G<Grv’ G)sin ededФ (57)

2 —J0 J0 |r|^0 r

From Eqs.(28) and (29) , following relation is derived.

ai = ai — ai = -1^1-!^ Ю |lrim0|10 {<k(1 — 2Фkr) + Фіт(1 — 2<J — 2(ФkФkr — Ф’kФ’kr)} sin ededФ (58)

In Eq.(57), the terms, <‘G/|r| and фGr /are related to the fluctuating term of interfacial area concentration. On the other hand, the terms, <Grv’ Gr and фGv’G are the fluctuating term of gas phase velocity at the location, x+r. and x. Therefore, the second term of right hand side of Eq.(57) is considered to correspond to turbulent transport term due to the turbulent velocity fluctuation. In analogous to the turbulent transport of momentum, energy (temperature) and mass, the correlation term described above is assumed to be proportional to the gradient of interfacial area concentration which is transported by turbulence (diffusion model). Then, one can assume following relation.

— 2-!Г!0 liPl^G^’ Gr +MGrv’ G)sin 6ded<H-Daigradai (59)

Here, the coefficient, Dai is considered to correspond to turbulent diffusion coefficient of interfacial area concentration. In analogy to the turbulent transport of momentum, energy (temperature) and mass, this coefficient is assumed to be given by

Dai « |vg|L (60)

Here, L is the length scale of turbulent mixing of gas liquid interface and |vG| is the turbulent velocity of gas phase. In bubbly flow, it is considered that turbulent mixing of gas liquid interface is proportional to bubble diameter, dB and the turbulent velocity of gas phase is proportional to the turbulent velocity of liquid phase. These assumptions were confirmed by experiment and analysis of turbulent diffusion of bubbles in bubbly flow (Kataoka and Serizawa (1991a)). Therefore, turbulent diffusion coefficient of interfacial area concentration is assumed by following equation.

Dai = K1KK = 6K1 =| v’l| (61)

ai

Here, a is the averaged void fraction and |v’L| is the turbulent velocity of liquid phase. K is

empirical coefficient. For the case of turbulent diffusion of bubble, experimental data were well predicted assuming K1=1/3 For the case of turbulent diffusion of interfacial area

concentration, there are no direct experimental data of turbulent diffusion. However, the diffusion of bubble is closely related to the diffusion of interfacial area (surface area of bubble). Therefore, as first approximation, the value of K1 for bubble diffusion can be applied to diffusion of interfacial area concentration in bubbly flow.

Equations (61) is based on the model of turbulent diffusion of interfacial area concentration. In this model, it is assumed that turbulence is isotropic. However, in the practical two-phase flow in the flow passages turbulence is not isotropic and averaged velocities and turbulent velocity have distribution in the radial direction of flow passage. In such non-isotropic turbulence, the correlation terms of turbulent fluctuation of velocity and interfacial area concentration given by Eq.(57) is largely dependent on anisotropy of turbulence field. Such non-isotropic turbulence is related to the various terms consisting of turbulent stress which appear in the right hand side of Eq.(55). Assuming that turbulent stress of gas phase is proportional to that of liquid phase and turbulence model in single phase flow, turbulent stress is given by

vLvL = — SLTpfVvL + t(^vL)} + |kSij (62)

Here, sLTP is the turbulent diffusivity of momentum in gas-liquid two-phase flow. For bubbly flow, this turbulent diffusivity is given by various researchers (Kataoka and Serizawa (1991b,1993)).

SLTP — 3adB I v’L (63)

Based on the model of turbulent stress in gas-liquid two-phase flow and Eq.(55), it is assumed that turbulent diffusion of interfacial area concentration due to non-isotropic turbulence is proportional to the velocity gradient of liquid phase. For the diffusion of bubble due to non-isotropic turbulence in bubbly flow in pipe, Kataoka and Serizawa (1991b,1993) proposed the following correlation based on the analysis of radial distributions of void fraction and bubble number density.

JB — K2adBn^^’L (64)

8y

Here, JB is the bubble flux in radial direction and nB is the number density of bubble. y is radial distance from wall of flow passage. K2 is empirical coefficient and experimental data were well predicted assuming K2=10. In analogous to Eq.(64), it is assumed that turbulent diffusion of interfacial area concentration due to non-isotropic turbulence is given by following equation.

Jai — K2adBarT-L (65)

8y

Here, Jai is the flux of interfacial area concentration in radial direction. Equation (64) can be interpreted as equation of bubble flux due to the lift force due to liquid velocity gradient.

As shown above, turbulent diffusion of interfacial area concentration due to non-isotropic turbulence is related to the gradient of averaged velocity of liquid phase and using analogy to the lift force of bubble, Eq.(58) can be rewritten in three dimensional form by

[Г [0 ,limn(^G^CrV’Gr +ФсФсУ G)sin0d0d(H

2nJ0 J0 |r|^°|r| (66)

-Daigradai + Caai(vG — VL) X rot(vL)

Empirical coefficient C in the right hand side of Eq. (66) should be determined based on the experimental data of spatial distribution of interfacial area concentration and averaged velocity of each phase. However, at present, there are not sufficient experimental data. Therefore, as first approximation, the value of coefficient C can be given by Eq.(67)

C — K2dB / uR (based on Eq.(65)) (67)

Using Eqs(55) and (66), transport equation of interfacial area concentration (Eq.(32),(34) or (36)) can be given by following equation for gas-liquid two-phase flow where gas phase is dispersed in liquid phase for bubbly flow.

Here, D/Dt denotes material derivative following the gas phase motion and turbulent diffusion coefficient of interfacial area concentration, Dai is given by Eq.(61). Coefficient of turbulent diffusion of interfacial area concentration due to non-isotropic turbulence, C is given by Eq.(67). The third term in the right hand side of Eq.(68) is source term of interfacial area concentration due to phase change and density change of gas phase due to pressure change. rG is the mass generation rate of gas phase per unit volume of two-phase flow due to evaporation. фот and ф^ are sink and source term due to bubble coalescence and break up

Similarly, the transport equation of interfacial area concentration for droplet flow is given by

) *rot(vG)} + ^k

" div(ai vL )

Here, D/Dt denotes material derivative following the liquid phase motion and фсо and фві< are sink and source terms due to droplet coalescence and break up and TL is the mass generation rate of liquid phase per unit volume of two-phase flow due to condensation. Here, turbulent diffusion coefficient of interfacial area concentration is approximated by turbulent diffusion coefficient of droplet (Cousins and Hewitt (1968)) as first approximation. The coefficient с for turbulent diffusion of interfacial area concentration due to non­isotropic turbulence (or lift force term) can be approximated by lift force coefficient of solid sphere as first approximation.

The research and development of source and sink terms in transport equation of interfacial area concentration have been carried out mainly for the bubbly flow based on detailed analysis and experiment of interfacial area concentration which is shown below.

Hibiki and Ishii(2000a,2002) developed the transport equation of interfacial area mentioned above and carried out detailed modeling of source and sink terms of interfacial area concentration. They assumed that the sink term of interfacial area concentration is mainly due to the coalescence of bubble. On the other hand, they assumed the source term is mainly contributed by the break up of bubble due to liquid phase turbulence. Based on detailed mechanistic modeling of bubble liquid interactions, they finally obtained the constitutive equations for sink and source terms of interfacial area transport.

The sink term of interfacial area concentration due to the coalescence of bubbles, фот is composed of number of collisions of bubbles per unit volume and the probability of coalescence at collision and given by

, , rca2s1/3

where the term —-tt-tS——-

volume and the term exp

collision. db, є and a are bubble diameter, turbulent dissipation and surface tension. amax is maximum permissible void fraction in bubbly flow and assumed to be 0.52. Гс and KC are empirical constants and following values are given

Гс =0.188, Kc =1.29 (71)

As for the source term due to the break up of bubble, it is assumed that bubble break up mainly occurs due to the collision between bubble and turbulence eddy of liquid phase. The constitutive equation is given based on the detailed mechanistic modeling of this phenomenon as

Гва(1 — а)є1/3

where the term 11/3/—————— Г represents the number of collisions of bubble and

db (amax — a)

turbulence eddy per unit volume and the term exp ^————— d s/з 2/3 J represents the

probability of break up at collision. ГB and Kb are empirical constants and following values are given.

Гв =0.264, Kb =1.37 (73)

The validity of transport equation of interfacial area concentration (Eq.(68)) and constitutive equations for sink term due to bubble coalescence (Eq.(70)) and source term due to bubble break up (Eq.(72)) are confirmed by experimental data as will be described later in details.

Hibiki and Ishii (2000b) further modified their model of interfacial area transport and applied to bubbly-to-slug flow transition. In bubbly-to-slug flow transition, bubbles are classified into two groups that are small spherical/distorted bubble (group I) and large cap/slug bubble (group II). They derived transport equations of interfacial area concentration and constitutive equations for sink and source terms for group I and group II bubbles based on the transport equation and constitutive equations for bubbly flow mentioned above.

Yao and More (2004) developed more practical transport equation of interfacial area concentration and constitutive correlations of source terms. They derived these equations based on the basic transport equation developed in CEA and models of source terms developed at Purdue University (Ishii ‘s group). They also developed sink term due to coalescence of bubbles which is given by

where g(a) and We is given by

We = 2PL(8db)2/3db

<j

On the other hand, source term due to break up of bubble by liquid phase turbulence is given by

The transport equation of interfacial area concentration and constitutive equations of source terms described above are implemented to CATHARE code which is developed at CEA using three-dimensional two-fluid model and k-є turbulence model. Predictions were carried out for thermal hydrodynamic structure of boiling and non-boiling (air-water) two- phase flow including of interfacial area concentration. Comparisons were made with experimental data of DEBORA experiment which is boiling experiment using R-12 carried out at CEA and DEDALE experiment which is air-water experiment carried out at EDF, Electricite de France. The predictions reasonably agreed with experimental data of boiling and non-boiling two-phase flow for distribution of void fraction, velocities of gas and liquid phase, turbulent velocity and interfacial area concentration and the validity of transport equation and constitutive equations described above was confirmed.

4.1 Experimental apparatus and procedure

The method of the mass increment was used for the investigations. Figure 7 shows a rough sketch of the apparatus. It consists of a test chamber, an electronic balance and a personal computer for data acquisition. The experimental procedure is mentioned in Sec. 2.1. Air enters the test chamber and the mass of the gas mixture in the test chamber increases.

The mass increment A m is automatically measured by the electronic balance with high accuracy. From mass increment data, the density increment of the gas mixture A pL = A m/V is calculated. The density increment means the difference of densities of the gas mixture from the density of pure helium in the test chamber. Then, volumetric exchange flow rate is evaluated by the following equation:

Q (2)

Ph Pl

The densimetric Froude number is defined by the following equation derived from the dimensional analysis suggested by Keulegan (Merzkirch, 1974):

Fr = Q I——- AjgDAp

In the above equations, V is the volume of test chamber, pH the density of air, pL the density of gas mixture in the test chamber, A pL ( =pH — pHe) = the density increment of the gas mixture, t the elapsed time, U(=Q/A) the exchange-velocity, p( =pH +pHe)/ 2, D the diameter and g the acceleration of gravity. The experiments are performed under atmospheric pressure and room temperature using the vertical and inclined round tubes, and using the vertical annular tube. The density of the gas mixture is close to that of helium in the present experiment. The sizes of the tubes are as follows. The diameter of the round tube D is 20 mm, which is much smaller than that of the test chamber. The inclination angle 0 ranges from 15 to 90 deg and the height L ranges from 0.5 to 200 mm.