A. Khedr: Thermal-Hydraulic Fortran Program for Steady-State Calculations of Plate-Type Fuel ...

19

THERMAL-HYDRAULIC FORTRAN PROGRAM FOR

STEADY-STATE CALCULATIONS OF PLATE-TYPE FUEL
RESEARCH REACTORS
by

Ahmed KHEDR

Received on February 26, 2008; accepted in revised form on June 11, 2008

The safety assessment of research and power reactors is a continuous process covering
their lifespan and requiring verified and validated codes. Power reactor codes all over
the world are well established and qualified against real measuring data and qualified
experimental facilities. These codes are usually sophisticated, require special skills and
consume a lot of running time. On the other hand, most research reactor codes still re-
quire much more data for validation and qualification. It is, therefore, of benefit to any
regulatory body to develop its own codes for the review and assessment of research re-
actors.
The present paper introduces a simple, one-dimensional Fortran program called
THDSN for steady-state thermal-hydraulic calculations of plate-type fuel research
reactors. Besides calculating the fuel and coolant temperature distributions and pres-
sure gradients in an average and hot channel, the program calculates the safety limits
and margins against the critical phenomena encountered in research reactors, such as
the onset of nucleate boiling, critical heat flux and flow instability. Well known ther-
mal-hydraulic correlations for calculating the safety parameters and several formulas
for the heat transfer coefficient have been used. The THDSN program was verified by
comparing its results for 2 and 10 MW benchmark reactors with those published in
IAEA publications and a good agreement was found. Also, the results of the program
are compared with those published for other programs, such as the PARET and
TERMIC.

Key words: research reactor, steady-state thermal-hydraulic calculation, safety parametars,

Fortran program

INTRODUCTION testing reactor (MTR) type fuel elements. They usu-

ally operate in the single-phase, liquid-water regime,
Research reactors (RR) are used all over the at low pressure and temperature. A safe operation of
world for various purposes, e. g. research, experi- these reactors requires some limits on the variables of
ments, education, training, radioisotope production, the major process in order to protect the reactor barri-
neutron radiography, material tests, etc. Most of these ers and prevent uncontrolled radioactivity releases.
RR are open-pool, light-water cooled and material These limits are established during the design stage
and should be verified at any reactor modification, up-
grading, core conversion, and so on.
Scientific paper The review and assessment process by the regu-
UDC: 621.039.572:519.876.5 latory body for nuclear facilities usually requires qual-
BIBLID: 1451-3994, 23 (2008), 1, pp. 19-30
DOI: 10.2298/NTRP0801019K ified codes covering the different fields of importance
to safety. The qualification of these codes depends on
Author's address:
National Center for Nuclear Safety and Radiation Control,
the hazard that may result in case of a transient or acci-
Atomic Energy Authority dent in the said facility. For example, huge amounts of
3, Ahmed El-Zomor Str., Nasr City, 11762 Cairo, financial and human resources have been invested into
P. O. Box 7551, Egypt
the improvement and development of qualified codes
Author's e-mail: for power reactors. On the other hand, available codes
ahmedkhedr111@yahoo.com for the review and assessment of research reactors still
20 Nuclear Technology & Radiation Protection –1/2008

require much more effort for their qualification. In ad- ant temperature is much lower than the saturation tem-
dition to this, the in-depth concept of defense requires perature, is the flow excursion or Ledinegg instability
that the safety analyses are independently assessed by [11].
the operating organization and by the regulatory body The proposed program is a Fortran-type program
[1]. Consequently, it is preferable for the regulatory called THDSN that calculates the TH and safety param-
body to have its own codes or to use internationally eters of plate-type RR during the build-up of
qualified codes, rather than those used by the designer steady-state conditions. In this program, the reactor
or the operating organization. core is represented by two channels, an average and a
The verification of RR design goals during nor- hot channel. An analytical solution of the energy equa-
mal operation is one of the regulatory body’s func- tion is used to calculate the temperature and power dis-
tions. Most international thermal-hydraulic (TH) tribution in the two channels. Well-known correlations
qualified codes applicable for RR safety analysis, such for core pressure drop, ONB, DNB, OFI, and critical
as PARET and the new versions of RELAP, are tran- velocity are also established in this program. There are
sient TH system codes [2, 3]. These codes are usually different correlations for the heat transfer coefficient
sophisticated and require skills during the preparation and an access for the user to choose between them. This
of their input deck, nodalization qualifications and program is simple, fast and provides an oversight of all
output processing [4, 5]. Although such packages are TH and safety parameters during the steady-state opera-
accurate and have capabilities to simulate a wide range tion of research reactors. For sake of verification, best
of TH transients, different research institutes have dif- estimate safety limits and margins for a typical 2 MW
ferent approaches to the build-up of simplest transient or 10 MW MTR IAEA reactor were calculated at
programs [6-9]. It is, therefore, not rational to use such steady-state conditions and the results are compared
large packages for steady-state calculations. The pres- with those published in IAEA TECDOC-233, Appen-
ent paper is focused on the build-up of a thermal-hy- dix A. The program was also verified against other TH
draulic Fortran program to be used in the safety codes, such as PARET and TERMIC.
assessment of plate-type fuel research reactors.
It is well known that the fuel clad represents the
most important barrier against the release of radioac- MODELING
tive materials [10]. To protect this barrier, from the
thermal point of view, some safety limits on the vari- In this section, the different correlations used in
ables of the said process have to be implemented. In the THDSN program to calculate the axial distribution
addition, safety margins against the appearance of any of heat flux, coolant temperature, clad and fuel tem-
critical/undesirable phenomena that may lead to a peratures and core pressure drop are presented. Also,
rapid increase in the cladding temperature and, conse- the correlations used to calculate the heat transfer co-
quently, damage in fuel elements are to be observed, efficient, the heat flux at ONB, the heat flux at DNB,
also. and the heat flux at OFI are introduced.
The most critical TH phenomena which, if ap-
proached too closely, will present a safety problem re-
garding fuel-plate integrity, are the departure from nu- Heat flux
cleate boiling (DNB) and the onset of flow instability
(OFI). In addition to this, coolant velocity beyond The core axial power distribution is considered
which the fuel plates will collapse (critical velocity) to follow a cosine shape given by
and the core pressure drop, are among the safety pa- pz
rameters that are of utmost importance. Although the q¢ ( i, z ) = q¢c ( i )cos (1)
le
onset of nucleate boiling (ONB) does not correspond
where q'(i, z) is the power density of fuel at axial loca-
to any critical phenomenon, for design purposes, it is
tion z in channel i (i equals 1 for an average channel
considered as a conservative constraint.
and 2 for a hot channel), and le is the extrapolated
Flow instability (FI) is a critical phenomenon be-
length. q¢c ( i ) is the maximum power density in chan-
cause the associated flow oscillations affect the local
nel i at half of the fuel plate length, where the axial co-
heat transfer characteristics and may induce premature
ordinate (z) is put equal to zero and is given by
burnout. Burnout heat flux occurs under unstable flow
conditions well below the burnout heat flux under sta- q¢c ( i ) = F ( i ) q¢a (2)
ble flow conditions (DNB). For practical purposes, the
heat flux that leads to the OFI may be more limiting in where q¢a is the core average power density which
the design of plate-type fuel elements than that of a sta- equals core power divided by the fuel meat volume. In
ble burnout [11]. In addition, there are different types an average channel F(i) equals the axial power factor
of FI encountered in heated, forced convection chan- (FA), in a hot channel F(i) equals the nuclear power
nels. The most common one encountered in RR, where factor (FNUC) which is the multiplication of the axial
the operating system pressure is low and the inlet cool- power factor FA and the radial power factor FR.
A. Khedr: Thermal-Hydraulic Fortran Program for Steady-State Calculations of Plate-Type Fuel ... 21

The heat flux axial distribution in channel i as a where FD is a factor less than unity, WT – the total core
function of the local power density is given by flow rate, and FPTN – the total number of fuel plates.
Water density r is evaluated at the core inlet tem-
q¢ ( i, z )
q( i, z ) = (3) perature.
SVR Channel velocity Vch is
where SVR is the fuel meat surface to volume ratio.
mch
Vch = 104 (8)
rA ch
Temperature distribution where Ach is the channel cross-sectional area. Water
density r is evaluated at the channel average tempera-
At steady-state, the heat generated in the fuel is ture. In some cases, the channel velocity instead of the
transferred to the coolant through the clad. The coolant core flow is given as input data. In such cases, the
axial temperature distribution in channel i is given by channel mass flow rate is calculated from eq. (8) and
[12] the volume flow rate is calculated from eq. (7).
T f ( i, z ) = T f1 +
Pressure drop
q¢ ( i ) A m le æ pz pl ö
+ 0.001 c ç ± sin + sin h ÷÷ (4)
pC p mch çè le 2le ø Most plate-type research reactors were designed
for single-phase flow under normal operation [11].
where Tf1 is the core coolant inlet temperature, Am – Core pressure losses are divided into two parts: pres-
the fuel meat cross-section area (= wmtm), wm – the sure losses in fuel element channels and pressure
meat width, tm – its thickness, mch – the channel mass losses in the fuel element nozzle. The channel pressure
flow rate, lh – the active fuel-plate length or heated losses (DPch) are: entrance losses DPen, friction losses
length of the channel, and Cp – the coolant specific DPf, and exit losses DPex i. e.
heat. The sign ± is for upward or downward core flow,
respectively. DPch = DPen + DPf + DPex (9)
The clad surface axial temperature distribution
(Tc) in channel i is given by [12] These three parts are given in bar by [11]

q¢c ( i )t m pz
Tc ( i, z ) = T f ( i, z ) + cos (5) rVch2
2h( i ) le DPen = 10–5 K
2
4 fL p rVch2
where h is the heat transfer coefficient. DPf = 10–5
The meat centerline axial temperature distribu- Dy 2
tion in channel i is [12] 2
–5 æV ö rVch2
DPex = 10 çç 1– o ÷
÷
T m ( i, z ) = Tc ( i, z ) + è Vch ø 2

where K is the entrance losses coefficient considered

pz æ t 2 t t ö to be 0.5, Lp – the total length of the fuel plate, Vo – the
+100.0 q¢c ( i )cos ç m + m c ÷ (6) water velocity in fuel element plenums, and Dy – the
le è 8k m 2k c ÷ø
ç
channel hydraulic diameter
where tc is the clad thickness, km – the meat thermal 4 × cross-section area
Dy =
conductivity, and kc – the clad thermal conductivity. weted perimeter

The nozzle pressure losses (DPnoz) consist of two

Mass flow rate parts: entrance losses (DPen) and friction losses (DPf) i. e.

The total core coolant flow is divided into two DPnoz = DPen + DPf (10)
parts. Active flow passes through the fuel elements
and non-active flow passes through the by-pass chan- These two parts are given in bar by
nels, such as the control rod channels. Active core flow
is determined by multiplying the total core flow in a 2
4 fL noz æ rV noz ö
factor less than one; in this study, this factor is 0.9. DPf = 10–5 ç ÷
Therefore, the mch is equal to Dy ç 2 ÷
è ø
FD WT r 2
æ rV noz ö
mch = (7) DPen = 10–5 K ç ÷
3600 FPTN ç 2 ÷
è ø
22 Nuclear Technology & Radiation Protection –1/2008

where Vnoz is the coolant velocity in the fuel element more conservative calculations, some references such
nozzle, K – a constant equal to one, and Lnoz – the noz- as TECDOC-233 [11] used correlation (15) for calcu-
zle length. lating the heat transfer coefficient into the core chan-
The friction factor f for turbulent flow in smooth nels.
channels is calculated as follows [13]:
– Sieder-Tate correlation [12]
f = 0.316 Re -0. 25 . × 105
Re < 10
0.14
(11) æm ö
1 Nu = 0.023 Re 0. 8 Pr 0. 4 çç w ÷÷ (16)
= 2 log 10 (Re f ) - 0.8 . × 105
Re > 10 è m ø
f
The total core pressure drop (DPcor) is given by All properties are evaluated at bulk temperature,
except mw which is evaluated at surface temperature.
DPcor = DPch + DPnoz ± DPst (12)
– Colburn correlation [14]
The plus sign in the downward flow and the neg-
Nu = 0.023 Re 0. 8 Pr 0. 33 (17)
ative sign in the upward flow. DPst is the static pressure
difference between the core inlet and outlet, given by
All properties are evaluated at film temperature
-7 (arithmetic mean between the fluid bulk and surface
DPst = 10 g ( L p + L noz ) (13)
temperatures), except for the specific heat which is
evaluated at fluid bulk temperature.
where g [Nm–3] is the water specific weight .
– Petukhov correlation [15]
f
Heat transfer coefficient Re Pr 0.11
8 æ mb ö
Nu = ç
çm
÷
÷ (18)
f 3 2
The program contains different correlations for 107
. + 12.7 ( Pr - 1) è w ø
calculating the heat transfer coefficient (h) and the 8
user has options to choose any of them. The flow in the
RR at normal operation is usually a turbulent flow; where
1
consequently, the correlations mentioned here cover f =
(182 . )2
. log 10 Re- 164
this flow regime. These correlations are summarized
as following: All properties are evaluated at film temperature
– Dittus-Boelter correlation [14] (arithmetic mean between the fluid bulk and surface
temperatures), except for mb and mw which are evalu-
There are two forms of this correlation. The first ated at bulk and wall temperature, respectively.
one is the usual formula which is used when the struc-
ture is cooled by a fluid, as in reactor core cooling, and
takes the form Onset of nucleate boiling
Nu = 0.023 Re 0. 8 Pr 0. 4 (14)
The ONB is not a limiting criterion in the design of
The second one is used when the structure is a fuel element, but is considered as a conservative state-
heated by a fluid, as in a heat exchanger primary side, ment for design purposes. The heat flux that initiates
and takes the form ONB is frequently used as a thermal design constraint.
Nu = 0.023 Re 0. 8 Pr 0. 3 (15) This constraint is determined by equating the clad sur-
where face temperature at the ONB which is calculated from the
Bergles and Rohsenow correlation with the clad surface
hD h G D mC p temperature calculated at the hot channel’s worst condi-
Nu = , Re = ch h , Pr = 103 , ¢¢ ) at the ONB, is
tion [11]. Therefore, the heat flux (qONB
k 100m k
calculated from the following correlation

k is the water thermal conductivity [Wm–1°C–1), m – the

water viscosity [Pa×s], and Gch – the channel mass flux P 0 . 0234
5 æ 9.23FA FR q av ö 2.16
(= rVch). Tsat + ç ÷ =
All the water properties in these two correlations 9è P 1.156 ø
are evaluated at bulk temperature. It is clear that the
heat transfer coefficient predicted by correlation (15) 20FR q av l h w h FA FR q av
= T f1 + + (19)
is lower than that predicted by correlation (14). So, for Gch t w C p w h
A. Khedr: Thermal-Hydraulic Fortran Program for Steady-State Calculations of Plate-Type Fuel ... 23

where wh is the channel heated width, w – the channel Mirshak and Labuntsov correlations must be
width, P – the absolute pressure, Tf1 – the core inlet used at positive sub-cooling. When the exit sub-cool-
temperature, and qav – an average heat flux. ing becomes negative, the burnout heat fluxes can be
Iterations are made on qav, starting from the aver- reasonably estimated by using these two correlations
age heat flux at nominal power (qa). The value of qav at extrapolated with DTsub equal to zero. The lower limit
which the two sides are equal is the average heat flux at for this extrapolation is the Rohsenow and Griffith
which the onset starts (qONB). The ratio between qONB pool-boiling critical heat flux correlation given by [11]
and qa is denoted by ONBR. This ratio is usually used
0. 6
to define the margin against ONB and its value must be -3 æ r - rv ö
q c = 121
. × 10 r v lçç l ÷
÷ (23)
greater than unity. The value of qONB results from eq.
è rv ø
(19) is lower than the actual value because it is calcu-
lated at the channel’s worst conditions. Consequently, where rl and rv [kgm–3], are the liquid and steam den-
its value will be multiplied by 1.15 [11]. sities, respectively.

Departure from nucleate boiling Onset of flow instability

Two correlations are recommended for the calcu- Flow instabilities are undesirable in heated
lation of burnout heat flux (critical heat flux) at research channels because flow oscillations affect the local heat
reactor conditions of low pressure and exit sub-cooling. transfer characteristics and may induce premature
These correlations are Mirshak and Labuntsov correla- burnout [11]. The burnout heat flux occurring under
tions [11]. These two correlations are: unstable flow conditions was well below the burnout
heat flux for the same channel under stable flow condi-
– Mirshak correlation tions. Consequently, the critical heat flux leading to
the onset of flow instability may be more limiting than
q c = 1510
. (1 + 01198
. Vch ) ×
that of stable burnout.
×(1+ 0.00914 DTsub )(1+ 019
. P) (20) The most common flow instabilities encoun-
tered in heated channels under nearly atmospheric
condition with forced convection are flow excursions
– Labuntsov correlation
of the Ledinegg-type. Experiments by Whittle and
0. 25
é 2.5Vch2 ù Forgan for sub-cooled water flowing (upward and
q c = 145.4 q( P ) ê1 + ú × downward) in narrow heated channels lead to their es-
êë q( P ) úû timate of flow instability occurring at an average heat
flux of [11]
æ 151. C p DTsub ö
×çç 1+ ÷
÷ (21) w t ch
è lP 0. 5 ø q FI = RrC p Vch (Tsat - T f1 ) (24)
wh lh
where
where
q( P ) = 0.99531P 1/ 3 (1- P / Pcr ) 4/ 3 , 1
R= ,
–1 D he
Pcr – the water critical pressure, and l [kJkg ] – the 1+ h
water latent heat. lh
The water exit sub-cooling (DTsub) is given by tch – the channel thickness, Dhe – the heated equivalent
diameter of the channel given by
20l h w h q c 2t ch w
DTsub = Tsat - T f1 + (22) channel flow area
rC p t w WAVch D he = 4 = ,
channel heated perimeter t ch + w h

The burnout, or critical heat flux, is determined Tsat – the water saturation temperature, and h – (ex-
by making iterations on the value of qc in eqs. (21) and perimental fit parameter) equal to 25.
(22), starting from the maximum heat flux in the hot Another correlation derived by Winkler for the
channel. When the difference between two subsequent average heat flux at onset of flow instability [11]
iterations is zero, the value of qc is equal to the critical
heat flux. This value should be greater than the maxi- q FI = -29.35 + (12815 . T f1 )Vch0. 8
. - 1104 (25)
mum heat flux in the hot channel and the ratio between
them represents the margin against DNB, denoted by
DNBR. DNB is a critical phenomenon and must be The peak heat flux qFI is obtained by multiplying
prevented to preserve the fuel element. q FI by the axial power factor (FA).
24 Nuclear Technology & Radiation Protection –1/2008

Flow instability is a critical phenomenon and ity. All the data required for these calculations are put
must be prevented by adequate TH design. The ratio into an input unit and the results are saved in an output
between qFI and the maximum heat flux in hot channel unit. The program is simple and easily applicable for
qc represents the margin against FI. This ratio must be calculating TH parameters and safety limits and mar-
greater than unity. gins during normal reactor operation.

Critical velocity PROGRAM VERIFICATION

For a given plate assembly there is a critical flow 2 MW and 10 MW MTR IAEA reactors
velocity at which the plates become unstable and large
deflections of the plates can occur. A critical velocity The 2 MW and 10 MW MTR reactors described
formula derived by Miller [11] is in the IAEA document TECDOC-233, Appendix-A
0. 5
[11] are typical plate-type RR present all over the
5 3 3
2 é15× 10 E ( t p - t m ) t ch ù world. All the data regarding these two types of reac-
Vcrit = ê ú (26) tors, such as core dimensions, neutron information,
3ê
ë rw 4 (1- v 2 ) úû
operating conditions, safety limits, margins and so on,
that may be required for the application and verifica-
where E [bar] is the Young’s modulus of elasticity, tp – tion of new programs are available and documented in
fuel plate thickness, and n – Poisson’s ratio. that reference. Consequently, these two reactors are
For design purposes, it is recommended that the used for our own program verification. Table 1 shows
coolant velocity be limited to 2/3 of the critical velocity. the IAEA reactors’ data required for this verification
[11]. Water properties such as polynomial correlations
in temperatures required for the calculations are in-
PROGRAM STRUCTURE
cluded in the THDSN program. These correlations are
derived from a curve fitted to tabulated water proper-
The THDSN is a Fortran program built to oper-
ties at an atmospheric pressure present in ref. [13].
ate on PCs. The program simulates the reactor core as
two parallel channels, an average and a hot channel.
The hot channel represents a single channel and the av- Safety parameters
erage one represents all other core channels. The
THDSN consists of a main part and a number of sub- The THDSN program results regarding safety pa-
routines, each one calculating a certain parameter, as rameters of benchmark reactors in comparison with
illustrated in the program structure shown in fig. 1. As those mentioned in TECDOC-233 are shown in tabs. 2
shown, there are subroutines to calculate the tempera- and 3. In our calculations, the channel flow is consid-
ture distribution in the coolant, clad, and fuel meat. ered to be downward and its exit pressure equal to the
Other subroutines for calculating the heat transfer co- core exit pressure given in tab. 1. Table 2 shows that, ex-
efficient, core pressure drop and power and heat flux cept for the average heat flux and the flux at the ONB,
distribution are included. There are also subroutines there is a good agreement between the other calculated
for calculating the ONB, FI, DNB, and critical veloc- and referenced safety parameters. The deviation in the
average heat flux returns to the heat transfer area used in
the calculation. In the present study, contrary to ref. [11]
that used the plate surface area, the meat surface area is
the one used. This result can be checked by hand calcu-
lation of the average heat flux. With respect to ONB,
there is no clear reason for the over-estimation appear-
ing in the present calculated value. But it is important to
mention that in the present study, the water properties
and the heat transfer coefficient are recalculated at the
new, updated temperature during the iterations on the
average heat flux using eq. (22). In ref. [11], the water
properties and, consequently, the heat transfer coeffi-
cient, were considered constant.

Safety margins

The safety margin is equal to the value of the

safety parameter at the critical phenomenon divided
Figure 1. THDSN program structure by the corresponding nominal value. Therefore, tab. 3
A. Khedr: Thermal-Hydraulic Fortran Program for Steady-State Calculations of Plate-Type Fuel ... 25

Table 1. Input values used in the TH calculations for 2 MW and 10 MW MTR IAEA reactors [11]
Property 2 MW 10 MW Property 2 MW 10 MW
No. of standard fuel elements (SFE) 19 23 No. of plates in SFE 19 23
No. of control fuel elements (CFE) 4 5 No. of plates in CFE 15 17
Channel flow width [cm] 6.64 6.64 Plate total length [cm] 62.6 62.5
Channel heated width [cm] 6.30 6.3 Channel thickness [cm] 0.2916 0.219
Fuel plate heated length [cm] 60.0 60.0 Clad thickness [cm] 0.0381 0.038
Core inlet temperature [°C] 38.0 38.0 Meat thickness [cm] 0.051 0.051
Core exit pressure (bar absolute) 1.961 1.566 Radial power factor 2.0 1.78
Clad thermal conductivity [W/m°C] 180.0 180.0 Axial power factor 1.58 1.4
Meat thermal conductivity [W/m°C] 53.6 53.6 Total core flow [m3/h] 300 1000
Fuel element nozzle diameter [cm] 6.0 6.0 Channel velocity [m/s] 0.94 2.97
Fuel element nozzle length [cm] 18.0 18.0

Table 2. Best estimate TH safety parameters

Pressude drop across Average heat flux Average heat flux Burnout heat flux [W/cm2] Heat flux at OFI
Reactor channel [bar] [W/cm2] at ONB [W/cm2] [W/cm2]
Labuntsov Mirshak
power
[MW] Present Present Present Present Present Present
Ref. [5] Ref. [5] Ref. [5] Ref. [5] Ref. [5] Ref. [5]
study study study study study study
2 0.019 0.0186 6.28 5.8 13.58 11.4 231.4 231 230 231 102.7 102.2
10 0.199 0.193 21.54 20.54 41.97 35.9 353.2 353 265.6 266 208.0 208.8

Table 3. Best estimate TH safety margins

Margin to ONB Margin to DNB, eq. (26) Margin to OFI Critical velocity
Reactor eq. (23) eq. (30) [m/s]
Labuntsov Mirshak
power
[MW] Present Present Present Present Present
Ref. [5] Ref. [5] Ref. [5] Ref. [5] Ref. [5]
study study study study study
2 2.18 1.94 11.75 12.6 11.71 12.6 5.216 5.58 12.64 ~13
10 1.95 1.75 6.58 6.9 4.95 5.2 3.88 4.08 10.94 ~11

shows that there is a small deviation between the cal- reactor core as average and hot channels. The
culated and referenced safety margins. Except for the RELAP5 steady-state results for the clad and coolant
ONB margin, the slight difference between the calcu- axial temperature distribution in comparison with
lated and referenced safety margins returns to the devi- those of the THDSN program are shown in figs. 2 and
ation in the average heat flux demonstrated in the 3, respectively. In this comparison, the THDSN pro-
above paragraph. The ONB margin of the present gram uses the Dittus Boelter correlation, eq. (14), for
study is nearly 12% higher than the reference value. calculating the heat transfer coefficient. The figures
This disagreement also returns to the deviation in the show that the THDSN temperature results are in good
ONB heat flux, as demonstrated previously. The ap- agreement with the RELAP5 results.
proximate sign that appears beside the referenced val-
ues of critical velocity in tab. 3 means that values are
determined by hand from figures in ref. [11]. Comparison with the PARET code

During the conversion of the Pakistan Research

Temperature distribution Reactor-1 (PARR-1) from HEU to LEU with an up-
grade from 5 MW to 10 MW, all the thermal-hydraulic
TECDOC-233 doesn’t contain the axial distri- and safety parameters were calculated at the elevated
bution of clad/coolant temperature in the hot or aver- power of 10 MW and by using the PARET code. For
age channel for the 2 MW or 10 MW RR. To verify the purposes of verification of the THMSN program, its
THDSN program results of temperature distribution, results regarding the safety limits and margins of
the RELAP5 model used in ref. [10] for 10 MW PARR-1 were compared with PARET results given in
benchmark reactor analysis is modified to simulate the refs. [16, 17]. The PARR-1 data required for the pres-
26 Nuclear Technology & Radiation Protection –1/2008

Figure 2. Clad temperature in the hot channel Figure 3. Coolant temperature in the hot channel

ent calculations are quoted from those references and tions. On the other hand, when the comparison
tabulated in tab. 4. between THDSN, RELAP5 and TECDOC-233 is
In order to take into consideration the engineer- considered, it can be concluded that PARET has over-
ing factor shown in tab. 4, in our calculations, the ra- estimated the heat transfer coefficient and that ad-
dial power factor used in THDSN is taken to be equal versely affects the peak temperatures and partially
to Fr × Feng. The core flow is considered to be down- overestimates the ONB heat flux.
ward. The Dittus-Boelter correlation (14) is chosen
to calculate the heat transfer coefficient. Therefore,
the results of THDSN compared to those of PARET Comparison with the TERMIC program
are shown in tab. 5. THDSN results, with and without
the engineering factor (Feng), are tabulated in tab. 5. The MTR reactor is an open-pool, plate-type
It appears that, except for the steady-state peak clad fuel element research reactor. The core is cooled and
and centerline temperatures and the ONB heat flux, moderated by light water and reflected by beryllium
there is a good agreement between the THDSN and and light water. Typical MPR 22 MW research reac-
PARET results. The agreement between the peak heat tor data are tabulated in tab. 6 [18]. Steady-state TH
flux and the coolant temperature rise in the hot chan- calculations of the safety limits and margins regard-
nel means that the treatment of the engineering factor ing this MTR reactor were performed by using the
in the two programs is similar. The THDSN predic- TERMIC program of the MTR_PC package, as tab-
tion of peak clad and centerline temperatures are, re- ulated in that reference. For the purpose of verifica-
spectively, 6.5% and 7.8% higher than the PARET tion, steady-state TH calculations for that reactor
ones. The ONB heat flux predicted by THDSN is are performed with THDSN and TERMIC programs
more conservative, being 12.7% less than that of using the data in tab. 6 and the results tabulated in
PARET, because it is calculated at worst core condi- tab. 7.
Table 4. Pakistan research reactor PARR-1 data [11, 12]
Property Value Propery Value
Core power [MW] 10 Core inlet temperature [°C] 38.0
No. of standard fuel elements (SFE) 29 Pressure at core end (bar absolute) 1.61
No. of control fuel elements (CFE) 5 Core pressure drop (bar) 0.2
No. of plates in SFE 23 Clad thermal conductivity [W/m°C] 180.0
No. of plates in CFE 13 Meat thermal conductivity [W/m°C] 53.6
Channel flow width [cm] 6.692 Radial power factor (Fr) 2.228
Channel heated width [cm]* 6.30 Axial power factor (Fa) 1.303
Fuel plate heated length [cm] 60.0 Engineering factor (Feng) 1.584
Fuel plate total length [cm] 62.6 Total flow rate [m3/h] 950
Channel thickness [cm]* 0.21 Channel coolant velocity [m/s] 2.46
Clad thickness [cm]* 0.0381 Fuel element nozzle diameter [cm] 6.0
Meat thickness [cm]* 0.051 Fuel element nozzle length [cm] 18.0
* These date not present in those references and quoted from ref. [5] for 10 MW benchmark reactor
A. Khedr: Thermal-Hydraulic Fortran Program for Steady-State Calculations of Plate-Type Fuel ... 27

Table 5. Comparison between PARET and THDSN

PARET, [15] THDSN
Parameter Feng = 1.584 Feng = 1.584 Feng = 1.0
Steady-state parameters
Average heat flux (AHF) [W/cm2] 18.1 18.07 18.07
Peak heat flux [W/cm2] = AHF Fa Fr Feng 83.4 83.094 52.458
Coolant temperature rise [°C]
– Average channel 9.4 9.54 9.54
– Hot channel 33.6 33.58 21.23
Peak clad surface temperature (hot channel) 102.47 109.167 84.644
Peak centerline temperature (hot channel) 104.64 112.796 86.935
Onset of nucleate boiling (ONB)
Average heat flux [W/cm2]
– Bergles and Rohsenow 25.52 20.548 35.28
Onset of flow instability (OFI)
Peak heat flux [W/cm2]
– Forgan (h = 48) 138.0 138.297 138.193
Departure from nucleate boiling (DNB)
Peak heat flux [W/cm2]
– Labunstov 326 325.469 325.469
– Mirshak 257 255.305 255.305
Critical velocity [m/s] 10.5 10.59 10.56
Safety margins
Margin to ONB 1.4 1.137 1.953
Margin to OFI (Forgan correlation) 1.6 1.664 2.915
Margin to DNB
– Labuntsov 3.9 3.917 6.204
– Mirshak 3.1 3.072 4.867

Table 6. MTR research reactor data [13]

Property Value Propery Value
Core power [MW] 22 Core inlet temperature [°C] 40.0
No. of fuel elements (FE) 29 Core exit pressure (bar absolute) 2.09
No. of plates in FE 19 Clad thermal conductivity [W/m°C] 180.0
Channel flow width [cm] 7.0 Meat thermal conductivity [W/m°C] 13.6
Channel heated width [cm] 6.40 Radial power factor 2.22
Fuel plate heated length [cm] 80.0 Axial power factor 1.35
Fuel plate total length [cm] 84.0 Total flow rate [m3/h] 1950
Channel thickness [cm] 0.27 Channel coolant velocity [m/s] 4.7
Clad thickness [cm] 0.04 Fuel element nozzle diameter [cm] 6.0
Meat thickness [cm] 0.07 Fuel element nozzle length [cm] 18.0

The comparison in tab. 7 shows that the THDSN CONCLUSIONS

results are more conservative than the TERMIC ones.
The THDSN results of clad temperature are 6.44% A steady-state, best estimate FORTRAN pro-
higher, while the ONB margin is 31.48% lower than gram for TH calculations of plate-type fuel RR called
corresponding TERMIC values. Also, THDSN THDSN has been built. The program simulates the re-
Mirshak calculated DNB is 21.67% under the actor core as two channels, an average and a hot chan-
TERMIC value. On the other hand, TERMIC overesti- nel. The program calculates the axial distribution of
mates the ONB margin because it uses a less conserva- coolant, clad and fuel meat temperatures, heat flux and
tive correlation, Foster and Greif, and overestimates core pressure drop. Also, it uses well-known correla-
the heat transfer coefficient because it calculates the tions for calculating safety parameters and margins
coolant properties in the Dittus-Boelter correlation at against the critical phenomena such as the onset of nu-
film temperature and not at bulk temperature [19]. cleate boiling, the departure of nucleate boiling and
Also, TERMIC considers the exit sub-cooling (DTsub) flow instability. The program gives the user possibility
in terms of the Mirshak correlation, eq. (20), inde- to choose one of six correlations for the heat transfer
pendent variable and that greatly overestimates the coefficient. All water properties are simulated in the
DNB margin. program by temperature and/or pressure dependant
28 Nuclear Technology & Radiation Protection –1/2008

Table 7. Comparison between TERMIC and THDSN Dhe – heated equivalent diameter of the channel, [cm]
Parameter THERIC THDSN Dy – channel hydraulic diameter, [cm]
FA – axial power factor, [–]
Steady-state parameters FD – factor less than unity, [–]
Coolant temperature rise [°C]
– Average channel 10.0 10.85 FNUC – nuclear power factor (=F1 F2), [–]
– Hot channel 22.78 24.11 FR – radial power factor, [–]
Peak clad surface temperature FPTN – fuel plate total number, [–]
jj (hot channel) 94.2 100.27 f – friction factor [–]
Peak centerline temperature Gch – channel mass flux, [kgm–2s–1]
jjj(hot channel) 117.79 h – heat transfer coefficient, [Wm–2°C–1]
Average heat flux (AHF)
jjj[W/cm2] 39.0 39.0 kc – clad thermal conductivity, [Wm–1°C–1]
Peak heat flux km – meat thermal conductivity, [Wm–1°C–1]
jjj[W/cm2] = AHF Fa Fr* 117.0 117.0 le – extrapolated length, [cm]
OFI peak heat flux [W/cm2] lh – fuel heated length, [cm]
Forgan 321.7 328.35 Lnoz – nozzle length, [cm]
Lp – total fuel plate length, [cm]
DNB critical heat flux [W/cm2] mch – channel mass flow rate, [kgs–1]
– Labunstov Nu – Nusselt number, [–]
468.047
– Mirshak P – absolute channel pressure, [bar]
445.3 348.819
– Bernarth Pcr – coolant critical pressure, [bar]
424.9
Safety margins DPch – channel pressure losses [bar]
Margin to ONB Pr – Prandtl number, [–]
– Bergles and Rohsenow 1.48 q – heat flux, [Wcm–2]
– Foster and Greif 2.16 qC – critical heat flux at DNB, [Wcm–2]
Margin to OFI qFI – peak heat flux at flow instability, [Wcm–2]
– Forgan 2.75 2.807
Margin to DNB q' – local power density, [Wcm–3]
– Labuntsov 4.001 qc¢ – maximum power density, [Wcm–3]
– Mirshak 4.235 2.982 Re – Reynolds number, [–]
– Bernarth 4.229 r – water density, [kgm–3]
SVR – fuel meat surface to volume ratio, [cm–1]
Tc – clad temperature, [°C]
polynomials produced from curve fitting to the pub- Tf – coolant temperature, [°C]
lished data. The program is constructed from a number Tm – meat centerline temperature, [°C]
of subroutines, each one calculating a certain distribu- Tsat – coolant saturation temperature, [°C]
tc – clad thickness, [cm]
tion or safety limit, all of them called from the program tch – channel thicknes, [cm]
main. tm – fuel meat thickness, [cm]
The THDSN program is verified by comparing tp – fuel plate thickness, [cm]
its results against commercial and customized pro- tw – water channel thickness, [cm]
grams. Also, the program is verified by comparing its Vch – channel coolant velocity, [ms–1]
Vnoz – coolant velocity in the fuel element nozzle, [ms–1]
results for 2 MW and 10 MW MTR IAEA reactors Vo – water velocity, [ms–1]
with those published in IAEA TECDOC-233, Appen- WA – active core flow rate (=FD WT), [m3h–1]
dix-A. A good agreement has been found with IAEA WT – total core flow rate, [m3h–1]
TECDOC, with the exception of the ONB limit which w – total channel width, [cm]
was overestimated by a narrow margin. Also, good wh – channel heated width, [cm]
wm – meat width, [cm]
agreement with PARET, except for the clad tempera- z – axial location, [cm]
ture and the ONB limit, was established. The THDSN
clad temperature turned out to be 6.5% over and the
ONB 12.7% under that calculated by PARET. This un- ABBREVIATION
derestimation in the ONB may be returned to the
THDSN if the estimated value is at the channel worst DNB – departure of nucleate boiling
condition. Another explanation may lie in the fact that DNBR – departure of nucleate boiling ratio
the PARET overestimates the heat transfer coefficient FE – fuel element
which affects the clad temperature and the ONB limit. FI – flow instability
HEU – highly enriched uranium
The comparison to TERMIC also shows that THDSN
LEU – low enriched uranium
uses more conservative correlations than those used MTR – material test reactor
by TERMIC. Also, TERMIC overestimates ONB and OFI – onset of flow instability
DNB values. ONB – onset of nucleate boiling
ONBR – onset of nucleate boiling ratio
RR – research reactor
TH – thermal hydraulic
NOMENCLATURE

Ach – channel cross-sectional area, [cm2]

Cp – coolant specific heat, [kJkg–1 °C–1]
A. Khedr: Thermal-Hydraulic Fortran Program for Steady-State Calculations of Plate-Type Fuel ... 29

ACKNOWLEDGEMENT [17] Bokhari, I. H., Israr, M., Pervez, S., Thermal Hydrau-
lic and Safety Analysis for Pakistan Research Reac-
tor-1, 22nd International Meeting on Reduced Enrich-
I wish to express my deep appreciation and sin- ment for Research and Test Reactors (RERTR),
cere gratitude to Prof. Galal El-Din Orabi of the De- Budapest, Hungary, October 4-8, 1999, http:// www.
partment of Fuel Cycle and Material Safety of the retr.anl.gov/Web1999/Abtracts/38bokharia99.html
[18] ***, INVAP MPR document Number
Egyptian Atomic Energy Authority for his valuable 0767-0750-3TPTH-359-1O, Core Thermo Hydrau-
comments, corrections, as well as his proposal regard- lics- Commissioning stage, 17 March 1998
ing the structure of this paper. [19] El-Morshdy, S., Tallha, K., Khedr, A., Determination
of Thermal-Hydraulic Safety Margins of ETRR-2 for
Normal Operation, Internal Report, Atomic Energy
Authority, Cairo, Egypt, 2005
REFERENCES

[1] ***, IAEA Safety Standards, Safety of Research Re-

actors, Safety Requirements No. NS-R-4, Vienna,
2005
[2] Obenchain, C. F., PARET-A Program for the Analysis
of Reactor Transient, ACE Research, IDO-17282,
Idaho National Laboratory, Idaho Fals, Id., USA, Jan-
uary 1969
[3] ***, RELAP5/MOD3 Code Manual, Volume 2: User’s
Guide and Input Requirements, NUREG/CR-5335,
INEL-95/0174, Idaho National Engineering Labora-
tory, Idaho Fals, Id., USA, June 1995
[4] Khedr, A., D’Auria, F., Nodalization Effects on
RELAP5 Results Related to MTR Research Reactor
Transient Scenarios, Nuclear Technology & Radia-
tion Protection, 20 (2005), 2, pp. 3-9
[5] Khedr, A., Adorni, M., D’Auria, F., The Effect of
Code User and Boundary Conditions on RELAP Cal-
culations of MTR Research Reactor Transient
Scnarios, Nuclear Technology & Radiation Protec-
tion, 20 (2005), 1, pp. 16-22
[6] Kater, H., El-Maty, T., El-Morshdy, S., Thermal-Hy-
draulic Modelling of Reactivity Accidents in MTR
Reactors, Nuclear Technology & Radiation Protec-
tion, 21 (2006), 2, pp. 21-32
[7] Housiadas, C., Lumped Parameters Analysis of Cou-
pled Kinetics and Thermal-Hydraulics for Small Re-
actors, Annals of Nuclear Energy, 29 (2002), 11, pp.
1315-1325
[8] El-Messiry, A., Accidents of Loss of Flow for the
ETRR-2 Reactor: Deterministic Analysis, Nukleonika,
45 (2000), 4, pp. 229-233.
[9] Louis, J., Simple Reactor Model Simulation of a
LOFT ATWS Events, Nuclear Technology, 6 (1983),
1, pp. 25-32
[10] ***, IAEA Safety Practice, Operational Limits and
Conditions for Research Reactors, Safety Series No.
35-P1, Vienna, March 1994
[11] ***, TECDOC-233, Research Reactor Core Conver-
sion from the Use of Highly Enriched Uranium to the
Use of Low Enriched Uranium Fuels Guidebook,
IAEA TECDOC-233, 1980
[12] El-Wakil, M. M., Nuclear Heat Transport, Interna-
tional Textbook Company, Cranton, Pa., USA, 1971
[13] Prasuhn, A. L., Fundamentals of Fluid Mechanics,
Prentice-Hall, Inc., Englewood Cliffs, N. J., USA
[14] Glasstone, S., Sasonske, A., Nuclear Reactor Engi-
neering, Van Nostrand Reinhold Company, New
York, 1967
[15] Holman, J. P., Heat Transfer, Mc Graw-Hill Book
Company, New York, USA, 1989
[16] Bokhari, I. H., Steady-State Thermal Hydraulic and
Safety Analysis of a Proposed Mixed Fuel (HEU&LEU)
Core for Pakistan Research Reactor, Annals of Nuclear
Energy, 31 (2004), 11, pp. 1265-1273
30 Nuclear Technology & Radiation Protection –1/2008

Ahmed KEDR

FORTRANSKI PROGRAM ZA TERMOHIDRAULI^KE PRORA^UNE STABILNOG

STAWA ISTRA@IVA^KIH REAKTORA SA PLO^ASTIM GORIVOM

Sigurnosna procena istra`iva~kih i energetskih reaktora obavqa se neprekidno tokom

wihove upotrebe i zahteva proverene i vaqane programe. U ~itavom svetu programi za energetske
reaktore dobro su utemeqeni i utvr|eni merewima i podacima sa odgovaraju}ih eksperimentalnih
postrojewa. Ovi kodovi obi~no su sofisticirani, a od korisnika zahtevaju posebne ve{tine i
tro{e mnogo ra~unarskog vremena. S druge strane, ve}ini programa namewenih istra`iva~kim
reaktorima jo{ uvek nedostaju podaci radi proveravawa i potvr|ivawa. Otuda je od koristi za
svako regulatorno telo da razvije svoje sopstvene kodove za pregled i ocenu istra`iva~kih
reaktora.
U ovom radu predstavqen je THDSN – jednostavan jednodimenzionalni fortranski pro-
gram za termohidrauli~ke prora~une stabilnog stawa istra`iva~kih reaktora sa gorivom
plo~astog tipa. Pored prora~una raspodela temperatura i gradijenata pritiska goriva i hladioca
u sredwem i vru}em kanalu, programom se prora~unavaju sigurnosne granice i margine na osnovu
pojava kriti~nosti koje se susre}u kod istra`iva~kih reaktora, kao {to su: fenomen kqu~awa,
kriti~ni toplotni fluks i nestabilan tok. Kori{}ewe su poznate termohidrauli~ne korelacije
za prora~un sigurnosnih parametara i vi{e formula za koeficijente prenosa toplote. Program
THDSN potvr|en je pore|ewem wegovih rezultata za ben~mark reaktore od 2 MW i 10 MW, sa
rezultatima objavqenim u MAAE izdawima, a sa kojima je na|ena dobra saglasnost. Tako|e,
upore|eni su rezultati ovog programa sa objavqenim rezultatima drugih programa, kao {to su
PARET i TERMIC.

Kqu~ne re~i: istra`iva~ki reaktor, termohidrauli~ki prora~un, stabilno stawe,

jjjjjjjjjjjjjjjjjjjjjjjjsigurnosni parametri, Fortran program