Nuclear Engineeririg and Design 52 (1979) 343-347 © North-Holland Publishing Company
T H E R M A L - H Y DRAULIC ANALYSIS TECHNIQUES FOR
AXISYMMETRIC PEBBLE BED REACTOR CORES K.R. STROH Reactor and Advanced Heat Transfer Technology Group, Energy Division, Los Alamos Scientific Laboratory, Los Alamos, NM 87545, USA R.J. JIACOLETTI Design Engineering Division, Los Alamos Scientific Laboratory, Los Alamos, NM 87545, USA and H.G. OLSON Mechanical Engineering Department, Colorado State University, Fort Collins, CO 80523, USA Received 30 June 1978
A mathematical model for the analysis of coupled thermal-hydraulic problems in steady-state pebble bed nuclear reactor cores is presented. The bed has been treated macroscopically as a generating, conducting porous medium. The model uses a nonlinear Forchheimer-type relation between the coolant pressure gradient and mass flux, and new coefficients of the viscous and inertial loss terms are presented. The remaining equations in the model make use of continuity and thermal energy balances on the solid and fluid phases. None of the usual simplifying assumptions such as constant properties, constant velocity flow or negligible conduction and/or radiation are used. A computer program based on this model has been constructed; it has been validated by comparing predictions with measured values of previous experiments. Validation of the nonlinear fluid flow model is reported in a companion paper.
fission energy from the pebbles and rejects heat to a secondary system before being returned by the loop for another pass through the core. For a given maximum fuel temperature, the gas outlet temperature for a PBR can be considerably greater than that obtained using the prismatic block high-temperature gas-cooled reactor (HTGR) core [1 ]. Incentives to develop nuclear reactors of this type are high because they represent an attractive, and perhaps singular, alternate high-temperature-process heat and power source for fossil fuel conversion, primary metal production and other large high-temperature-process applications. Assembly o f a mathematical model, and the derivation of equations that describe the coupled t h e r m a l hydraulic behavior o f a PBR core, are presented in this paper. In a companion paper [2], predictions by
1. Introduction it is sta_ndard practice in the petroleum and chemical industries to use packed beds of catalyst-containing particulates to improve yields in flowing reactants. Because of the random nature o f the bed's constituents, it is not possible to formulate a detailed physical description when predicting the reactants' flow behavior. The typical calculational approach is to treat the entire bed macroscopically as a porous medium, usually assuming constant velocity fluid flow, and then to use empirically derived relationships between variables. In the nuclear pebble bed reactor (PBR) concept, the cylindrical core similarly contains a random bed o f small, spherical, fuel-moderator elements. The downward-flowing helium primary coolant, under high pressure to increase its heat transfer capabilities, extracts 343
344
K.R. Stroh et al. / Thermal-hydraulic analysis for reactor cores
the computer code based on this model are compared with data obtained from a full-scale mockup of the Oak Ridge Pebble Bed Reactor Experiment (PBRE).
where the products f l V and f2 V2 represent the vis6ous and inertial loss terms, respectively. Koida [21 ] has developed expressions for fl and f2 so that
2. Previous work
fl =Ao/eR~,
(2)
Considerable literature exists regarding modeling of heat transfer phenomena in chemical catalytic packed bed reactors [3-8], which can assist when assembling a model for the thermal analysis of a nuclear PBR. The previous work has limited applicability however, because it generally assumes that the fluid velocity is constant, one-dimensional and known. Recently there has been a growing interest in developing modeling equations to represent flow maldistribution in packed beds. The mathematical models l~or these systems have used the differential, vectorial form of the Ergun equation [9-15]. Szekely and Povermo [16] have presented direct experimental evidence to establish the validity of this approach. To date, none of these flow analyses have been coupled to an accurate thermal model, and most are constrained by a narrower range of Reynolds number, and other parameter variations, than is required for analysis of a nuclear PBR. A limitation of the Ergun equation [17,18] is that it overpredicts the pressure drop across a random bed of uniform spheres at the high Reynolds numbers of interest in PBR applications.
f2 = Bp/e:Rh,
(3)
where A and B are the empirical coefficients that are constant for a given particle shape, R h is the effective hydraulic radius for the flow within the packed bed (m), e is the bed void fraction,/a is the fluid dynamic viscosity (Pa s), and p is the fluid density (kg/m3). Rh for a randomly packed bed of spheres is given by [22] Rn = ~dp e/(1 - e),
where dp = pebble diameter (m). Numerical values for the coefficients A and B of 24.5 and 0.1754, respectively, were calculated using Barthels' [23] friction factor correlation for high Reynolds number flows in packed sphere beds. If the fluid pressure drop and temperature rise across the bed result in an appreciable change in fluid density, it is more convenient to calculate the mass flux, G = p V, rather than the fluid velocity. Using G, eq. (1) becomes VP =
-G(g
1
+g2 IGI),
(5)
where
3. Hydraulic model In formulating this model, the core volume was treated as a heat generating, conducting porous medium characterized by its uniform pebble diameter, the core diameter to pebble diameter ratio and its void fraction distribution. In a PBR, the bed volume and diameter ratio are both sufficiently large to allow characterization by a statistically averaged void fraction distribution, which allows calculation of local macroscopic values of the thermal-hydraulic variables. When the Reynolds number of the flow in porous media is high, the nonlinear relationship between fluid velocity and pressure drop can be described by a Forchheimer-type equation [19,20], relating the pressure, P, and the velocity, V, so that VP = - V ( f x +I'2 IVI),
(4)
(1)
gl = f l / P ,
(6)
gz = f 2 / o 2 .
(7)
Since curl • grad = O, the pressure variable can be eliminated from eq. (5); thus V X [ - G ( g l +g=lGI)] = 0 .
(8)
It should be noted, however, that the pressure distribution still affects the mass flux through the fluid density part o f g I andg 2. The model's equations have been formulated in axisymmetric cylindrical coordinates because they allow a realistic representation of actual physical conditions and facilitate coupling of this model to the current neutronics model. Manipulation of eq. (8)
K.R. $troh et al. / Thermal-hydraulic analysisfor reactor cores yields the following scalar equation:
345
of the nondimensional variables in the following form: a
aq/
~z,[-(gl +g2lGl)Gr] - ~r, [-(gl +g2lGDGz] =O ,
(9) where z' and r' are the axial and radial coordinates. The equation is made nondimensional by the application of the following definitions: G* _ IGI
el
=
g2GxN ' z' R z =-L, a = L ,
g~ _ g2
GtN '
g2rN '
r : - -r' R' (10)
where R is the bed radius (m), L is the bed height (m), and subscript IN denotes the inlet plenum values. Introducing these definitions into eq. (9) gives a nondimensional scalar equation in terms of G* and Gz, the unknown radial and axial components of G *, D , , , a"-z [(~ + a ) g : a r ]
la
+a
oz
=0.
aa¢/ r az'
z =0.
(12)
lamb G; - r Or .
az Lr
+G )g2
~z
-
(151
where ~ is the dependent variable, ~b is the stream function, and a m, b m, c m and d m are functions used as required to make eq. (15) and the equation transformed identical. For eq. (14), the following functions are used: O =if,
am = 0 ,
b m = ( ~ + G )g2,
cm = 1 .
ar Lra:
+G )g2
(16)
Eqs. (14) and (15) were expanded and set equal to each other to obtain the expression for dm. The resulting general nondimensional form of eq. (14) becomes
+
.
=0
(/j+G*)g~-~-r = 0 .
(17)
Once the stream function field is obtained, the pressure distribution follows. Computing the divergence of eq. (5), with appfication of the continuity equation, and applying the nondimensional definitions yields the following scalar equation:
ar /
azz~r-~--z]- L
P,~f
/a a
(13)
When so defined, the stream function inherently satisfies the continuity equation. Combining eqs. (1 l) and (13) gives a nonlinear, elliptic partial differential equation
_Zrl
+rd¢=O,
(11)
For the purposes of computation, the coupled system of eqs. (11) and (12) is replaced by a single equation of higher order for the nondimensional stream function, ~k, defined by
a;-
m~ r r (cm
•
-aOr [(~+ a*)g~az]
The continuity equation, V. G = 0, in nondimensional form and axisymmetric cylindrical coordinates is
r or
-
a
X [~-r ~ [(~ + G*)g;] - az-- at-- [(~ + G ' ) g ;
1}
= O,
(18) where P* = P/Pra, the nondimensional pressure. Eq. (18) is in the form ofeq. (15) where ¢=P*,
a re:O,
bm = c ¢ = l ,
"
(14)
Because the coefficients are a function of temperature, and the temperatures in turn depend on the mass flux, a complete thermal-hydraulic description results in a set of coupled elliptic partial differential equations like eq. (14). The numerical solution method selected [24] requires that the equations be expressed in terms
dm
L rPref
J[ar az [(~ + G*)g~]
a~ka, } az ~ [(~ +G*)g~] .
(19)
Using the coupled set of eqs. (17) and (18), one can solve for the pressure and velocity distributions within the conf'mes of the bed.
346
K.R. Stroh et aL / Thermal-hydraulic analysis for reactor cores
4. Thermal model The rate of heat transfer in a generating, conducting porous medium with a flowing fluid phase is controlled by a number of mechanisms, including bulk movement of the fluid, conduction in both the solid and fluid phases, convective transfer between phases, dispersion of the fluid in the interstices of the porous medium and radiant exchange. When different mechanisms have a common driving potential, they may be combined by applying an effective heat transfer coefficient. The macroscopic temperature gradient in the solid phase drives a complex group of heat transfer mechanisms, which operate both in parallel and in series. These include thermal conduction through the solid pebbles, conduction through the stagnant fluid film near the contact point of two adjacent pebbles, thermal conduction by contact between pebbles and, for a gaseous fluid, radiant heat exchange. The net effect of these mechanisms can be expressed by an effective thermal conductivity coefficient for the solid phase, kse. The Static term, k ° , derived by Kunii and Smith [25] is used in the model. In the fluid phase, heat is transferred between fluid regions at different temperatures by two mechanisms: molecular conduction and turbulent dispersion in the interconnected interstitial voids of the solid phase. These combined effects are approximated by a single effective thermal conductivity coefficient for the fluid phase [3] given by
kfe = ekf + p c E ,
mal energy balances on the two phases. These balances are based on the fluid superficial velocity (a velocity based on a void fraction of 1.0) and the total crosssectional area. Convective coefficients, generation rates and thermal conductivities are corrected for the proper fraction of the bed. A thermal energy balance on the solid phase results in the following equation: V • (-kseVts) + h a v ( t s
(21)
tf) - q = O,
where h is the convective heat transfer coefficient (W/m e K), av is the pebble surface area per unit volume of the bed (m-l), tf is the local value of the bulk temperature of the interstitial fluid (K), t s is the local value of the pebble average surface temperature (K) and q is the energy release rate per unit volume of the bed (W/m3). The terms in eq. (21) represent, from left to right, effective conduction in the solid phase, convective transfer between the phases and internal generation from fission and decay. Introducing the following nondimensional variables: Ts-
Ks
=
ts
tflN
tf
,
Tf . . . . . . , if iN
kse --,
c
. av =avL,
kflN
. _ CGINL
,
kflN h*
-
hL
kflN
,
Kf-
q
.
-
kfe
kflN
,
qL 2
kflN/fiN (22)
allows eq. (21) to be transformed by the same methods and definitions used to obtain eq. (17). The result, in axisymmetric cylindrical coordinates, is
(20)
where kf is the fluid molecular thermal conductivity (W/m K), c is the fluid specific heat capacity at constant pressure (J/kg K), and E is the turbulent thermal diffusivity of the flowing fluid phase (m2/s). Because the turbulent thermal diffusivity is anisotropic, the fluid effective thermal conductivity is also anisotropic. The diffusivity values used in the computer program are those reported in Oak Ridge PBRE data [26]. It is assumed that the required turbulent thermal diffusivity is numerically equal to the turbulent nrass diffusivity. The distributions of local bulk fluid temperature and average pebble surface temperature can be obtained by the solution of equations derived from ther-
-
ar -l + r[h*av(Ts
K -
T~) -
q*] = 0 ,
(23)
which is in the appropriate general nondimensional elliptic form required by the computer program. A similar treatment of the thermal energy balance for the fluid phase yields
-- L[rKfz ~Tfl -
0 r ~[DKfr ~Tf~ ~ *
*Tf ~-z -
rh*av(Ts
- Tf) = 0 ,
~
(24)
where it should be noted that the effective thermal
,
K.R. Stroh et al. / Thermal-hydraulic analysis for reactor cores conductivity coefficients are anisotropic. Values for the convective heat transfer coefficient, h, are obtained from the Jeschar correlation [27] for the Nusselt number, which is applied locally:
NNu-hdp
kf - 2 " 0 +
I ~
Nn
e] °'s
+0.005
NRe e ,(25)
where the Reynolds number, based on the superficial velocity and the particle diameter, dp, is given by NRe = Gdp/la.
(26)
The Nusselt number correlation is reported to be valid for Reynolds numbers between 250 and 5.5 × 104. The effective thermal conductivities of the solid and fluid phases are functionally related to the local film temperature, usually the arithmetic average of t s and tf. This dependence of the coefficients on the dependent variables makes eqs. (23) and (24) nonlinear. The computer code based on this model solves eqs. (17), (18), (23) and (24) as a coupled system o f nonlinear elliptic partial differential equations. Boundary conditions for the solution of a selected physical system are specified at r = 0, r = 1 and both z limits (inlet and outlet), which need not be parallel. The numerical solution uses finite difference forms of the coupled system of equations, with the method of successive relaxation used to solve the resulting nonlinear finite difference equations. Details o f this method and a listing of the computer program, which was developed by the first-named author, can be found in ref. [28].
5. Conclusions A mathematical model has been developed for the analysis of coupled t h e r m a l - h y d r a u l i c problems in steady-state pebble bed nuclear reactor cores. A computer program based on this model has been constructed and has been validated by comparing predictions with measured values o f previous experiments. Validation of the nonlinear fluid flow model is reported in a companion paper [2], using data from the Oak Ridge PBRE full-scale mockup [26]. The mathematical model (and resulting computer program) does not use any of the usual simplifying assumptions such as constant velocity flow, constant properties, or
347
negligible conduction and/or radiation. Coupled therm a l - h y d r a u l i c test problems are now being run using computed two-dimensional fission power profiles, and will be reported as they are completed and analyzed.
References [ 1] K. Kugeler, M. Kugeler, H.F. Niessen and H. Hohn, Nucl. Eng. Des. 34 (1975) 15. [2] K.R. Stroll, H.G. Olson and R.J. Jiacoletti, Nucl. Eng. Des. 52 (1979) 349. [3] E. Singer and R.H. Wilhelm, Chem. Eng. Progr. 46 (1950) 343. [4] L. Green, Jr., J. Appl. Mech. 19 (1952) 173. [5] L.S. Dzung, Int. J. Heat Mass Transfer 1 (1960) 236. [6] W.B. Argo and J.M. Smith, Chem. Eng. Progr. 49 (1953) 443. [7] N.G. Karanth and P. Hughes, Chem. Eng. Sci. 29 (1974) 197. [8] W.U. Choudhury, dissertation, University of Wisconsin (1969). [9] S. Ergun, Chem. Eng. Progr. 48 (1957) 89. [10] V. Stanek and J. Szekely, Can. J. Chem. Eng. 50 (1972) 9. [l 1] V. Stanek and J. Szekely, Can. J. Chem. Eng. 51 (1973) 22. [12] V. Stanek and J. Szekely, AIChE J. 20 (1974) 974. [13] V.S. Shvydkii, Y.M. Gordon, Y.G. Yaroshenko and V.B. Shcherbatskii, Steel USSR 4 (1974) 662. [14] M. Choudhary, M. Propster and J. Szekely, AIChE J. 22 (1976) 600. [15] M. Choudhary, J. Szekely and S.W. Weller, AIChE J. 22 (1976) 1021. [16] J. Szekely and J.J. Povermo, AIChE J. 21 (1975) 769. [17] W.H. Denton, C.H. Robinson and R.S. Tibbs, AERE (Harwell) Research Report AERE-R4346 (1963). [18] R.D. Bundy and H.W. Hoffmann, Oak Ridge National Laboratory Report ORNL-3523 (1963) p. 414. [191 N. Ahmed and D.K. Sunada, Proc. ASCE-HY6 (1969) 1847. [20] A.E. Scheidegger, The Physics of Flow Through Porous Media (Macmillan, New York, 1960). [21] N.U. Koida, Russ. J. Phys. Chem. 34 (1960) 375. [22] S. Debbas and H. Rumpf, Chem. Eng. Sci. 21 (1966) 583. [23] H. Barthels, Brennstoff-Wffrme-Kraft 24 (1972) 233. [24] A.D. Gosman, W.M. Pun, A.K. Runchal, D.B. Spalding and M. Wolfshtein, Heat and Mass Transfer in Recirculating Flows (Academic Press, London, 1969). [25] D. Kunii and J.M. Smith, AIChE J. 6 (1960) 71. [26] R.D. Bundy, Oak Ridge National Laboratory Report ORNL-TM-1975 (1966). [27] R. Jeschar, Arch. Eisenhfittenwesen 35 (1964)517. [28] K.R. Stroh, dissertation, Colorado State University (1978).