- Email: [email protected]
Three-dimensional simulation of the coupled convective, conductive, and radiative heat transfer during decay heat removal in an HTR
Nuclear Engineering and Design 237 (2007) 1923–1937
Three-dimensional simulation of the coupled convective, conductive, and radiative heat transfer during decay heat removal in an HTR Jan-Patrice Simoneau a,∗ , Julien Champigny a , Brian Mays a , Lewis Lommers b a
AREVA, 10 rue Juliette R´ecamier, 69456 Lyon Cedex 06, France b AREVA, 3315 Old Forest Road, Lynchburg, VA 24506, USA
Received 16 June 2006; received in revised form 14 March 2007; accepted 14 March 2007
Abstract A three-dimensional model has been constructed to simulate the passive heat removal in a modular prismatic-block high temperature reactor during a loss of active cooling accident. This model, developed using the STAR-CD general computational fluid dynamics code, solves the combined conductive, convective, and radiative heat transfer within a 30◦ section of the core and reactor vessel. To accommodate the different spatial scales, it uses homogeneous equivalent media to represent the coolant flow and the prismatic fuel blocks. A customized procedure that manages solving alternatively the dynamic and thermal fields permits the computation of very long transients, which typically are performed for 100 or more hours of simulated time. The global methodology and specific modeling procedures are explained, and key points of the CFD analysis are highlighted. Next, the results of several calculations are presented, and the physical phenomena represented are described. Two commonly investigated loss of active cooling scenarios are considered: depressurized conduction cooldown and pressurized conduction cooldown. The results for these two scenarios are compared to assess the effect of heat transfer via internal natural convection – which is negligible during the depressurized event – on the thermal behavior of the system. In addition, the evolution of the natural convection flow through the core and in the annular spaces is examined and discussed. © 2007 Elsevier B.V. All rights reserved.
1. Introduction Recent years have seen a renewed interest in the high temperature reactor (HTR) as a high-temperature nuclear heat source. This interest is driven largely by the possibility of using nuclear energy for the production of hydrogen. The most efficient means of hydrogen production, such as high-temperature electrolysis or thermochemical water splitting, require very high temperatures, which are only achieved by HTRs. Since greater efficiency is possible as the process temperature is increased, current HTR designs are being extended to produce even higher outlet temperatures, increasing from 850 ◦ C to temperatures approaching 1000 ◦ C. The 600 MW HTR considered here is a prismatic-block concept. The power conversion system is an indirect combined Brayton and Rankine cycle. This HTR design is capable of pas-
∗
Corresponding author. E-mail address: [email protected] (J.-P. Simoneau).
0029-5493/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2007.03.010
sive decay heat removal for enhanced inherent safety. Decay heat can be removed using only conduction and radiation from the fuel to the vessel and radiation and natural convection from the vessel to a water circuit located on the vault walls. The reactor is configured to ensure acceptable temperatures for the fuel particles under these conditions. In addition to conduction and thermal radiation phenomena, natural convection can also occur within the core during passive cooling. However, the convective heat transfer remains negligible if the circuit is depressurized (referred to as depressurized conduction cooldown, or DCC). On the contrary, if the circuit remains pressurized (referred to as pressurized conduction cooldown, or PCC), the internal natural convection provides noticeable heat transfer, and a different thermal response is observed in the core structures. Thus, PCC implies combined convective/conductive/radiative heat transfer. The combination of these heat transfer modes, each with very different dynamical characteristics, can pose a computational challenge. The objective of this paper is to present the implementation of a 3-D PCC computation and the resulting thermal response of
1924
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
the reactor. The computations presented are general exploratory calculations, not design calculations. 2. Key features of the reactor design 2.1. General characteristics The HTR concept considered here uses three circuits, operating on a combined Brayton-Rankine cycle (Petit, 2005). The primary circuit, consisting of pressurized helium, transfers the heat from the reactor core to an intermediate heat exchanger (IHX). The secondary circuit, working with a nitrogen-based fluid, drives a gas turbine to form a Brayton cycle that uses a turbo compressor coupled with an alternator. The final bottoming circuit is a water-steam system that works as a conventional Rankine cycle. The key parameters of the reactor are given in Table 1. The core consists of hexagonal graphite blocks, which are divided into an active core of 1020 hexagonal prismatic fuel elements, stacked in 102 columns of 10 blocks each, and a set of solid graphite reflector blocks. The active core is arranged in an annulus, with some of the reflector blocks (the inner reflectors) located inside the annulus and the remaining blocks (the outer reflectors) located between the active core and the reactor vessel. The core is held by a steel core barrel, and all structures are enclosed in a pressure vessel, which is connected to an adjacent vessel containing the IHX via a short cross-vessel. Figs. 1 and 2 present axial and radial cross-sections of the reactor vessel and core. 2.2. Description of the primary coolant flow Under normal operating conditions, cold helium travels through the outer annulus of the cross-vessel and enters the bottom part of the reactor vessel, where it cools the metallic support structures. It then flows up to the top of the core in a riser between the core and the pressure vessel, to limit heat losses and protect the pressure vessel from excessive temperatures. Starting from a plenum at the top of the vessel, the coolant flows downward through the core, where it is heated, and is collected in a second plenum below the core, before exiting through the inner annulus of the cross-vessel to the IHX vessel. Most of the helium flows through the coolant channels in the blocks that form the active core, but part of this flow bypasses
Fig. 1. Vertical cross-section of the reactor indicating the flow path of the helium coolant during normal operation.
Table 1 Key reactor parameters Parameter
Value
Core thermal power (MW) Helium mass flow rate (kg/s) Helium pressure (bar) Reactor inlet temperature (◦ C) Reactor outlet temperature (◦ C) Core inner diameter (m) Core outer diameter (m) Core height (m) Vessel outer diameter (m)
600 226 50 490 1000 2.96 4.84 8 7.7
Fig. 2. Horizontal cross-section of the reactor.
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
the channels by flowing through the gaps between the blocks. This bypass flow cools the reflectors, but does not contribute directly to the nuclear heat removal. It mixes with the main flow in the lower plenum before exiting the reactor vessel. Very complex flows (Simoneau and Champigny, 2004) occur in this outlet plenum due to three effects: injection of “hot” helium from the coolant channels, injection of colder helium from the gaps, and highly turbulent flow that forms in the wakes resulting from the presence of graphite columns, which support the core, in the transverse flow toward the plenum exit. The gaps between the blocks play a significant role in the flow distribution and hence the thermal behavior of the reactor. During normal operation, the bypass flow has undesired consequences, since it does not directly cool the fuel. Therefore, the temperatures at the coolant channel outlets are locally higher than the reactor outlet temperature due to this bypass. Decay heat removal without any primary flow, on the other hand, is a transient situation where the thermal inertia of the different regions of the reactor dominate. During the beginning of this transient, the reflectors act as a heat sink, which limits the maximum temperatures and protects the fuel blocks. The temperature distribution within the core at the beginning of the transient – and hence the amount of heat that the reflectors can absorb – is mainly a result of the amount of gap flow. For this study, the gaps between the graphite blocks are assumed to be uniformly 2 mm wide, the value for new blocks at ambient temperature. During reactor operation, however, the blocks experience size changes due to the effects of thermal expansion, which reduces the gap size, and irradiation of the graphite, which tends to shrink the graphite, increasing the gap size. Thus, the bypass flow behavior changes greatly from fresh to used blocks. The loss of primary flow with fresh fuel and a core of new graphite blocks leads to higher reflector temperatures and at the same time to lower residual power. In the calculations presented in this paper, a bypass flow of 10% was selected and prescribed for the model; however, other calculations (not reported here) with different bypass flow values have been performed. 2.3. The fuel blocks The fuel blocks are hexagonal right-regular prisms composed of graphite with holes that form the coolant channels and smaller holes into which the fuel compacts are inserted. The TRISOcoated spherical fuel particles are incorporated with a graphite matrix to form the cylindrical compacts. These components are depicted in Fig. 3. 2.4. The decay heat removal system In the event of a loss of primary flow, the decay heat can be removed by the reactor cavity cooling system (RCCS). This system consists of vertical water pipes, which are set on the concrete walls of the cavity surrounding the reactor vessel. These pipes, which are connected to a tank, are designed to maintain an average temperature of 65 ◦ C and to protect the concrete walls during normal operation.
1925
Fig. 3. The prismatic fuel block, compact, and TRISO particle.
The system described here spans several very different spatial scales, which a CFD calculation has to manage. In terms of scale, the solid parts range from the fuel particles (less than 1 mm in size) to the reactor vessel (several meters long), and the fluid volumes range from the channel diameter (16 mm) and the gap width (2 mm) to the cavity height (about 20 m). 3. Assumptions for decay heat removal transients The transient operating conditions considered here for a loss of primary flow are as follows. The initiating event is a loss of the primary forced flow. In addition, all heat removal mechanisms except for the reactor cavity cooling system are assumed to be unavailable. Two scenarios are then considered. In the first scenario, depressurized conduction cooldown (DCC), the primary circuit is depressurized, either via a large break or intentional depressurization. (The term “conduction” is adopted here because the behavior of the reactor is similar to the cooling of a solid body by conduction; nevertheless, radiation and, to a much more limited extent, convection both contribute to the heat removal process.) Since the density of the helium gas is low – on the order of atmospheric pressure, due to equilibrium with the air in the reactor cavity volume – the heat transfer by convection is relatively small compared to other phenomena. The other scenario, pressurized conduction cooldown (PCC), assumes that the coolant remains in the primary circuit. Thus, natural convection is much more important than in the previous case.
1926
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
The phenomenon of decay heat removal occurs on very large time scales (tens of hours), but the development of convection in the core channels occurs within a few minutes. Moreover, the characteristic times of the convection in some regions inside the reactor are smaller than one second (e.g., T ∼ L/U, with U ∼ 0.5 m/s and L ∼ 10 mm). Thus, to simulate such a system, a model must manage a wide variety of time scales. 4. The numerical model and code 4.1. Computational techniques
Fig. 4. The decay heat curve (main plot) and axial distribution of the power (inset).
The conditions in the DCC scenario result in higher fuel temperatures, since the transfer of thermal energy from the core to the heat sink is limited to the radiative and conductive processes in the blocks and the gaps between the components. The PCC scenario can be the more limiting case for the pressure vessel, since the heat transfer from the core to the vessel is increased by the natural convection. Furthermore, this situation can lead to hot spots in the upper part of the reactor. Both scenarios follow the same general timeline. Before the beginning of the transient, the reactor is operating at full power (600 MW), the RCCS is in operation, and 10% of the total coolant flow bypasses the active core. The power generated in the annular core is assumed to be distributed uniformly in the radial direction. The axial power profile is shown in Fig. 4. At the beginning of the transient (time t = 0 s), the primary helium circulation is stopped, and the protection system automatically drops the rods in response to the loss of main loop circulation. In the DCC scenario, the system also rapidly depressurizes, so that it can be considered completely depressurized when the transient begins. During the transient, the decay heat is dissipated in the core, while the RCCS continues to operate and serves as the ultimate heat sink for the system. Beginning at 30 MW (5% of full power), the decay heat decreases according to the curve in Fig. 4, which is a standard radioactive decrease for a core at end of cycle. The axial power profile is assumed to remain homothetic to its initial profile. After the loss of forced flow, natural circulation develops in the reactor vessel as a result of the buoyancy forces arising from the temperature distribution at normal operation, which serves as the initial condition of the transient. This natural convection occurs through the core channels and gaps, in the top collector, in the bottom collector, and in the annular spaces along the pressure vessel. A shutoff valve, which automatically closes when forced convection is lost, prevents circulation through the primary loop. Natural convection occurs in both the DCC and PCC scenarios, but as mentioned above, the effect of this convection is ignored in the DCC case, since is it very weak and the amount of heat transfered by it is negligible.
The general-purpose fluid-mechanics and heat-transfer code STAR-CD, version 3.2, (Computational Dynamics Limited, 2004) was used to simulate the decay heat removal transients. This code solves the equations of fluid flow and heat transfer on an unstructured mesh, which for this study consists mainly of hexahedral and prismatic cell volumes. The partial differential equations are discretized by the finite volume method, using first-order upwind differences for the spatial terms and a firstorder, fully implicit scheme for the temporal discretization. The matrices that appear in the solution of the discrete finite volume equations are inverted by the double conjugate gradient method with preconditioning provided by the algebraic multigrid method. This approach has proved to be very useful for solving the pressure correction equation in tall geometries, such as those that exist in the model described here. Two algorithms are used to solve the pressure linked equations. For the steady-state calculations that are used to provide the initial conditions for the transient, the classic SIMPLE algorithm (Patankar, 1980) is used. During the transient, the more sophisticated PISO algorithm (Issa et al., 1991) is used. The convergence of the latter algorithm is ensured by selecting sufficiently small time steps. Turbulence is modeled only during the steady-state calculation for normal operation, which has high flow velocities due to the forced circulation. In this case, the model uses STAR-CD’s standard linear k–ε model and wall functions for high Reynolds number flow, which are well suited for forced convection. Heat transfer by thermal radiation is modeled by face-to-face exchanges of energy using a discrete beam method to calculate the view factors. This approach assumes that all surfaces are gray and diffuse and that the helium gas does not participate in the radiative heat transfer. The equations that are solved for each type of calculation described here are summarized in Table 2. Table 2 Equations solved by the STAR-CD code Physical quantity
Initial steady-state
PCC transient
DCC transient
Flow field
Solved with SIMPLE Solved Solved with k–ε high Re Solved
Solved with PISO Solved Not solved (Laminar flow) Solved
Not solved
Temperature field Turbulence Density
Solved Not solved Not solved
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
4.2. Modeling philosophy Instead of explicitly modeling the entire reactor core, only the minimum section required to represent the core threedimensionally via symmetry – a 30◦ sector – was included in the model. Although this sector is sufficient to represent the core itself, a rigorous modeling of the outlet plenum would require at least a 180◦ sector. Therefore, the modeling approach used here leads to an approximation of the outlet flow; however, a detailed representation of this flow is not the objective of this study. The metallic structures of the reactor are explicitly modeled, but some details are neglected: the shroud tubes (for control rods) in the upper plenum, the vertical ribs of the lower supporting plates, and other very small parts. Similarly, the graphite posts in the outlet plenum are not modeled. The blocks are represented by a homogeneous equivalent solid medium, which includes both the fuel compacts and the
Fig. 5. The grids representing the fluid parts (top) and solid parts (bottom) of the explicit model.
1927
graphite web. This technique is described in further detail in Section 4.4. The gaps between the core blocks are explicitly modeled by mesh cells, with the following assumptions. The solid-fluid heat transfer and the pressure drop across the core are prescribed via coefficients (i.e., the boundary layer is not resolved). The radiation across gaps is explicitly modeled. Finally, the gap size is assumed to be uniform and equal to the “cold” value of 2 mm. As mentioned above, the space between blocks in an actual core would be expected to vary throughout the core and would depend on both the temperature (expansion) and the irradiation (shrinkage) of the graphite. The details of the water-cooled RCCS are not included. Instead, the heat sink provided by the RCCS is represented as a boundary condition. For this purpose, RCCS is assumed to have a constant, uniform temperature of 65 ◦ C. The radiation and convection in the air cavity between the outer reactor vessel and the RCCS are modeled by an effective heat transfer coefficient. To determine the heat transfer by radiation, the geometry is assumed to simplify to a pair of con-
Fig. 6. The grids representing the fluid parts (top) and solid parts (bottom) of the equivalent homogeneous model.
1928
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
Fig. 7. The temperature field for normal operating conditions, calculated using the explicit model.
centric surfaces. To determine the heat transfer by conduction and convection, an appropriate correlation is used. The method used to calculate the solid-fluid heat transfer varies across the model. In some locations, such as the flow inside the coolant channels, the heat transfer is calculated using a standard set of heat transfer correlations. In other locations, where the spatial resolution of the model is sufficient, it is calculated by the STAR-CD code using its standard algorithms. 4.3. Physical properties The physical properties of the steel components are assumed to not vary with temperature. The steel surfaces are all assumed to have an emissivity of 0.6. The helium coolant is modeled as an ideal gas. Its conductivity, molecular viscosity, and specific heat are all assumed to be temperature independent. The thermal conductivity of graphite depends in general on its temperature and irradiation history. For the results presented here, the graphite fuel blocks and reflectors are assumed to be irradiated to a fluence that is consistent with an end of life situation. Thus, its thermal conductivity is significantly lower than the conductivity of “fresh” non-irradiated graphite. Although the conductivity of graphite is generally non-isotropic, the model presented here uses an isotropic conductivity consisting of the lowest possible value. In addition, the physical annealing effect that is observed in graphite, whereby irradiation damage to the graphite is repaired at high temperatures and better conductivity values are recovered, is not included in the model. These simplifications result in higher peak fuel temperatures, and therefore, are conservative for fuel performance considerations. The graphite surfaces are assumed to have an emissivity of 0.9.
model (i.e., the two mesh components occupy the same geometrical space in the model). This arrangement simplifies meshing considerations, especially with respect to connecting this region of the model to the other parts for fluid or heat transfer. The fluid cells (representing the coolant channels) are connected to the upper and lower plenums, while the solid cells (representing the graphite) exchange heat with the surrounding material by conduction and radiation. The overlapping fluid and solid cells are coupled by calculating the heat flux from the fluid to the solid: φ = K(Tf − Ts ), where Tf and Ts are the temperatures of the corresponding fluid and solid cells. The coefficient K represents the heat transfer (by a combination of conduction, radiation, and convection) between the helium and the graphite. Its value is obtained from a separate detailed, explicit model of a single fuel block. To obtain the bulk homogenized equivalent heat transfer coefficient, an internal heat source P and inlet coolant temperature are prescribed, and the average temperatures of the helium Tf and the graphite web Ts are obtained from the computation. Thus, K = P/(Tf − Ts ).
4.4. The homogeneous model for the fuel blocks A distributed resistance model (porous medium) is used to represent the coolant channels in the fuel blocks. The graphite of the fuel blocks themselves is represented by a separate set of mesh cells. Because the solid parts and fluid parts of the model in this region are coupled, the two sets of cells are overlaid in the
Fig. 8. Comparison of the temperatures calculated using the homogeneous model and the explicit model.
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
1929
Fig. 10. The fluid cells in the core, which represent the coolant channels and the gaps between the graphite blocks.
determined from eq = φ/T. This calculation is performed for several (average) temperatures of the fuel block to obtain the equivalent conductivity as a function of temperature eq (T). To validate the equivalent homogeneous model, a comparison was made between the detailed fuel block model in Fig. 5 and its equivalent counterpart, shown in Fig. 6, for a transient representative of a pressurized conduction cooldown. Starting from a steady state in which a forced flow of helium through the channels removes the internally generated power (the thermal field calculated by the explicit model is shown in Fig. 7), a transient calculation was performed in which the forced flow was stopped, hot and cold temperatures were imposed on opposite vertical faces of the model, and residual power continued to be generated by the fuel compacts or their homogenized equivalent. The results, shown in Fig. 8, demonstrate a good agreement between the two models.
Fig. 9. The fluid cells (left) and solid cells (right) of the full computational domain.
Even if the conductivity of the graphite is assumed to be isotropic, the non-isotropic geometry of the fuel blocks means that the rate of heat transfer in the radial (horizontal) direction is not the same as the rate of heat transfer in the axial direction. Since the rate of heat transfer in the radial direction is the smaller value of the two, it is used to determine the equivalent conductivity of the homogenized fuel blocks. Using the detailed model of the fuel block, shown in Fig. 5, a temperature gradient T is imposed across the block. From the temperature field, the heat flux φ is calculated, and the equivalent conductivity is
Fig. 11. Horizontal cross-section of the model showing the regions of the core.
1930
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
Fig. 13. The coolant flow during normal operation. Fig. 12. The temperature of the solid structures during normal operation.
4.5. The mesh
4.6. Modeling of convection over long times
The mesh for the entire domain is shown in Fig. 9. This mesh consists of 104,000 cells, divided into 36,000 solid cells and 68,000 fluid cells. Fig. 10 shows the fluid cells in the core, which represent the gaps between the graphite blocks and the homogenized coolant channels. A cross-section of the different regions in the core from the inner reflectors to the pressure vessel is presented in Fig. 11. Although the homogeneous approach described in Section 4.4 is used to represent the fuel blocks, the gaps between the blocks are modeled explicitly, since a porous media model would not be precise enough for resolving the natural convection during PCC. Because of this explicit approach using a three-dimensional model, the axial, radial, and azimuthal flow in the core are determined, and the radiative heat transfer between the blocks is calculated exactly for each pair of blocks. Since the temperature distribution in the core and the pressure vessel are the main objectives of these calculations, a finer mesh is used in these regions. The other parts of the model – mostly the top and bottom parts of the reactor – are meshed more coarsely to provide better computational efficiency.
A significant challenge for using this model to simulate PCC is that solving the convection field in the transient regime requires very small time steps, whereas the total time of the transient is driven by the decay heat curve and the core graphite thermal inertia. The first criterion is linked to the Courant number, and the STAR-CD methodology allows time steps corresponding to Courant numbers in the range of 10–50. This yields time steps of 10−3 to 10−1 s. The second criterion leads to transients that require tens of hours of simulated time, with 100 h being typical. Therefore, this model uses a method to alternate small time steps, during which both the flow field and the temperature field are calculated, and large time steps, during which the flow field is frozen and only the heat transfer is solved. Since the natural convection flow develops quickly and then does not change much until the overall temperature distribution in the reactor changes, the velocity field can be considered to be approximately constant over a period that is small compared to the overall evolution of the transient, but large compared to the time steps used to calculate the convective flow.
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
1931
Fig. 14. Temperature distribution of the solid structures during DCC.
The initial temperature field is taken from the steady-state calculation for normal operation, and zero flow is used for the initial velocity field for the PCC transient (corresponding to a loss of forced convection). During the first 1000 s of simulated time, the flow field is calculated every time step – i.e., the full PISO algorithm is run and both the NavierStokes and the energy equations are solved. The main convective flow develops during the beginning of the transient, and small time steps are needed to accurately capture the correct flow configuration. This process is complicated, however, by the coupling with radiative heat transfers. Specifically, a rapid increase in temperature occurs locally in some parts of the
model during the first few time steps, which changes the surface temperatures and hence greatly affects the radiative heat transfer, since the temperatures are already relatively high. After this initial period of the calculation, the global flow pattern stabilizes and evolves relatively slowly due to changes in the temperature distribution in the core, and less frequent updates of the flow field are sufficient to follow the evolution of this flow. Thus, the flow field is solved for five consecutive time steps (each lasting 10−3 to 10−1 s, depending on the Courant number) and then is frozen for 100 time steps of 15 s each. Such a procedure allows a transient lasting 100 h of simulated time
1932
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
Fig. 15. Temperature distribution of the solid structures during PCC.
to be computed in a reasonable amount of time (approximately one month of CPU time on a modern workstation). To test this scheme, two calculations were performed for 10 h of simulated time. For one calculation, the flow field was determined at every time step by solving the full set of coupled equations; for the other calculation, the approach described above was used. Comparison of the temperature fields for these two calculations showed very little difference in the results, with the maximum local difference being less than 5 ◦ C at the end of the 10 h transient. 5. Results For each of the calculations presented here, the model is first run at nominal operating conditions in a steady-state regime, and the resulting solution is used as the starting point of the transient calculation. Figs. 12 and 13 show the temperature and velocity maps in the reactor during normal operation. The core support structure, the core barrel, the duct shell, and the upper structures
are cooled at the inlet temperature of 490 ◦ C. The pressure vessel is partially insulated from this inlet flow by an annulus containing helium undergoing natural convection, so the vessel remains below 400 ◦ C. Although the flow reaches speeds of up to 70 m/s in some places, the fluid can still be considered incompressible because it is traveling much slower than the speed of sound (at least 1500 m/s under these conditions). The general flow pattern remains quasi-two-dimensional (axisymmetric), except for the azimuthal loops that form in the natural convection in the outer annulus. The reflector blocks (both inner and outer) are heated by the radial transfer of thermal energy from the active core and are cooled by the 10% of the helium flow that bypasses the coolant channels and flows in the gaps between the blocks. The “hot” regions (i.e., with temperatures T > 500 ◦ C) are found only in the active core annulus and in the outlet plenum. It is expected that future design studies will need to include a set of parametric studies for this flow in the gaps, which varies the amount of bypass from the no-bypass limit to 30% of the total flow.
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
1933
Fig. 16. Temperature distribution of the helium fluid during PCC.
The main objective of this paper is to study the event in which a loss of the forced flow occurs while the system remains pressured (pressurized conduction cooldown, PCC). However, to assess the effect of natural convection, the loss of flow event that is accompanied by depressurization of the system (depressurized conduction cooldown, DCC) is also computed. Comparison of the two scenarios will highlight the effects of the natural convection inside the reactor vessel. 5.1. Depressurized conduction cooldown This situation has been studied previously in detail with twodimensional models for exploratory purposes (Mays et al., 2004; Haque et al., 2004). The results of the three-dimensional model presented here are not directly comparable to the earlier results of the two-dimensional models, however, because of several differences in the design and the proposed operating conditions.
To understand the phenomena governing the passive removal of decay heat, the evolution of the temperature pattern is followed, including the migration of the hot spot from the bottom of the active core and its effect on the metallic structures. Fig. 14 presents the time history of the temperature of the solid structures. The hot region moves from the bottom of the core to the center due to equilibrium heat losses. The upper and lower parts are cooled, and the core maximum temperature (determined with the current postulated operating conditions to be 1450 ◦ C) occurs at around 80 h after the shutdown. During the transient, the internal power decreases. The heat losses of the vessel increase in the central part, but decrease in the upper and lower parts, so that the total losses do not change much during the transient. After approximately 80 h, the heat losses become greater than the internal power and the reactor globally cools down slowly. The annealing effect of the graphite, which has not been taken into account
1934
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
in these calculations, could reduce slightly the peak fuel temperatures. At the beginning of the transient, the central part of the pressure vessel (the cylindrical part) is first cooled because of the heat removal to the RCCS circuit and because the convective heat transfer from helium at the inlet temperature of 490 ◦ C vanishes. After 10 h, the heat from the core produces a rise in the pressure vessel temperature, and the maximum is reached at nearly 480 ◦ C a few hours after the end of the 100-h transient. 5.2. Pressurized conduction cooldown The time history of temperature during the PCC transient is shown in Fig. 15 for the solid domain and in Fig. 16 for the fluid domain. Natural convection develops quickly in the initial times. The flow develops first in the active core channels: helium flows upward in the central row and downward in the inner and outer rows. After about 10 h of simulated time, the flow develops into a different configuration, with helium flowing down in the channels of the outer row of blocks of the active core and in the gaps between the outer reflector blocks. It rises in the two inner rows of the active core and in the gaps between the inner reflector blocks. This flow pattern is diagrammed in Fig. 17. The hot spot, located during normal operation at the bottom part of the active core, moves slowly to the upper quarter of the core, as the axial temperature profiles in Fig. 18 demonstrate. The maximum core temperature decreases at first, since the hot spot moves to colder regions until about 10 h after the beginning of the transient due to the development of natural convection in the core. After this, the situation remains the same during the rest of the transient: the decay heat is removed via conduction and radiation to the pressure vessel, and also by the convection field in both the channels and the gaps, which is shown in Fig. 19. The peak fuel temperatures (at the mesh scale) are not found to be much higher during PCC than during normal operation. Nevertheless, the more important consequence is the redistribution of temperature during the event (shown axially in Fig. 18) and the large temperature increase in the fuel in the upper part of the reactor. Moreover, the temperatures computed are mesh-averaged temperatures. The actual peak temperature of the fuel at the center of the compact is almost identical to the mesh-averaged temperature after the reactor trip, but this peak temperature is about 100 ◦ C hotter than then mesh-averaged temperature during full-power operation. The lower structures cool while the upper plenum and its neighboring structures heat up during the first 50 h of the transient. The temperature of the helium in the upper plenum rises from 950 to 1000 ◦ C after 50 h. In the riser, the flow develops into radial and azimuthal loops. This situation is not very stable and may be linked to the choice of a 30◦ sector as the computational domain; a full 360◦ model, or at least a 180◦ model of half of the reactor, would be necessary for a precise study of this region. Nevertheless, the power transferred across the riser is correctly assessed. The core barrel reaches a temperature of 750 ◦ C after 80 h, as shown in Fig. 20, with
Fig. 17. The helium flow pattern during PCC after 10 h of simulated time.
this maximum temperature located at the very upper part of the barrel. Azimuthal convection loops develop in the annuli due to the narrowness of these regions and the three-dimensional dissymmetries. Some of these convection loops are visible in Fig. 21, which presents the velocity field in the annular space of the riser at the end of the 100 h PCC transient. Figs. 22 and 23 present the time evolution of the maximum temperatures of the core and the pressure vessel for both the PCC and DCC scenarios. The DCC is the more severe scenario for the maximum temperature. Nevertheless, the PCC is more penalizing for the upper structures, including the upper part of the pressure vessels (e.g., the flanges where the vessel head attaches). During the PCC, the maximum temperatures occur between 60 h (for the core) and 80 h (for the vessel), which corresponds to the time when the maximum amount of heat is stored
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
1935
Fig. 18. Axial temperature profiles for selected times during PCC.
in the reactor, as shown in Fig. 24. After this period, the power removed from the system is greater than the power produced by the decay heat. Finally, the PCC behavior is, after the first 10 h of the transient, rather similar to a DCC scenario with a different initial temperature distribution inside the reactor (i.e., the temperature pattern is shifted toward the top of the core). Fig. 22
illustrates the two periods of the transient: before 10 h, when the hot spot moves upward due to convection, and after 10 h, which is similar to a DCC transient. The velocity field (Fig. 19) shows the characteristic shape of the convection inside the gaps. This convection acts globally, after it has been established, to increase the radial heat transfer and to vertically shift hot spots toward the top of the vessel. The convective field
Fig. 19. The helium flow in the gaps and the core channels during PCC at 100 h of simulated time.
1936
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
Fig. 20. The maximum temperature of the core barrel during the PCC transient.
Fig. 22. The maximum temperature of the core during the PCC transient (solid line) and during the DCC transient (dashed line).
bounds the energy stored in the early part of the transient by increasing the radial losses, but it limits the peak fuel temperature mainly because it moves the hot spot to previously “cold” areas.
Fig. 23. The maximum temperature of the pressure vessel during the PCC transient (solid line) and during the DCC transient (dashed line).
Fig. 24. The energy stored in the core during the PCC transient.
6. Conclusions
Fig. 21. The flow field in the riser annulus during PCC at 100 h of simulated time.
An important safety characteristic of modern HTRs is passive heat removal in case of a loss of forced convection under either pressurized or depressurized conditions. The depressurized case (DCC) is a classical situation that has been examined by many investigators. The pressurized case (PCC), when the primary loop remains pressurized, has not been investigated as widely, however. A three-dimensional model, using the STAR-CD software, has been built. It solves the combined conductive, convective, and radiative heat transfer using an equivalent homogeneous media for transient events with a duration of 100 h or more.
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
Special techniques were employed to simulate the internal natural circulation flows during PCC without slowing the simulation excessively. These are scoping calculations, which do not include the variability of properties or design margins. The objectives of these calculations were both to set up a methodology and to explore the general behavior of HTRs in such situations. Compared to the depressurized situation, the extra convective heat transfer reduces the peak temperature of both the fuel and the pressure vessel during the PCC event. However, it leads to thermal stratification and hence to much higher temperatures in the upper reactor structures. The next planned extension of this model will be to couple it with a detailed model of convection inside the reactor cavity, in order to include the effects of thermal stratification in the surrounding air. The second planned extension will be to explicitly include the remainder of the primary loop outside the reactor vessel. Concerning comparisons with experimental data, separate effects data are available for individual phenomena, but integral test data of the entire facility are not currently available. The need for integral testing is still being evaluated in the industry (given that the phenomena are simpler and fewer in number than for LWR LOCA analysis). Additional separate effect test-
1937
ing could be of value. The greatest benefit would likely come from additional testing of material properties. References Computational Dynamics Limited, 2004. STAR-CD Version 3.24 User Guide. Computational Dynamics Limited, London. Haque, H., Feltes, W., Brinkmann, G., 2004. Thermal response of a high temperature reactor during passive cooling under pressurized and depressurized conditions. In: Proceedings of the Second International Topical Meeting on High Temperature Reactor Technology, Paper No. F02, Beijing, China, September 22–24. Issa, R.I., Ahmadi Befui, B., Beshay, K., Gosman, A.D., 1991. Solution of the implicit discretised reacting flow equations by operator splitting. J. Comp. Phys. 93, 388–410. Mays, B.E., Woaye-Hune, A., Simoneau, J.-P., Gabeloteau, T., Lefort, F., Haque, H., Lommers, L., 2004. The effect of operating temperature on depressurized conduction cooldown for a high temperature reactor. In: Proceedings of the ICAPP’04, Paper No. 4202, Pittsburgh, PA, USA, June 13–17. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation, New York. Petit, D., 2005. Overall simulation of a HTGR plant with the gas adapted MANTA code. In: Proceedings of the NURETH-11 Conference, Avignon, France, October 2–6. Simoneau, J.-P., Champigny, J., 2004. Large eddy simulation of mixing in the outlet plenum of a high temperature reactor: a benchmark exercise. In: Proceedings of the ICONE 12 Conference, Washington, DC, USA, April 25–29.
Three-dimensional simulation of the coupled convective, conductive, and radiative heat transfer during decay heat removal in an HTR Jan-Patrice Simoneau a,∗ , Julien Champigny a , Brian Mays a , Lewis Lommers b a
AREVA, 10 rue Juliette R´ecamier, 69456 Lyon Cedex 06, France b AREVA, 3315 Old Forest Road, Lynchburg, VA 24506, USA
Received 16 June 2006; received in revised form 14 March 2007; accepted 14 March 2007
Abstract A three-dimensional model has been constructed to simulate the passive heat removal in a modular prismatic-block high temperature reactor during a loss of active cooling accident. This model, developed using the STAR-CD general computational fluid dynamics code, solves the combined conductive, convective, and radiative heat transfer within a 30◦ section of the core and reactor vessel. To accommodate the different spatial scales, it uses homogeneous equivalent media to represent the coolant flow and the prismatic fuel blocks. A customized procedure that manages solving alternatively the dynamic and thermal fields permits the computation of very long transients, which typically are performed for 100 or more hours of simulated time. The global methodology and specific modeling procedures are explained, and key points of the CFD analysis are highlighted. Next, the results of several calculations are presented, and the physical phenomena represented are described. Two commonly investigated loss of active cooling scenarios are considered: depressurized conduction cooldown and pressurized conduction cooldown. The results for these two scenarios are compared to assess the effect of heat transfer via internal natural convection – which is negligible during the depressurized event – on the thermal behavior of the system. In addition, the evolution of the natural convection flow through the core and in the annular spaces is examined and discussed. © 2007 Elsevier B.V. All rights reserved.
1. Introduction Recent years have seen a renewed interest in the high temperature reactor (HTR) as a high-temperature nuclear heat source. This interest is driven largely by the possibility of using nuclear energy for the production of hydrogen. The most efficient means of hydrogen production, such as high-temperature electrolysis or thermochemical water splitting, require very high temperatures, which are only achieved by HTRs. Since greater efficiency is possible as the process temperature is increased, current HTR designs are being extended to produce even higher outlet temperatures, increasing from 850 ◦ C to temperatures approaching 1000 ◦ C. The 600 MW HTR considered here is a prismatic-block concept. The power conversion system is an indirect combined Brayton and Rankine cycle. This HTR design is capable of pas-
∗
Corresponding author. E-mail address: [email protected] (J.-P. Simoneau).
0029-5493/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2007.03.010
sive decay heat removal for enhanced inherent safety. Decay heat can be removed using only conduction and radiation from the fuel to the vessel and radiation and natural convection from the vessel to a water circuit located on the vault walls. The reactor is configured to ensure acceptable temperatures for the fuel particles under these conditions. In addition to conduction and thermal radiation phenomena, natural convection can also occur within the core during passive cooling. However, the convective heat transfer remains negligible if the circuit is depressurized (referred to as depressurized conduction cooldown, or DCC). On the contrary, if the circuit remains pressurized (referred to as pressurized conduction cooldown, or PCC), the internal natural convection provides noticeable heat transfer, and a different thermal response is observed in the core structures. Thus, PCC implies combined convective/conductive/radiative heat transfer. The combination of these heat transfer modes, each with very different dynamical characteristics, can pose a computational challenge. The objective of this paper is to present the implementation of a 3-D PCC computation and the resulting thermal response of
1924
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
the reactor. The computations presented are general exploratory calculations, not design calculations. 2. Key features of the reactor design 2.1. General characteristics The HTR concept considered here uses three circuits, operating on a combined Brayton-Rankine cycle (Petit, 2005). The primary circuit, consisting of pressurized helium, transfers the heat from the reactor core to an intermediate heat exchanger (IHX). The secondary circuit, working with a nitrogen-based fluid, drives a gas turbine to form a Brayton cycle that uses a turbo compressor coupled with an alternator. The final bottoming circuit is a water-steam system that works as a conventional Rankine cycle. The key parameters of the reactor are given in Table 1. The core consists of hexagonal graphite blocks, which are divided into an active core of 1020 hexagonal prismatic fuel elements, stacked in 102 columns of 10 blocks each, and a set of solid graphite reflector blocks. The active core is arranged in an annulus, with some of the reflector blocks (the inner reflectors) located inside the annulus and the remaining blocks (the outer reflectors) located between the active core and the reactor vessel. The core is held by a steel core barrel, and all structures are enclosed in a pressure vessel, which is connected to an adjacent vessel containing the IHX via a short cross-vessel. Figs. 1 and 2 present axial and radial cross-sections of the reactor vessel and core. 2.2. Description of the primary coolant flow Under normal operating conditions, cold helium travels through the outer annulus of the cross-vessel and enters the bottom part of the reactor vessel, where it cools the metallic support structures. It then flows up to the top of the core in a riser between the core and the pressure vessel, to limit heat losses and protect the pressure vessel from excessive temperatures. Starting from a plenum at the top of the vessel, the coolant flows downward through the core, where it is heated, and is collected in a second plenum below the core, before exiting through the inner annulus of the cross-vessel to the IHX vessel. Most of the helium flows through the coolant channels in the blocks that form the active core, but part of this flow bypasses
Fig. 1. Vertical cross-section of the reactor indicating the flow path of the helium coolant during normal operation.
Table 1 Key reactor parameters Parameter
Value
Core thermal power (MW) Helium mass flow rate (kg/s) Helium pressure (bar) Reactor inlet temperature (◦ C) Reactor outlet temperature (◦ C) Core inner diameter (m) Core outer diameter (m) Core height (m) Vessel outer diameter (m)
600 226 50 490 1000 2.96 4.84 8 7.7
Fig. 2. Horizontal cross-section of the reactor.
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
the channels by flowing through the gaps between the blocks. This bypass flow cools the reflectors, but does not contribute directly to the nuclear heat removal. It mixes with the main flow in the lower plenum before exiting the reactor vessel. Very complex flows (Simoneau and Champigny, 2004) occur in this outlet plenum due to three effects: injection of “hot” helium from the coolant channels, injection of colder helium from the gaps, and highly turbulent flow that forms in the wakes resulting from the presence of graphite columns, which support the core, in the transverse flow toward the plenum exit. The gaps between the blocks play a significant role in the flow distribution and hence the thermal behavior of the reactor. During normal operation, the bypass flow has undesired consequences, since it does not directly cool the fuel. Therefore, the temperatures at the coolant channel outlets are locally higher than the reactor outlet temperature due to this bypass. Decay heat removal without any primary flow, on the other hand, is a transient situation where the thermal inertia of the different regions of the reactor dominate. During the beginning of this transient, the reflectors act as a heat sink, which limits the maximum temperatures and protects the fuel blocks. The temperature distribution within the core at the beginning of the transient – and hence the amount of heat that the reflectors can absorb – is mainly a result of the amount of gap flow. For this study, the gaps between the graphite blocks are assumed to be uniformly 2 mm wide, the value for new blocks at ambient temperature. During reactor operation, however, the blocks experience size changes due to the effects of thermal expansion, which reduces the gap size, and irradiation of the graphite, which tends to shrink the graphite, increasing the gap size. Thus, the bypass flow behavior changes greatly from fresh to used blocks. The loss of primary flow with fresh fuel and a core of new graphite blocks leads to higher reflector temperatures and at the same time to lower residual power. In the calculations presented in this paper, a bypass flow of 10% was selected and prescribed for the model; however, other calculations (not reported here) with different bypass flow values have been performed. 2.3. The fuel blocks The fuel blocks are hexagonal right-regular prisms composed of graphite with holes that form the coolant channels and smaller holes into which the fuel compacts are inserted. The TRISOcoated spherical fuel particles are incorporated with a graphite matrix to form the cylindrical compacts. These components are depicted in Fig. 3. 2.4. The decay heat removal system In the event of a loss of primary flow, the decay heat can be removed by the reactor cavity cooling system (RCCS). This system consists of vertical water pipes, which are set on the concrete walls of the cavity surrounding the reactor vessel. These pipes, which are connected to a tank, are designed to maintain an average temperature of 65 ◦ C and to protect the concrete walls during normal operation.
1925
Fig. 3. The prismatic fuel block, compact, and TRISO particle.
The system described here spans several very different spatial scales, which a CFD calculation has to manage. In terms of scale, the solid parts range from the fuel particles (less than 1 mm in size) to the reactor vessel (several meters long), and the fluid volumes range from the channel diameter (16 mm) and the gap width (2 mm) to the cavity height (about 20 m). 3. Assumptions for decay heat removal transients The transient operating conditions considered here for a loss of primary flow are as follows. The initiating event is a loss of the primary forced flow. In addition, all heat removal mechanisms except for the reactor cavity cooling system are assumed to be unavailable. Two scenarios are then considered. In the first scenario, depressurized conduction cooldown (DCC), the primary circuit is depressurized, either via a large break or intentional depressurization. (The term “conduction” is adopted here because the behavior of the reactor is similar to the cooling of a solid body by conduction; nevertheless, radiation and, to a much more limited extent, convection both contribute to the heat removal process.) Since the density of the helium gas is low – on the order of atmospheric pressure, due to equilibrium with the air in the reactor cavity volume – the heat transfer by convection is relatively small compared to other phenomena. The other scenario, pressurized conduction cooldown (PCC), assumes that the coolant remains in the primary circuit. Thus, natural convection is much more important than in the previous case.
1926
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
The phenomenon of decay heat removal occurs on very large time scales (tens of hours), but the development of convection in the core channels occurs within a few minutes. Moreover, the characteristic times of the convection in some regions inside the reactor are smaller than one second (e.g., T ∼ L/U, with U ∼ 0.5 m/s and L ∼ 10 mm). Thus, to simulate such a system, a model must manage a wide variety of time scales. 4. The numerical model and code 4.1. Computational techniques
Fig. 4. The decay heat curve (main plot) and axial distribution of the power (inset).
The conditions in the DCC scenario result in higher fuel temperatures, since the transfer of thermal energy from the core to the heat sink is limited to the radiative and conductive processes in the blocks and the gaps between the components. The PCC scenario can be the more limiting case for the pressure vessel, since the heat transfer from the core to the vessel is increased by the natural convection. Furthermore, this situation can lead to hot spots in the upper part of the reactor. Both scenarios follow the same general timeline. Before the beginning of the transient, the reactor is operating at full power (600 MW), the RCCS is in operation, and 10% of the total coolant flow bypasses the active core. The power generated in the annular core is assumed to be distributed uniformly in the radial direction. The axial power profile is shown in Fig. 4. At the beginning of the transient (time t = 0 s), the primary helium circulation is stopped, and the protection system automatically drops the rods in response to the loss of main loop circulation. In the DCC scenario, the system also rapidly depressurizes, so that it can be considered completely depressurized when the transient begins. During the transient, the decay heat is dissipated in the core, while the RCCS continues to operate and serves as the ultimate heat sink for the system. Beginning at 30 MW (5% of full power), the decay heat decreases according to the curve in Fig. 4, which is a standard radioactive decrease for a core at end of cycle. The axial power profile is assumed to remain homothetic to its initial profile. After the loss of forced flow, natural circulation develops in the reactor vessel as a result of the buoyancy forces arising from the temperature distribution at normal operation, which serves as the initial condition of the transient. This natural convection occurs through the core channels and gaps, in the top collector, in the bottom collector, and in the annular spaces along the pressure vessel. A shutoff valve, which automatically closes when forced convection is lost, prevents circulation through the primary loop. Natural convection occurs in both the DCC and PCC scenarios, but as mentioned above, the effect of this convection is ignored in the DCC case, since is it very weak and the amount of heat transfered by it is negligible.
The general-purpose fluid-mechanics and heat-transfer code STAR-CD, version 3.2, (Computational Dynamics Limited, 2004) was used to simulate the decay heat removal transients. This code solves the equations of fluid flow and heat transfer on an unstructured mesh, which for this study consists mainly of hexahedral and prismatic cell volumes. The partial differential equations are discretized by the finite volume method, using first-order upwind differences for the spatial terms and a firstorder, fully implicit scheme for the temporal discretization. The matrices that appear in the solution of the discrete finite volume equations are inverted by the double conjugate gradient method with preconditioning provided by the algebraic multigrid method. This approach has proved to be very useful for solving the pressure correction equation in tall geometries, such as those that exist in the model described here. Two algorithms are used to solve the pressure linked equations. For the steady-state calculations that are used to provide the initial conditions for the transient, the classic SIMPLE algorithm (Patankar, 1980) is used. During the transient, the more sophisticated PISO algorithm (Issa et al., 1991) is used. The convergence of the latter algorithm is ensured by selecting sufficiently small time steps. Turbulence is modeled only during the steady-state calculation for normal operation, which has high flow velocities due to the forced circulation. In this case, the model uses STAR-CD’s standard linear k–ε model and wall functions for high Reynolds number flow, which are well suited for forced convection. Heat transfer by thermal radiation is modeled by face-to-face exchanges of energy using a discrete beam method to calculate the view factors. This approach assumes that all surfaces are gray and diffuse and that the helium gas does not participate in the radiative heat transfer. The equations that are solved for each type of calculation described here are summarized in Table 2. Table 2 Equations solved by the STAR-CD code Physical quantity
Initial steady-state
PCC transient
DCC transient
Flow field
Solved with SIMPLE Solved Solved with k–ε high Re Solved
Solved with PISO Solved Not solved (Laminar flow) Solved
Not solved
Temperature field Turbulence Density
Solved Not solved Not solved
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
4.2. Modeling philosophy Instead of explicitly modeling the entire reactor core, only the minimum section required to represent the core threedimensionally via symmetry – a 30◦ sector – was included in the model. Although this sector is sufficient to represent the core itself, a rigorous modeling of the outlet plenum would require at least a 180◦ sector. Therefore, the modeling approach used here leads to an approximation of the outlet flow; however, a detailed representation of this flow is not the objective of this study. The metallic structures of the reactor are explicitly modeled, but some details are neglected: the shroud tubes (for control rods) in the upper plenum, the vertical ribs of the lower supporting plates, and other very small parts. Similarly, the graphite posts in the outlet plenum are not modeled. The blocks are represented by a homogeneous equivalent solid medium, which includes both the fuel compacts and the
Fig. 5. The grids representing the fluid parts (top) and solid parts (bottom) of the explicit model.
1927
graphite web. This technique is described in further detail in Section 4.4. The gaps between the core blocks are explicitly modeled by mesh cells, with the following assumptions. The solid-fluid heat transfer and the pressure drop across the core are prescribed via coefficients (i.e., the boundary layer is not resolved). The radiation across gaps is explicitly modeled. Finally, the gap size is assumed to be uniform and equal to the “cold” value of 2 mm. As mentioned above, the space between blocks in an actual core would be expected to vary throughout the core and would depend on both the temperature (expansion) and the irradiation (shrinkage) of the graphite. The details of the water-cooled RCCS are not included. Instead, the heat sink provided by the RCCS is represented as a boundary condition. For this purpose, RCCS is assumed to have a constant, uniform temperature of 65 ◦ C. The radiation and convection in the air cavity between the outer reactor vessel and the RCCS are modeled by an effective heat transfer coefficient. To determine the heat transfer by radiation, the geometry is assumed to simplify to a pair of con-
Fig. 6. The grids representing the fluid parts (top) and solid parts (bottom) of the equivalent homogeneous model.
1928
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
Fig. 7. The temperature field for normal operating conditions, calculated using the explicit model.
centric surfaces. To determine the heat transfer by conduction and convection, an appropriate correlation is used. The method used to calculate the solid-fluid heat transfer varies across the model. In some locations, such as the flow inside the coolant channels, the heat transfer is calculated using a standard set of heat transfer correlations. In other locations, where the spatial resolution of the model is sufficient, it is calculated by the STAR-CD code using its standard algorithms. 4.3. Physical properties The physical properties of the steel components are assumed to not vary with temperature. The steel surfaces are all assumed to have an emissivity of 0.6. The helium coolant is modeled as an ideal gas. Its conductivity, molecular viscosity, and specific heat are all assumed to be temperature independent. The thermal conductivity of graphite depends in general on its temperature and irradiation history. For the results presented here, the graphite fuel blocks and reflectors are assumed to be irradiated to a fluence that is consistent with an end of life situation. Thus, its thermal conductivity is significantly lower than the conductivity of “fresh” non-irradiated graphite. Although the conductivity of graphite is generally non-isotropic, the model presented here uses an isotropic conductivity consisting of the lowest possible value. In addition, the physical annealing effect that is observed in graphite, whereby irradiation damage to the graphite is repaired at high temperatures and better conductivity values are recovered, is not included in the model. These simplifications result in higher peak fuel temperatures, and therefore, are conservative for fuel performance considerations. The graphite surfaces are assumed to have an emissivity of 0.9.
model (i.e., the two mesh components occupy the same geometrical space in the model). This arrangement simplifies meshing considerations, especially with respect to connecting this region of the model to the other parts for fluid or heat transfer. The fluid cells (representing the coolant channels) are connected to the upper and lower plenums, while the solid cells (representing the graphite) exchange heat with the surrounding material by conduction and radiation. The overlapping fluid and solid cells are coupled by calculating the heat flux from the fluid to the solid: φ = K(Tf − Ts ), where Tf and Ts are the temperatures of the corresponding fluid and solid cells. The coefficient K represents the heat transfer (by a combination of conduction, radiation, and convection) between the helium and the graphite. Its value is obtained from a separate detailed, explicit model of a single fuel block. To obtain the bulk homogenized equivalent heat transfer coefficient, an internal heat source P and inlet coolant temperature are prescribed, and the average temperatures of the helium Tf and the graphite web Ts are obtained from the computation. Thus, K = P/(Tf − Ts ).
4.4. The homogeneous model for the fuel blocks A distributed resistance model (porous medium) is used to represent the coolant channels in the fuel blocks. The graphite of the fuel blocks themselves is represented by a separate set of mesh cells. Because the solid parts and fluid parts of the model in this region are coupled, the two sets of cells are overlaid in the
Fig. 8. Comparison of the temperatures calculated using the homogeneous model and the explicit model.
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
1929
Fig. 10. The fluid cells in the core, which represent the coolant channels and the gaps between the graphite blocks.
determined from eq = φ/T. This calculation is performed for several (average) temperatures of the fuel block to obtain the equivalent conductivity as a function of temperature eq (T). To validate the equivalent homogeneous model, a comparison was made between the detailed fuel block model in Fig. 5 and its equivalent counterpart, shown in Fig. 6, for a transient representative of a pressurized conduction cooldown. Starting from a steady state in which a forced flow of helium through the channels removes the internally generated power (the thermal field calculated by the explicit model is shown in Fig. 7), a transient calculation was performed in which the forced flow was stopped, hot and cold temperatures were imposed on opposite vertical faces of the model, and residual power continued to be generated by the fuel compacts or their homogenized equivalent. The results, shown in Fig. 8, demonstrate a good agreement between the two models.
Fig. 9. The fluid cells (left) and solid cells (right) of the full computational domain.
Even if the conductivity of the graphite is assumed to be isotropic, the non-isotropic geometry of the fuel blocks means that the rate of heat transfer in the radial (horizontal) direction is not the same as the rate of heat transfer in the axial direction. Since the rate of heat transfer in the radial direction is the smaller value of the two, it is used to determine the equivalent conductivity of the homogenized fuel blocks. Using the detailed model of the fuel block, shown in Fig. 5, a temperature gradient T is imposed across the block. From the temperature field, the heat flux φ is calculated, and the equivalent conductivity is
Fig. 11. Horizontal cross-section of the model showing the regions of the core.
1930
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
Fig. 13. The coolant flow during normal operation. Fig. 12. The temperature of the solid structures during normal operation.
4.5. The mesh
4.6. Modeling of convection over long times
The mesh for the entire domain is shown in Fig. 9. This mesh consists of 104,000 cells, divided into 36,000 solid cells and 68,000 fluid cells. Fig. 10 shows the fluid cells in the core, which represent the gaps between the graphite blocks and the homogenized coolant channels. A cross-section of the different regions in the core from the inner reflectors to the pressure vessel is presented in Fig. 11. Although the homogeneous approach described in Section 4.4 is used to represent the fuel blocks, the gaps between the blocks are modeled explicitly, since a porous media model would not be precise enough for resolving the natural convection during PCC. Because of this explicit approach using a three-dimensional model, the axial, radial, and azimuthal flow in the core are determined, and the radiative heat transfer between the blocks is calculated exactly for each pair of blocks. Since the temperature distribution in the core and the pressure vessel are the main objectives of these calculations, a finer mesh is used in these regions. The other parts of the model – mostly the top and bottom parts of the reactor – are meshed more coarsely to provide better computational efficiency.
A significant challenge for using this model to simulate PCC is that solving the convection field in the transient regime requires very small time steps, whereas the total time of the transient is driven by the decay heat curve and the core graphite thermal inertia. The first criterion is linked to the Courant number, and the STAR-CD methodology allows time steps corresponding to Courant numbers in the range of 10–50. This yields time steps of 10−3 to 10−1 s. The second criterion leads to transients that require tens of hours of simulated time, with 100 h being typical. Therefore, this model uses a method to alternate small time steps, during which both the flow field and the temperature field are calculated, and large time steps, during which the flow field is frozen and only the heat transfer is solved. Since the natural convection flow develops quickly and then does not change much until the overall temperature distribution in the reactor changes, the velocity field can be considered to be approximately constant over a period that is small compared to the overall evolution of the transient, but large compared to the time steps used to calculate the convective flow.
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
1931
Fig. 14. Temperature distribution of the solid structures during DCC.
The initial temperature field is taken from the steady-state calculation for normal operation, and zero flow is used for the initial velocity field for the PCC transient (corresponding to a loss of forced convection). During the first 1000 s of simulated time, the flow field is calculated every time step – i.e., the full PISO algorithm is run and both the NavierStokes and the energy equations are solved. The main convective flow develops during the beginning of the transient, and small time steps are needed to accurately capture the correct flow configuration. This process is complicated, however, by the coupling with radiative heat transfers. Specifically, a rapid increase in temperature occurs locally in some parts of the
model during the first few time steps, which changes the surface temperatures and hence greatly affects the radiative heat transfer, since the temperatures are already relatively high. After this initial period of the calculation, the global flow pattern stabilizes and evolves relatively slowly due to changes in the temperature distribution in the core, and less frequent updates of the flow field are sufficient to follow the evolution of this flow. Thus, the flow field is solved for five consecutive time steps (each lasting 10−3 to 10−1 s, depending on the Courant number) and then is frozen for 100 time steps of 15 s each. Such a procedure allows a transient lasting 100 h of simulated time
1932
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
Fig. 15. Temperature distribution of the solid structures during PCC.
to be computed in a reasonable amount of time (approximately one month of CPU time on a modern workstation). To test this scheme, two calculations were performed for 10 h of simulated time. For one calculation, the flow field was determined at every time step by solving the full set of coupled equations; for the other calculation, the approach described above was used. Comparison of the temperature fields for these two calculations showed very little difference in the results, with the maximum local difference being less than 5 ◦ C at the end of the 10 h transient. 5. Results For each of the calculations presented here, the model is first run at nominal operating conditions in a steady-state regime, and the resulting solution is used as the starting point of the transient calculation. Figs. 12 and 13 show the temperature and velocity maps in the reactor during normal operation. The core support structure, the core barrel, the duct shell, and the upper structures
are cooled at the inlet temperature of 490 ◦ C. The pressure vessel is partially insulated from this inlet flow by an annulus containing helium undergoing natural convection, so the vessel remains below 400 ◦ C. Although the flow reaches speeds of up to 70 m/s in some places, the fluid can still be considered incompressible because it is traveling much slower than the speed of sound (at least 1500 m/s under these conditions). The general flow pattern remains quasi-two-dimensional (axisymmetric), except for the azimuthal loops that form in the natural convection in the outer annulus. The reflector blocks (both inner and outer) are heated by the radial transfer of thermal energy from the active core and are cooled by the 10% of the helium flow that bypasses the coolant channels and flows in the gaps between the blocks. The “hot” regions (i.e., with temperatures T > 500 ◦ C) are found only in the active core annulus and in the outlet plenum. It is expected that future design studies will need to include a set of parametric studies for this flow in the gaps, which varies the amount of bypass from the no-bypass limit to 30% of the total flow.
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
1933
Fig. 16. Temperature distribution of the helium fluid during PCC.
The main objective of this paper is to study the event in which a loss of the forced flow occurs while the system remains pressured (pressurized conduction cooldown, PCC). However, to assess the effect of natural convection, the loss of flow event that is accompanied by depressurization of the system (depressurized conduction cooldown, DCC) is also computed. Comparison of the two scenarios will highlight the effects of the natural convection inside the reactor vessel. 5.1. Depressurized conduction cooldown This situation has been studied previously in detail with twodimensional models for exploratory purposes (Mays et al., 2004; Haque et al., 2004). The results of the three-dimensional model presented here are not directly comparable to the earlier results of the two-dimensional models, however, because of several differences in the design and the proposed operating conditions.
To understand the phenomena governing the passive removal of decay heat, the evolution of the temperature pattern is followed, including the migration of the hot spot from the bottom of the active core and its effect on the metallic structures. Fig. 14 presents the time history of the temperature of the solid structures. The hot region moves from the bottom of the core to the center due to equilibrium heat losses. The upper and lower parts are cooled, and the core maximum temperature (determined with the current postulated operating conditions to be 1450 ◦ C) occurs at around 80 h after the shutdown. During the transient, the internal power decreases. The heat losses of the vessel increase in the central part, but decrease in the upper and lower parts, so that the total losses do not change much during the transient. After approximately 80 h, the heat losses become greater than the internal power and the reactor globally cools down slowly. The annealing effect of the graphite, which has not been taken into account
1934
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
in these calculations, could reduce slightly the peak fuel temperatures. At the beginning of the transient, the central part of the pressure vessel (the cylindrical part) is first cooled because of the heat removal to the RCCS circuit and because the convective heat transfer from helium at the inlet temperature of 490 ◦ C vanishes. After 10 h, the heat from the core produces a rise in the pressure vessel temperature, and the maximum is reached at nearly 480 ◦ C a few hours after the end of the 100-h transient. 5.2. Pressurized conduction cooldown The time history of temperature during the PCC transient is shown in Fig. 15 for the solid domain and in Fig. 16 for the fluid domain. Natural convection develops quickly in the initial times. The flow develops first in the active core channels: helium flows upward in the central row and downward in the inner and outer rows. After about 10 h of simulated time, the flow develops into a different configuration, with helium flowing down in the channels of the outer row of blocks of the active core and in the gaps between the outer reflector blocks. It rises in the two inner rows of the active core and in the gaps between the inner reflector blocks. This flow pattern is diagrammed in Fig. 17. The hot spot, located during normal operation at the bottom part of the active core, moves slowly to the upper quarter of the core, as the axial temperature profiles in Fig. 18 demonstrate. The maximum core temperature decreases at first, since the hot spot moves to colder regions until about 10 h after the beginning of the transient due to the development of natural convection in the core. After this, the situation remains the same during the rest of the transient: the decay heat is removed via conduction and radiation to the pressure vessel, and also by the convection field in both the channels and the gaps, which is shown in Fig. 19. The peak fuel temperatures (at the mesh scale) are not found to be much higher during PCC than during normal operation. Nevertheless, the more important consequence is the redistribution of temperature during the event (shown axially in Fig. 18) and the large temperature increase in the fuel in the upper part of the reactor. Moreover, the temperatures computed are mesh-averaged temperatures. The actual peak temperature of the fuel at the center of the compact is almost identical to the mesh-averaged temperature after the reactor trip, but this peak temperature is about 100 ◦ C hotter than then mesh-averaged temperature during full-power operation. The lower structures cool while the upper plenum and its neighboring structures heat up during the first 50 h of the transient. The temperature of the helium in the upper plenum rises from 950 to 1000 ◦ C after 50 h. In the riser, the flow develops into radial and azimuthal loops. This situation is not very stable and may be linked to the choice of a 30◦ sector as the computational domain; a full 360◦ model, or at least a 180◦ model of half of the reactor, would be necessary for a precise study of this region. Nevertheless, the power transferred across the riser is correctly assessed. The core barrel reaches a temperature of 750 ◦ C after 80 h, as shown in Fig. 20, with
Fig. 17. The helium flow pattern during PCC after 10 h of simulated time.
this maximum temperature located at the very upper part of the barrel. Azimuthal convection loops develop in the annuli due to the narrowness of these regions and the three-dimensional dissymmetries. Some of these convection loops are visible in Fig. 21, which presents the velocity field in the annular space of the riser at the end of the 100 h PCC transient. Figs. 22 and 23 present the time evolution of the maximum temperatures of the core and the pressure vessel for both the PCC and DCC scenarios. The DCC is the more severe scenario for the maximum temperature. Nevertheless, the PCC is more penalizing for the upper structures, including the upper part of the pressure vessels (e.g., the flanges where the vessel head attaches). During the PCC, the maximum temperatures occur between 60 h (for the core) and 80 h (for the vessel), which corresponds to the time when the maximum amount of heat is stored
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
1935
Fig. 18. Axial temperature profiles for selected times during PCC.
in the reactor, as shown in Fig. 24. After this period, the power removed from the system is greater than the power produced by the decay heat. Finally, the PCC behavior is, after the first 10 h of the transient, rather similar to a DCC scenario with a different initial temperature distribution inside the reactor (i.e., the temperature pattern is shifted toward the top of the core). Fig. 22
illustrates the two periods of the transient: before 10 h, when the hot spot moves upward due to convection, and after 10 h, which is similar to a DCC transient. The velocity field (Fig. 19) shows the characteristic shape of the convection inside the gaps. This convection acts globally, after it has been established, to increase the radial heat transfer and to vertically shift hot spots toward the top of the vessel. The convective field
Fig. 19. The helium flow in the gaps and the core channels during PCC at 100 h of simulated time.
1936
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
Fig. 20. The maximum temperature of the core barrel during the PCC transient.
Fig. 22. The maximum temperature of the core during the PCC transient (solid line) and during the DCC transient (dashed line).
bounds the energy stored in the early part of the transient by increasing the radial losses, but it limits the peak fuel temperature mainly because it moves the hot spot to previously “cold” areas.
Fig. 23. The maximum temperature of the pressure vessel during the PCC transient (solid line) and during the DCC transient (dashed line).
Fig. 24. The energy stored in the core during the PCC transient.
6. Conclusions
Fig. 21. The flow field in the riser annulus during PCC at 100 h of simulated time.
An important safety characteristic of modern HTRs is passive heat removal in case of a loss of forced convection under either pressurized or depressurized conditions. The depressurized case (DCC) is a classical situation that has been examined by many investigators. The pressurized case (PCC), when the primary loop remains pressurized, has not been investigated as widely, however. A three-dimensional model, using the STAR-CD software, has been built. It solves the combined conductive, convective, and radiative heat transfer using an equivalent homogeneous media for transient events with a duration of 100 h or more.
J.-P. Simoneau et al. / Nuclear Engineering and Design 237 (2007) 1923–1937
Special techniques were employed to simulate the internal natural circulation flows during PCC without slowing the simulation excessively. These are scoping calculations, which do not include the variability of properties or design margins. The objectives of these calculations were both to set up a methodology and to explore the general behavior of HTRs in such situations. Compared to the depressurized situation, the extra convective heat transfer reduces the peak temperature of both the fuel and the pressure vessel during the PCC event. However, it leads to thermal stratification and hence to much higher temperatures in the upper reactor structures. The next planned extension of this model will be to couple it with a detailed model of convection inside the reactor cavity, in order to include the effects of thermal stratification in the surrounding air. The second planned extension will be to explicitly include the remainder of the primary loop outside the reactor vessel. Concerning comparisons with experimental data, separate effects data are available for individual phenomena, but integral test data of the entire facility are not currently available. The need for integral testing is still being evaluated in the industry (given that the phenomena are simpler and fewer in number than for LWR LOCA analysis). Additional separate effect test-
1937
ing could be of value. The greatest benefit would likely come from additional testing of material properties. References Computational Dynamics Limited, 2004. STAR-CD Version 3.24 User Guide. Computational Dynamics Limited, London. Haque, H., Feltes, W., Brinkmann, G., 2004. Thermal response of a high temperature reactor during passive cooling under pressurized and depressurized conditions. In: Proceedings of the Second International Topical Meeting on High Temperature Reactor Technology, Paper No. F02, Beijing, China, September 22–24. Issa, R.I., Ahmadi Befui, B., Beshay, K., Gosman, A.D., 1991. Solution of the implicit discretised reacting flow equations by operator splitting. J. Comp. Phys. 93, 388–410. Mays, B.E., Woaye-Hune, A., Simoneau, J.-P., Gabeloteau, T., Lefort, F., Haque, H., Lommers, L., 2004. The effect of operating temperature on depressurized conduction cooldown for a high temperature reactor. In: Proceedings of the ICAPP’04, Paper No. 4202, Pittsburgh, PA, USA, June 13–17. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation, New York. Petit, D., 2005. Overall simulation of a HTGR plant with the gas adapted MANTA code. In: Proceedings of the NURETH-11 Conference, Avignon, France, October 2–6. Simoneau, J.-P., Champigny, J., 2004. Large eddy simulation of mixing in the outlet plenum of a high temperature reactor: a benchmark exercise. In: Proceedings of the ICONE 12 Conference, Washington, DC, USA, April 25–29.