- Email: [email protected]
CFD study on inlet flow blockage accidents in rectangular fuel assembly
Nuclear Engineering and Design 292 (2015) 177–186
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
CFD study on inlet flow blockage accidents in rectangular fuel assembly Wenyuan Fan, Changhong Peng, Yun Guo ∗ School of Nuclear Science and Technology, University of Science and Technology of China, Huangshan Road, Hefei 230026, Anhui, China
h i g h l i g h t s • • • •
3D CFD and Relap5 simulations on inlet flow blockage are performed. Transient effects are investigated by dynamic mesh technique. Similar flow and power redistributions are predicted in both methods. Local effects of the blockage are captured by CFD method and analyzed.
a r t i c l e
i n f o
Article history: Received 6 March 2015 Received in revised form 17 June 2015 Accepted 19 June 2015
a b s t r a c t Three-dimensional transient CFD simulation of 90% inlet flow blockage accidents in rectangular fuel assembly is performed, using the dynamic mesh technique. One-dimensional steady calculation is done for comparison, using Relap5 code. Similar mass flow rate redistributions and asymmetric power redistributions of the plate in the blocked scenario are obtained. No boiling is predicted in both simulations, however, CFD approach provides more in-depth investigations of flow transients and the thermalhydraulic interaction. The development of flow blockage transients is so fast that the rapid redistribution of mass flow rates occurs in only 0.015 s after the formation of the blockage. As a sequence of the inlet flow blockage, jet-flows and reversed flows occur in the blocked channel. This leads to complex temperature distributions of coolants and fuel plates, in which, the highest coolant temperature no longer occurs around the channel outlet. The present study shows the advantage and significance of the application of three-dimensional transient CFD technique in investigating flow blockage accidents. © 2015 Elsevier B.V. All rights reserved.
1. Introduction The research reactor safety (NRS) analysis is usually conducted on two basic categories of accidents, namely, reactivity insertion accident (RIA) and loss of flow accident (LOFA). A benchmark problem in an idealized generic 10 MW MTR light-water pool-type reactor was defined by IAEA (1980). Though a lot of researches are about RIA (El-Messiry, 2000; Bousbia-salah and Hamidouche, 2005; Varvayanni et al., 2005; Hamidouche and Bousbia-salah, 2010), many investigations on LOFA have been done based on the benchmark problem (Hainouna et al., 2010; El-Morshedy, 2011; Chatzidakisa et al., 2013), most of which use one-dimensional codes. Salama and El-Morshedy (2011) and Salama (2012) showed the importance and advantages of using three-dimensional
∗ Corresponding author. Tel.: +86 13965106983. E-mail addresses: [email protected] (W. Fan), [email protected] (C. Peng), [email protected] (Y. Guo). http://dx.doi.org/10.1016/j.nucengdes.2015.06.016 0029-5493/© 2015 Elsevier B.V. All rights reserved.
computational fluid dynamics (CFD) method to simulate LOFA. In these investigations, mass flow rate of the whole core drops dramatically, which will trigger the protect system to shut down the reactor. However, flow blockage of a coolant channel only leads to LOFA in the blocked channel, and it affects little on the total mass flow rate and outlet temperature of the core. Therefore, it is difficult for the protect system to detect flow blockages. The reasons for flow blockage can be categorized mainly into two groups. One is the buckling of the fuel plate, and the other is obstruction caused by debris that are carried by the coolant. The former is usually caused by irradiation damage and thermal stress, so the formation of this kind of blockage always lasts a relatively long time, while the latter may occur in a very short time due to the fast circulation of the coolant. In the second situation, temperature in the blocked channel will increase rapidly due to the fast drop in mass flow rate, and the integrity of the fuel is, therefore, threatened. In 1963, a plate-type fuel assembly of the Oak Ridge Research Reactor was blocked by the rubber gasket, and ended up in partial melting of three fuel plates (Sim and Tabor, 1964).
178
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
There are two ways to simulate the flow blockage accident. System analysis codes (Adorni et al., 2005; Lu et al., 2009) assume that all the calculated quantities can only change along the flow direction, and use lots of empirical formulas to solve problems. Though motor valve components in these codes can be used to simulate the formation of the blockage, which makes these simulations more practical, they fail to investigate local effects of the obstruction. CFD analysis, which can fully capture the required information of the investigated system, is more applicable to this simulation and becomes the alternative. CFD is increasingly being used in NRS analyses because a number of important phenomena cannot be predicted by traditional one-dimensional system codes with required accuracy and spatial resolution (Bestion, 2004; Höhne et al., 2010). However, computing facilities limit the extension of CFD analysis to the complex and large scale (Yadigaroglu et al., 2003). Fortunately, on the basis of proper selection and simplification of the investigated domain, one can make the simulation reasonable and amenable. For instance, Salama and El-Morshedy (2012) performed a two-dimensional steady CFD investigation on the flow blockage of the channel in the IAEA 10MW MTR research reactor due to the buckling of the fuel plate during normal operation. And the result shows that, when the blockage ratio of the hot channel is larger than 85%, boiling of the coolant may occur. Salama (2012) went a step further by investigating behaviors in the blocked channel when the fast loss of flow accident occurs. And the result indicates that, boiling may occur when the blockage ratio in the average channel is larger than 80%. Recently, Salama et al. (2015) performed a three-dimensional CFD flow blockage investigation, and no boiling is predicted for the 90% blockage. Few CFD investigations have been done toward flow blockages which are caused by debris, because transient simulations are necessary for the formation of the blockage. In addition, dynamic mesh technique is needed to simulate motions of blockages. This investigation process is quite complicated and time-consuming. The present study aims at developing transient modeling method for debris blockages. In the present paper, using dynamic mesh technique, three-dimensional transient CFD simulation on the partial (90%) blockage of the hot channel in the IAEA 10MW research reactor is performed and compared with Relap5 results.
W
G
G
G
G
W
W
SFA 1
SFA 2
SFA 3
SFA 4
W
SFA 5
CFA 22
SFA 6
SFA 7
CFA 23
SFA 8
SFA 9
SFA 10
SFA 12
SFA 13
SFA 14
CFA 24
SFA 15
SFA 16
CFA 25
SFA 17
W
SFA 18
SFA 19
SFA 20
SFA 21
W
W
G
G
G
G
W
SFA 11
H2O
+ AL
SFA 11
Fig. 1. Core configuration of the IAEA 10 MW MTR research reactor (IAEA, 1992).
Fig. 2. Cross-sectional view of a SFA (dimensions in mm) (IAEA, 1992).
2. Core configuration The IAEA MTR research reactor is a benchmark system mainly used to validate institution computer programs. The configuration of the core of IAEA 10 MW MTR research reactor (IAEA, 1992) is depicted in Fig. 1. There are 21 Standard Fuel Assemblies (SFA) and 4 Control Fuel Assemblies (CFA) in the core. The core is reflected by graphite (G) on two opposite faces and surrounded by water (W). Fig. 2 shows the cross-sectional view of a SFA. The related key features of the core are shown in Table 1.
3. Modeling methodology 3.1. Assumptions and simplifications When flow blockage accident occurs, channels in several assemblies may be influenced. However, it is not wise to investigate the accident in such a big scale, even in the scale of one full fuel assembly. According to the work of Lu et al. (2009) and Salama and El-Morshedy (2012), parameters in the blocked channel and its adjacent channel are affected mostly, while those in the third channel are slightly influenced. Since the geometry of the channels and plates is symmetric, only a quarter of the domain is considered to reduce the mesh size. Therefore, the simulation domain is set to
Table 1 Specifications of the IAEA 10 MW MTR research reactor (IAEA, 1992). Nuclear fuel Fuel element Coolant Coolant flow direction Moderator Fuel thermal conductivity(W/cm K) Cladding thermal conductivity (W/cm K) Fuel specific heat (J/g K) Cladding specific heat (J/g K) Fuel density (g/cm3 ) Cladding density (g/cm3 ) Radial peaking factor Axial peaking factor Engineering peaking factor Inlet coolant temperature (◦ C) Operating pressure (bar) Length (mm) Width (mm) Height (mm) Number of fuel elements in SFA/CFA Plate meat (mm) Width (mm) active/total Height (mm) Water channel thickness (mm) Plate clad thickness (mm)
MTR Plate-type clad in Al Light water Downward Graphite-light water 1.58 1.8 0.728 0.892 0.68 2.7 1.4 1.5 1.2 38.0 1.7 80.5 76 600 23/17 0.51 63/66.5 600 2.23 0.38
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
three channels and corresponding fuel plates, which is the red line surrounded area in Fig. 2. In the present study, no buckling of the fuel plate is considered during the calculation, and the blockage is caused by debris only. Therefore, the flow blockage can only occur at the inlet of the channel. However, debris carried by the coolant may be of different kinds of forms and shapes. In this paper, it is assumed that the blockage is caused by a zero-thickness plane, and the area of the plane is 90% of the blocked channel inlet size. However, it should be emphasized that blockages of any given shapes can be modeled and investigated by means of CFD methods. The investigated domain is depicted in Fig. 3, in which, an upper plenum and a lower plenum are added to give reasonable inlet and outlet conditions conveniently. During the initial calculation of the simulation, the blockage plane is stationary at the upper plenum. When transient calculation begins, the blockage plane starts to move with a speed of 2.3 m/s, which is the main stream velocity in the upper plenum. Finally, the blockage plane arrives at the inlet of the blocked channel and stops there. During this process, the heat generation rate of the fuel is set to a chopped cosine function with the effective length set to 0.6 m, and the extrapolated length is 0.76 m (Salama and El-Morshedy, 2012). Due to the potential big rise in temperature, the density, specific heat, thermal conductivity and viscosity of the coolant are all considered to change with temperature at the form of quartic polynomials.
accurate near-wall treatment is the wall functions approach. Wall functions are a series of semi-empirical formulas, linking the viscosity-affected region and the main stream without solving the viscosity-affected inner region. In these functions, there is a significant non-dimensional parameter y* defined as follows: 1/4 1/2
y∗ =
C kP yP
Direct Numerical Simulation (DNS) of turbulent flows is too computationally expensive to be applied to complex geometries. Reynolds-averaging and filtering are the alternative methods to simulate turbulent flow, but the latter also requires a significant amount of computational resources. The Reynolds-averaged Navier–Stokes (RANS) equations are time-averaged equations of motion for fluid flow, which are computationally less expensive to solve and the results can be acceptable for engineering. The RANS technique therefore is widely adopted for practical engineering applications. But RANS equations introduce additional unknown variables (Reynolds Stress Tensor) and turbulence models are established in order to achieve the problem closure and determine these variables in terms of known quantities. The realizable k–ε mode is applied in this study, and the modeled transport equations for k and ε in the realizable k–ε model (ANSYS, Inc., 2009) are
3.3. Meshing strategies Structured mesh is used to discretize the simulated domain which is of high aspect ratio. The combination of wall functions and structured mesh reduces the mesh quantity significantly, making it reasonably accurate and economical for the transient calculation. To simulate the moving plane, dynamic mesh technique must be applied. The layering type dynamic mesh in ANSYS FLUENT is designed for structured mesh, allowing motion along a specific straight line, which makes it suitable for the current simulation. Moving zones and stationary zones in the layering method are connected by sliding interfaces, through which data from one zone are interpolated and transferred to the other. Errors occur during
⎧
t ∂k ∂ ∂ ∂ ⎪ ⎪ ku = + + Gk + Gb − εA − YM − Sk + (k) ⎪ j ⎨ ∂t k ∂xj ∂xj ∂xj
⎪ ∂ ∂ ∂ ε2 ε t ∂ε ⎪ ⎪ εu = + + C1 Sε − C2 + √ + C1ε C3ε Gb + Sε ⎩ (ε) j ∂xj
where
C1 max 0.43,
∂xj
+5
,
k =S , ε
ε
S=
∂xj
k+
2Sij Sij
In these equations, Gk represents the generation of turbulence kinetic energy due to the mean velocity gradients. Gb is the generation of turbulence kinetic energy due to buoyancy. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. C1ε , C2ε , and C3ε are constants. k and ε are the turbulent Prandtl numbers for k and ε, respectively. Sk and Sε are source terms. On the other hand, like other k–ε models, the realizable k–ε model is primarily valid for turbulent flows in the main stream. Special attention, therefore, needs to be paid to simulate turbulent flows with confining walls through the use of an appropriate near-wall treatment. An economical, robust, and reasonably
(2)
where P is the near-wall node; C is a constant for standard and RNG k–ε models, but is a variable in realizable k–ε model; kP is turbulence kinetic energy at P; yP is distance from point P to the wall; is kinematic viscosity of the fluid. Among these functions, standard wall functions are popular for its good performances in a wide range of wall-bounded flows. But in unsteady flows, in which the mass flow rate changes with time, y* will change. Especially when the mass flow rate decays dramatically, y* may become smaller than 11.225, which is the intersection of the logarithmic and the linear near-wall profile, and this will lead to deterioration of standard wall functions. When flow blockage occurs, the mass flow rate of the channel will decrease significantly, resulting from the increase in flow resistance. Standard wall functions therefore cannot be used to simulate the channel flow blockage transient. On the basic assumption that the viscous sub-layer is so thin that it has the same boundary with the wall, scalable wall functions (ANSYS, Inc., 2012) overcome the drawback by introducing a limiter in the y* calculations such that y˜ ∗ = max(y∗ , 11.225), and y* used for any formula in standard wall functions is replaced by y˜ ∗ . This feature of scalable wall functions makes them suitable for the present study.
3.2. CFD techniques and turbulence models
∂t
179
ε
(1)
k
interpolations, which also require extra computing time. Therefore, the size of moving zone should be minimized to lessen interface size, which can reduce numerical errors and accelerate the simulation. In this study, to guarantee the accuracy of heat transfer calculation between claddings and coolants, all interfaces are set between fluid zones and fluid zones only. Obviously, the zone between the initial plane location and the inlet of the blocked channel is being compressed during the calculation. However, this zone becomes zero-volume when the plane finally arrives at the inlet, leading to change in topology structure of domains, which is not allowed by the solver. To avoid the alteration in topology structure, an L-grid division is conducted in the blocked channel block, as illustrated in Fig. 4. The short branch near the inlet and the compressed zone in the upper plenum are combined to a new compressed zone, and the long branch beside the fuel plate is set to boundary layer mesh.
180
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
Fig. 3. Modeling of the investigated domain. (a) Front view. (b) Oblique view.
According to Best Practice Guidelines for the use of CFD in NRA (Mahaffy and Farm, 2010), grid sensitivity analysis should be conducted to make sure that the simulation is independent of mesh size. Velocity magnitude and relative static pressure of the centroid of channel-1 are compared, as shown in Table 3. Several grids have been examined; when the node number is increased to 1797,576 and 2205,604, both velocity magnitude and relative static pressure differ by less than 0.8%. When a more dense mesh with 3861,684 nodes is taken into consideration, the differences in the two finer meshes are less than 0.05%. Therefore, the 2205,604 nodes mesh is chosen for further calculations, a cross-sectional view of the mesh is depicted in Fig. 5. As advised by ANSYS FLUENT, the time step size of transient calculations is chosen from 5 × 10−6 s to 4 × 10−4 s during different simulation periods to satisfy the convergence requirements. Fig. 4. The L-grid division strategy.
Therefore, even though the plane arrives at the inlet, the new compressed zone still has a size of the short branch, which is accepted by the solver. 4. Results and discussions Because of its good performance in unsteady simulations, the PISO method is employed for pressure–velocity coupling. The gradient term and the pressure term are calculated by interpolation using the least squares cell based scheme and the PRESTO! scheme, respectively. The second order upwind scheme is used for all the variables. The first order implicit scheme is applied for the discretization of the transient formulation. According to Salama and El-Morshedy (2012), the average inlet velocity of the upper plenum is set to 2.3 m/s. In order to define an accurate flow boundary for the inlet, fully developed velocity and turbulence profiles are added to the inlet, and these quantities are fixed during transient calculations. For the lower plenum outlet, the boundary type is set to pressure-outlet. To summarize, main boundary and initial conditions are listed in Table 2.
4.1. Validation consideration According to Best Practice Guidelines for the use of CFD in NRA (Mahaffy and Farm, 2010), model validation should be performed before application to real case studies. However, due to the compact structure and extremely high heat flux, it is difficult to conduct accurate experiment on flow blockage accidents. Besides, transient experiment data are needed to validate results of the current study. Therefore, it is hard to validate results of the present study. In order to compare with the CFD results, system analysis code Relap5(Mod3.4) is used to calculate the steady results in both
Table 2 Boundary and initial conditions. Inlet of the upper plenum Outlet of the lower plenum Heat generation rate of fuels Moving speed of the blockage plane Duration of the plane moving transient
Fully developed velocity-inlet, average velocity is 2.3 m/s Pressure-outlet Chopped cosine function 2.3 m/s 0.085 s
Fig. 5. Cross-sectional view of the mesh.
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
181
Table 3 Grid sensitivity analysis. Mesh
Cells
Nodes
Velocity
1 2 3 4
1019,208 1470,198 1868,337 3246,042
1316,064 1797,576 2205,604 3861,684
4.126431 4.137421 4.140962 4.139675
non-blocked and blocked states of the domain, the corresponding model is shown in Fig. 6. System codes are widely used for their simplicity and relatively low requirements on computing resources, however there are some disadvantages in the use of such system analysis codes. As mentioned before, empirical formulas are used in one-dimensional codes. For instance, during the calculating process, Tw (z), the coolant temperature, is calculated according to the law of energy conservation first, and then Tcs (z), the cladding surface temperature, is calculated as follows: Tcs (z) = Tw (z) +
q (z) h (z)
(3)
where q (z) is the heat flux of the cladding surface, and h(z) is the convection heat transfer coefficient, which is derived from empirical formulas. However, there are many different empirical formulas that can be used for h(z) calculation, therefore different Tcs (z) are predicted. For instance, in Salama (2011), the Dittus–Bolter equation (Winterton, 1998) provides consistent results with the CFD simulation on one channel, while in Salama et al. (2015), the Sieder and Tate (1936) equation gives consistent predictions with the CFD calculation on more channels. Therefore, it is the choice of empirical equations that determines the Tcs (z) calculation. On the other hand, different modeling strategies toward the assembly structure may also influence the result. As depicted in Fig. 2, along the width direction, there are active zones and inactive zones in fuel plates. In CFD approaches, models are precisely built according to the assembly structure, while during Relap5 modeling, due to the limitation of the code, fuels and claddings are assumed
Fig. 6. The Relap5 model.
Error (%) 0.266 0.086 −0.031
Relative static pressure 10,241.457 10,298.699 10,377.161 10,382.158
Error (%) 0.559 0.762 0.048
to have the same width. This simplification changes the plate size, and the temperature field of the plate is changed consequently. Due to the differences in modeling and solving methods, like other system codes, Relap5 results are only used to compare with CFD predictions. For the non-blocked scenario, the coolant, cladding surface and fuel temperature profiles along the channel-1 in both simulations are depicted in Fig. 7. For the CFD approach, temperature is calculated for each volume. In order to compare with the onedimensional results, temperatures should be averaged in every cross section along the channel. However, proper weight should be selected when calculating the average temperature of the coolant. The area-weighted temperature and the z-velocity-weighted temperature represent the average temperature of the cross section statically and dynamically, respectively. As shown in Fig. 7, for the non-blocked scenario, these two kinds of average values are quite close. When the comparison is done between results calculated by CFD and Relap5, good closeness is also obtained for coolant temperatures. Though the cladding surface and fuel temperatures calculated by CFD method are higher than those obtained by Relap5, differences in the models can explain differences in results. 4.2. Velocity transients Like system analysis, mass flow rate of a channel can be obtained by CFD method. Besides, there many more details of the simulated domain can be obtained through CFD method. 4.2.1. Mass flow rate transients It is obvious that when the plane stops at the inlet of the blocked channel, flow resistance increases significantly, resulting in fast changes in mass flow rate of the investigated channels, especially in the blocked one. In this paper, for the sake of comparisons, all mass flow rates are converted to those in the half-channel correspondingly. Transient mass flow rates in each channel are shown in Fig. 8. When the plane is moving with the coolant, mass flow rate varies little. While when the plane arrives at the inlet, mass flow rate in
Fig. 7. Temperature profiles along channel-1(non-blocked scenario).
182
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
Fig. 8. Transitional mass flow rate in different channels.
the blocked channel drops dramatically, with those in the other two channels both increasing. The plane stops at 0.085 s after the simulation, when flow time comes to 0.1 s, the rapid change in mass flow rate of all the channels converts to a quite slow one. Therefore, the redistribution of the mass flow rate can be considered to complete in only 0.015 s. However, the result of the redistribution is different from that acquired by Salama and El-Morshedy (2012), which suggests that flow in channel-3 almost remains unchanged in different blockage ratio scenarios. The difference results from the different blockage modeling methods. In their investigation, the blockage is caused by buckling of the plate. The sizes of all the inlets and outlets are consistent with those in the non-blocked scenario, so the mass flow rate redistribution is mainly driven by flow resistances along different channels. Considering there is a flow area expansion in channel-2 corresponding to the flow area contraction in channel-1, and the flow area in channel-3 is unchanged, therefore, the flow resistance of channel-3 is larger than that in channel-2, but lower than that in channel-1. Being the channel with the smallest flow resistance, channel-2 has the highest mass flow rate in blocked scenarios. In the present study, the blockage occurs at the inlet of the flow channel, and it is the local flow resistances of different inlets that determine the mass flow rate redistribution. Compared to channel3, channel-2 is closer to the blockage, and is more affected by the blockage. Therefore the mass flow rate of channel-2 is lower than that of channel-3. The dissimilarity in mass flow rate redistribution prediction shows that the blockage form affects the flow field in the assembly strongly. For blockages caused by buckling, more coolant flows through channel-2 than that in inlet flow blockage cases, so the fuel plate adjacent to the channel-1 is cooled better. Therefore, the inlet flow blockage may lead to a severer consequence than the buckling blockage does. 4.2.2. Reversed flow around the outlet of the blocked channel A small reversed flow region occurs around the outlet of the blocked channel at about 0.1 s. Fig. 9 shows the z-velocity distribution at the outlet of the blocked channel. This tiny region then becomes larger, after that reaches a relatively steady size. However, it is hard to explain the reversed flow around the outlet, especially for channels through which coolant flows downward. As will be seen later, due to the presence of the blockage, quite a large area of coolant in the blocked channel are influenced, and flow field is not fully developed in channel-1. In addition, when coolant from channels flows to the lower plenum, the flow field is quite
Fig. 9. Profiles of z-velocity on the outlet of the blocked channel.
complex, even in normal scenario. For the blocked condition, mass flow rate in the blocked channel is much lower than other channels, which makes the flow field more complicated. Therefore, the large mass flow rate difference between the blocked and unblocked channels should be another possible reason. 4.2.3. Flow field in the blocked channel The flow field in channel-1 changes rapidly during the calculation, especially when the plane is about to arrive at the inlet and a very short time after its arrival. As shown in Fig. 10, because of the shrink of the flow area and the inertia of the upstream coolant, a high speed jet-flow occurs when the plane stops at the inlet, which leads to backflow behind the plane at the same time. The maximum velocity in the blocked channel rockets to as high as 50 m/s, and the corresponding maximum backflow velocity also reaches its peak, which is almost 13 m/s, and maximum and minimum values of them during the simulation are depicted in Fig. 11. Again, the fast change takes place in about 0.015 s. A problem therefore arouses, that the high speed jet-flow may threaten the assembly structure. As depicted in Fig. 12, the high speed region spreads instantly from the inlet to the downstream of the blocked channel after the formation of the blockage. Due to the viscosity of the coolant, a larger region is affected by the blockage during the spread, with the maximum velocity decaying sharply both in flow direction and backflow direction. Under the influence of both the high speed region and the small reversed flow region around the outlet, a new backflow region occurs around the middle part of the channel, in about 1.5 s after the calculation. After that, the new region grows larger, and the velocity field structure becomes stable. 4.3. Temperature transients As a consequence of the redistribution of the flow field, the temperature field of the domain changes with time. The maximum temperatures of the cladding and the fuel during calculation are depicted in Fig. 13. Since the melting point of the cladding is 855.15 K, the peak temperatures during calculation are far from their corresponding melting points, therefore the integrity of the fuel plate is maintained in this simulation. The boiling point
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
Fig. 10. Profiles of flow direction velocity in the blocked channel (flow time form 0.085 s to 0.0999 s).
of the coolant is 388.32 K under the operating pressure, therefore no boiling of the coolant occurs during the transient. The peak temperature around 0.8 s can be explained as follow. As the power distribution is a chopped cosine function along the channel, the peak always appears around the middle of the plate. The mass flow rate of the blocked channel decreases sharply after the formation of the blockage. Temperatures, except those in high speed regions, rise dramatically due to the bad heat transfer condition. As mentioned before, the high speed regions will spread to the
Fig. 11. The maximum and minimum z-velocity in the blocked channel.
183
Fig. 12. Profiles of flow direction velocity in the blocked channel (flow time form 0.299 s to 1.998 s).
down stream; at about 0.8 s, the middle part is involved, as shown in Fig. 12. High speed coolant with lower temperature enhances heat transfer in this region. As a consequence, the maximum temperature starts to drop at about 0.8 s. The previous discussion also indicates that, the high speed regions will slow down and then a backflow region will occur around the middle part of the channel. This leads to the rebound of the highest temperatures of the cladding and fuel. This explains the closeness between the time of the formation of the backflow region and that of the start of the bouncing.
Fig. 13. Maximum temperatures during calculation.
184
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
Fig. 14. Temperature and velocity field in the blocked channel.
4.4. Steady results of blocked scenario 4.4.1. Details calculated by CFD approach As the calculation continues, both flow field and temperature field of the investigated domain gradually transform to a new steady distribution. Fig. 14 shows the temperature and velocity fields in the blocked channel. A big vortex occurs behind the blockage, extending to the middle part. There is another vortex below the big one, but with small scale. As mentioned before, there is a reversed flow region around the outlet. Actually, there are two tiny vortex regions below the outlet. The highest coolant temperature does not occur in the middle of the channel where is the hottest area of the plate, as shown in Fig. 15, because coolant around the middle part is of high speed. It is the region where the largest vortex occurs that the maximum coolant temperature appears, because of the relatively slow circle motion of the coolant. For the cladding, besides the highest temperature in the middle, another small area is of a relatively high temperature a little past the middle, because of the corresponding vortex. 4.4.2. Comparisons with Relap5 To simulate the blocked scenario in Relap5, the flow area of the single junction, which is used to connect the upper plenum and the blocked channel, is reduced to 10% of its initial value. The mass flow rates in all channels are depicted in Fig. 16, in which mass flow rate in channel-1 calculated by FLUENT is 19.8% higher than that simulated by Relap5. This large difference results from the quite small mass flow rates in channel-1; actually, the difference between two predictions is about 0.00812 kg/s, which is only 1.2% of the total mass flow rate. Close predictions are obtained for mass flow rates in channel-2 and channel-3, and the dissimilarity is 1.8% and 0.5%, respectively. However, different profiles for channel-2 and channel-3 in two treatments appear.
Fig. 15. Cladding temperature field of the blocked channel.
In CFD approach, as the adjacent channel of the blocked channel, channel-2 is more affected by the blockage than channel-3. More coolants, therefore, pass through channel-3 than channel-2. While in Relap5, channel-2 and channel-3 are not recognizable in geometric positions. Therefore, compared to channel-3, a larger mass flow
Fig. 16. Mass flow rate comparisons.
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
185
Fig. 17. Power redistribution comparisons. Fig. 19. Temperature profiles along channel-1 (blocked scenario).
rate for channel-2 is predicted, as a sequence of the higher heating power in channel-2, which will be discussed in the following. Because of the big mass flow rate differences of the blocked channel and its adjacent one, power generated in the fuel plate between them redistributes, as shown in Fig. 17, and quite similar profiles are obtained. Compared to coolants in channel-2, those in channel-1 are of very low Re number. Therefore, the heat transfer coefficients on both sides of fuel Plate 1are quite different. For the wall side toward channel-1, the heat transfer coefficient is much lower than the coefficient on the wall side toward channel-2. Consequently, huge asymmetric heat transfer occurs in fuel Plate 1. While for fuel Plate 2, due to the little difference in mass flow rates of channel-2 and channel-3, the distributions are almost symmetric for both simulations. Fig. 18 shows temperature profiles of different outlet in both normal and blocked cases. The outlet temperature of the blocked channel dramatically rises due to the huge drop in mass flow rate of the channel. As mentioned above, channel-1 are almost equally heated in both calculations, however a higher mass flow rate of channe-1 is predicted by CFD approach. The outlet temperature, therefore, is a little lower than that calculated by Relap5 codes. Due to the total mass flow rate is fixed in both simulations, the mass flow rates in channel-2 and channel-3 both become larger. Channel-2 gets more energy form fuel Plate 1, therefore, the outlet temperature rises a little. For channel-3, the energy source is almost
same with the normal case, increase in mass flow rate leads to a slight drop in outlet temperature. Both the comparisons of heat transfer ratios of fuel plates and outlet temperature profiles indicate that, when considering heat transfer effects of the flow blockage, only the blocked channel and plates that constitute the blocked channel are mostly affected. In another word, the blocked channel and its adjacent ones are significant for investigation. The coolant, cladding surface and fuel temperature profiles along the flow direction are compared in Fig. 19. For Relap5 treatment, except for the detailed values, tendency of profiles are quite similar to those in the normal scenario. This indicates that flow details in the blocked channel are totally ignored. For CFD approach, three-dimensional results are converted to one-dimensional profiles, and the shapes of the cladding surface temperature and fuel temperature profiles are different from those in the non-blocked case, especially around the middle of the plate. This difference indicates that power redistribution occurs in the flow direction, which is not considered in Relap5 treatment. Since a relatively higher temperature is predicted in the CFD approach, the z-velocity-weighted temperature is lower than that calculated by Relap5 with a similar profile shape. However, the area-weighted temperature is quite irregular due to the complex flow and temperature fields. In Relap5, outlet temperature is the highest temperature of the coolant, as shown in Fig. 18, and it is higher than the outlet temperature that is predicted by CFD method. However, the value is still lower than the maximum area-weighted temperature in CFD approach. Therefore, due to the neglect of flow details, conservative perdition in the mass flow rate cannot guarantee a conservative temperature result. 5. Conclusions and outlooks
Fig. 18. Temperature profiles of different outlet in both normal and blocked cases.
Using dynamic mesh technique, the 90% inlet flow blockage of an MTR plate-type assembly is investigated by three dimensional CFD methods. Relap5 simulation is conducted for comparison. No boiling of the coolant is predicted, and both CFD and Relap5 simulations indicate that, when concerning heat transfer effects, the blocked channel and its adjacent ones are mostly affected. Besides, quite similar asymmetry power redistributions are predicted by these two methods. The transient CFD simulation shows that the development of flow blockage transients is so fast that the rapid redistribution of mass flow rates occurs in only 0.015 s after the formation of the blockage. Besides, more significant details of the domain are captured by the CFD method, such as the jet-flow and vortexes in the blocked channel, and the irregular temperature profiles caused by the blockage.
186
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
The present study shows advantages of using transient CFD method to investigate flow blockages which are caused by debris. However, relevant experiments are needed for validations. As mentioned before, it is quite difficult to conduct experiment. Using current results, targeted measurements may make the experiment more practical. Besides, further simulations can be conducted on larger domain and toward different forms of debris. Acknowledgements The authors would like to thank the Science and Technology on Reactor System Design Technology Laboratory of Nuclear Power Institute of China for financial and experiential support of this work. This study is also supported by National Natural Science Foundation of China (Grant no. 11305169). References Adorni, M., Bousbia-Salah, A., Hamidouche, T., et al., 2005. Analysis of partial and total blockage of a single fuel assembly of an MTR research reactor core. Ann. Nucl. Energy 32, 1679–1692. ANSYS, Inc., 2009. ANASYS FLUENT 12.0 Theory Guide. ANSYS, Inc., Canonsburg, PA, USA. ANSYS, Inc., 2012. ANASYS 14.5 Help. ANSYS, Inc., Canonsburg, PA, USA. Bestion, D., 2004. Recommendation on use of CFD codes for nuclear reactor safety analysis. In: EVOL-ECORA-D14. Bousbia-salah, A., Hamidouche, T., 2005. Analysis of the IAEA research reactor benchmark problem by the RETRAC-PC code. Nucl. Eng. Des. 235, 661–674. Chatzidakisa, S., Ikonomopoulosa, A., Ridikasb, D., 2013. Evaluation of RELAP5/MOD3 behavior against loss of flow experimental results from two research reactor facilities. Nucl. Eng. Des. 255, 321–329. El-Messiry, A.M., 2000. Cobalt irradiation box ejection accident of ETRR-2. Nucl. Eng. Des. 198, 287–293. El-Morshedy, S.E., 2011. Prediction, analysis and solution of the flow inversion phenomenon in a typical MTR-reactor with upward core cooling. Nucl. Eng. Des. 241, 226–235.
Hainouna, A., Ghazia, N., Mansour Abdul-Moaizb, B., 2010. Safety analysis of the IAEA reference research reactor during loss of flow accident using the code MERSAT. Nucl. Eng. Des. 240, 1132–1138. Hamidouche, T., Bousbia-salah, A., 2010. Assessment of RELAP5 point kinetic model against reactivity insertion transient in the IAEA 10MW MTR research reactor. Nucl. Eng. Des. 240, 672–677. Höhne, T., Krepper, E., Rohde, U., 2010. Application of CFD codes in nuclear reactor safety analysis. Sci. Technol. Nucl. Instal. 2010, 1–8. IAEA, 1980. IAEA Research Reactor Core Conversion from the use of high-enriched uranium to the use of low enriched uranium fuels Guidebook. In: IAEA-TECDOC233. IAEA. IAEA, 1992. Research reactor core conversion guidebook, volume 3: analytical verification (Appendices G and H). In: IAEA-TECDOC-643. IAEA. Lu, Q., Qiu, S., Su, G.H., 2009. Flow blockage analysis of a channel in a typical material test reactor core. Nucl. Eng. Des. 239, 45–50. Mahaffy, J., Farm, W.S., 2010. Development of best practice guidelines for CFD in nuclear reactor safety. Nucl. Eng. Des. 42, 377–381. Salama, A., 2011. CFD investigation of flow inversion in typical MTR research reactor undergoing thermal–hydraulic transients. Ann. Nucl. Energy 38, 1578–1592. Salama, A., 2012. CFD analysis of fast loss of flow accident in typical MTR reactor under going partial and full blockage: the average channel scenario. Prog. Nucl. Energy 60, 1–13. Salama, A., El-Amin, M.F., Sun, S., 2015. Three-dimensional, numerical investigation of flow and heat transfer in rectangular channels subject to partial blockage. Heat Transfer Eng. 36, 152–165. Salama, A., El-Morshedy, S.E., 2011. CFD simulation of the IAEA 10 MW generic MTR reactor under loss of flow transient. Ann. Nucl. Energy 38, 564–577. Salama, A., El-Morshedy, S.E., 2012. CFD simulation of flow blockage through a coolant channel of a typical material testing reactor core. Ann. Nucl. Energy 41, 26–39. Sieder, E.N., Tate, G.E., 1936. Heat transfer and pressure drop of liquids in tubes. Ind. Eng. Chem. 28, 1429–1435. Sim, T.M., Tabor, W.H., 1964. Report on Fuel-plate Melting at the Oak Ridge Research Reactor: ORNL-TM-627. ORNL, USA, Oak Ridge. Varvayanni, M., Stakakis, E., Catsaros, N., Antonopoulos-Domis, M., 2005. Void induced reactivity in a mixed MTR core. Nucl. Eng. Des. 235, 855–865. Winterton, R.H.S., 1998. Where did the Dittus and Boelter equation come from? Int. J. Heat Mass Transfer 41, 809–810. Yadigaroglu, G., Andreani, M., Dreier, J., Coddington, P., 2003. Trends and needs in experimentation and numerical simulation for LWR safety. Nucl. Eng. Des. 221, 205–223.
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
CFD study on inlet flow blockage accidents in rectangular fuel assembly Wenyuan Fan, Changhong Peng, Yun Guo ∗ School of Nuclear Science and Technology, University of Science and Technology of China, Huangshan Road, Hefei 230026, Anhui, China
h i g h l i g h t s • • • •
3D CFD and Relap5 simulations on inlet flow blockage are performed. Transient effects are investigated by dynamic mesh technique. Similar flow and power redistributions are predicted in both methods. Local effects of the blockage are captured by CFD method and analyzed.
a r t i c l e
i n f o
Article history: Received 6 March 2015 Received in revised form 17 June 2015 Accepted 19 June 2015
a b s t r a c t Three-dimensional transient CFD simulation of 90% inlet flow blockage accidents in rectangular fuel assembly is performed, using the dynamic mesh technique. One-dimensional steady calculation is done for comparison, using Relap5 code. Similar mass flow rate redistributions and asymmetric power redistributions of the plate in the blocked scenario are obtained. No boiling is predicted in both simulations, however, CFD approach provides more in-depth investigations of flow transients and the thermalhydraulic interaction. The development of flow blockage transients is so fast that the rapid redistribution of mass flow rates occurs in only 0.015 s after the formation of the blockage. As a sequence of the inlet flow blockage, jet-flows and reversed flows occur in the blocked channel. This leads to complex temperature distributions of coolants and fuel plates, in which, the highest coolant temperature no longer occurs around the channel outlet. The present study shows the advantage and significance of the application of three-dimensional transient CFD technique in investigating flow blockage accidents. © 2015 Elsevier B.V. All rights reserved.
1. Introduction The research reactor safety (NRS) analysis is usually conducted on two basic categories of accidents, namely, reactivity insertion accident (RIA) and loss of flow accident (LOFA). A benchmark problem in an idealized generic 10 MW MTR light-water pool-type reactor was defined by IAEA (1980). Though a lot of researches are about RIA (El-Messiry, 2000; Bousbia-salah and Hamidouche, 2005; Varvayanni et al., 2005; Hamidouche and Bousbia-salah, 2010), many investigations on LOFA have been done based on the benchmark problem (Hainouna et al., 2010; El-Morshedy, 2011; Chatzidakisa et al., 2013), most of which use one-dimensional codes. Salama and El-Morshedy (2011) and Salama (2012) showed the importance and advantages of using three-dimensional
∗ Corresponding author. Tel.: +86 13965106983. E-mail addresses: [email protected] (W. Fan), [email protected] (C. Peng), [email protected] (Y. Guo). http://dx.doi.org/10.1016/j.nucengdes.2015.06.016 0029-5493/© 2015 Elsevier B.V. All rights reserved.
computational fluid dynamics (CFD) method to simulate LOFA. In these investigations, mass flow rate of the whole core drops dramatically, which will trigger the protect system to shut down the reactor. However, flow blockage of a coolant channel only leads to LOFA in the blocked channel, and it affects little on the total mass flow rate and outlet temperature of the core. Therefore, it is difficult for the protect system to detect flow blockages. The reasons for flow blockage can be categorized mainly into two groups. One is the buckling of the fuel plate, and the other is obstruction caused by debris that are carried by the coolant. The former is usually caused by irradiation damage and thermal stress, so the formation of this kind of blockage always lasts a relatively long time, while the latter may occur in a very short time due to the fast circulation of the coolant. In the second situation, temperature in the blocked channel will increase rapidly due to the fast drop in mass flow rate, and the integrity of the fuel is, therefore, threatened. In 1963, a plate-type fuel assembly of the Oak Ridge Research Reactor was blocked by the rubber gasket, and ended up in partial melting of three fuel plates (Sim and Tabor, 1964).
178
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
There are two ways to simulate the flow blockage accident. System analysis codes (Adorni et al., 2005; Lu et al., 2009) assume that all the calculated quantities can only change along the flow direction, and use lots of empirical formulas to solve problems. Though motor valve components in these codes can be used to simulate the formation of the blockage, which makes these simulations more practical, they fail to investigate local effects of the obstruction. CFD analysis, which can fully capture the required information of the investigated system, is more applicable to this simulation and becomes the alternative. CFD is increasingly being used in NRS analyses because a number of important phenomena cannot be predicted by traditional one-dimensional system codes with required accuracy and spatial resolution (Bestion, 2004; Höhne et al., 2010). However, computing facilities limit the extension of CFD analysis to the complex and large scale (Yadigaroglu et al., 2003). Fortunately, on the basis of proper selection and simplification of the investigated domain, one can make the simulation reasonable and amenable. For instance, Salama and El-Morshedy (2012) performed a two-dimensional steady CFD investigation on the flow blockage of the channel in the IAEA 10MW MTR research reactor due to the buckling of the fuel plate during normal operation. And the result shows that, when the blockage ratio of the hot channel is larger than 85%, boiling of the coolant may occur. Salama (2012) went a step further by investigating behaviors in the blocked channel when the fast loss of flow accident occurs. And the result indicates that, boiling may occur when the blockage ratio in the average channel is larger than 80%. Recently, Salama et al. (2015) performed a three-dimensional CFD flow blockage investigation, and no boiling is predicted for the 90% blockage. Few CFD investigations have been done toward flow blockages which are caused by debris, because transient simulations are necessary for the formation of the blockage. In addition, dynamic mesh technique is needed to simulate motions of blockages. This investigation process is quite complicated and time-consuming. The present study aims at developing transient modeling method for debris blockages. In the present paper, using dynamic mesh technique, three-dimensional transient CFD simulation on the partial (90%) blockage of the hot channel in the IAEA 10MW research reactor is performed and compared with Relap5 results.
W
G
G
G
G
W
W
SFA 1
SFA 2
SFA 3
SFA 4
W
SFA 5
CFA 22
SFA 6
SFA 7
CFA 23
SFA 8
SFA 9
SFA 10
SFA 12
SFA 13
SFA 14
CFA 24
SFA 15
SFA 16
CFA 25
SFA 17
W
SFA 18
SFA 19
SFA 20
SFA 21
W
W
G
G
G
G
W
SFA 11
H2O
+ AL
SFA 11
Fig. 1. Core configuration of the IAEA 10 MW MTR research reactor (IAEA, 1992).
Fig. 2. Cross-sectional view of a SFA (dimensions in mm) (IAEA, 1992).
2. Core configuration The IAEA MTR research reactor is a benchmark system mainly used to validate institution computer programs. The configuration of the core of IAEA 10 MW MTR research reactor (IAEA, 1992) is depicted in Fig. 1. There are 21 Standard Fuel Assemblies (SFA) and 4 Control Fuel Assemblies (CFA) in the core. The core is reflected by graphite (G) on two opposite faces and surrounded by water (W). Fig. 2 shows the cross-sectional view of a SFA. The related key features of the core are shown in Table 1.
3. Modeling methodology 3.1. Assumptions and simplifications When flow blockage accident occurs, channels in several assemblies may be influenced. However, it is not wise to investigate the accident in such a big scale, even in the scale of one full fuel assembly. According to the work of Lu et al. (2009) and Salama and El-Morshedy (2012), parameters in the blocked channel and its adjacent channel are affected mostly, while those in the third channel are slightly influenced. Since the geometry of the channels and plates is symmetric, only a quarter of the domain is considered to reduce the mesh size. Therefore, the simulation domain is set to
Table 1 Specifications of the IAEA 10 MW MTR research reactor (IAEA, 1992). Nuclear fuel Fuel element Coolant Coolant flow direction Moderator Fuel thermal conductivity(W/cm K) Cladding thermal conductivity (W/cm K) Fuel specific heat (J/g K) Cladding specific heat (J/g K) Fuel density (g/cm3 ) Cladding density (g/cm3 ) Radial peaking factor Axial peaking factor Engineering peaking factor Inlet coolant temperature (◦ C) Operating pressure (bar) Length (mm) Width (mm) Height (mm) Number of fuel elements in SFA/CFA Plate meat (mm) Width (mm) active/total Height (mm) Water channel thickness (mm) Plate clad thickness (mm)
MTR Plate-type clad in Al Light water Downward Graphite-light water 1.58 1.8 0.728 0.892 0.68 2.7 1.4 1.5 1.2 38.0 1.7 80.5 76 600 23/17 0.51 63/66.5 600 2.23 0.38
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
three channels and corresponding fuel plates, which is the red line surrounded area in Fig. 2. In the present study, no buckling of the fuel plate is considered during the calculation, and the blockage is caused by debris only. Therefore, the flow blockage can only occur at the inlet of the channel. However, debris carried by the coolant may be of different kinds of forms and shapes. In this paper, it is assumed that the blockage is caused by a zero-thickness plane, and the area of the plane is 90% of the blocked channel inlet size. However, it should be emphasized that blockages of any given shapes can be modeled and investigated by means of CFD methods. The investigated domain is depicted in Fig. 3, in which, an upper plenum and a lower plenum are added to give reasonable inlet and outlet conditions conveniently. During the initial calculation of the simulation, the blockage plane is stationary at the upper plenum. When transient calculation begins, the blockage plane starts to move with a speed of 2.3 m/s, which is the main stream velocity in the upper plenum. Finally, the blockage plane arrives at the inlet of the blocked channel and stops there. During this process, the heat generation rate of the fuel is set to a chopped cosine function with the effective length set to 0.6 m, and the extrapolated length is 0.76 m (Salama and El-Morshedy, 2012). Due to the potential big rise in temperature, the density, specific heat, thermal conductivity and viscosity of the coolant are all considered to change with temperature at the form of quartic polynomials.
accurate near-wall treatment is the wall functions approach. Wall functions are a series of semi-empirical formulas, linking the viscosity-affected region and the main stream without solving the viscosity-affected inner region. In these functions, there is a significant non-dimensional parameter y* defined as follows: 1/4 1/2
y∗ =
C kP yP
Direct Numerical Simulation (DNS) of turbulent flows is too computationally expensive to be applied to complex geometries. Reynolds-averaging and filtering are the alternative methods to simulate turbulent flow, but the latter also requires a significant amount of computational resources. The Reynolds-averaged Navier–Stokes (RANS) equations are time-averaged equations of motion for fluid flow, which are computationally less expensive to solve and the results can be acceptable for engineering. The RANS technique therefore is widely adopted for practical engineering applications. But RANS equations introduce additional unknown variables (Reynolds Stress Tensor) and turbulence models are established in order to achieve the problem closure and determine these variables in terms of known quantities. The realizable k–ε mode is applied in this study, and the modeled transport equations for k and ε in the realizable k–ε model (ANSYS, Inc., 2009) are
3.3. Meshing strategies Structured mesh is used to discretize the simulated domain which is of high aspect ratio. The combination of wall functions and structured mesh reduces the mesh quantity significantly, making it reasonably accurate and economical for the transient calculation. To simulate the moving plane, dynamic mesh technique must be applied. The layering type dynamic mesh in ANSYS FLUENT is designed for structured mesh, allowing motion along a specific straight line, which makes it suitable for the current simulation. Moving zones and stationary zones in the layering method are connected by sliding interfaces, through which data from one zone are interpolated and transferred to the other. Errors occur during
⎧
t ∂k ∂ ∂ ∂ ⎪ ⎪ ku = + + Gk + Gb − εA − YM − Sk + (k) ⎪ j ⎨ ∂t k ∂xj ∂xj ∂xj
⎪ ∂ ∂ ∂ ε2 ε t ∂ε ⎪ ⎪ εu = + + C1 Sε − C2 + √ + C1ε C3ε Gb + Sε ⎩ (ε) j ∂xj
where
C1 max 0.43,
∂xj
+5
,
k =S , ε
ε
S=
∂xj
k+
2Sij Sij
In these equations, Gk represents the generation of turbulence kinetic energy due to the mean velocity gradients. Gb is the generation of turbulence kinetic energy due to buoyancy. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. C1ε , C2ε , and C3ε are constants. k and ε are the turbulent Prandtl numbers for k and ε, respectively. Sk and Sε are source terms. On the other hand, like other k–ε models, the realizable k–ε model is primarily valid for turbulent flows in the main stream. Special attention, therefore, needs to be paid to simulate turbulent flows with confining walls through the use of an appropriate near-wall treatment. An economical, robust, and reasonably
(2)
where P is the near-wall node; C is a constant for standard and RNG k–ε models, but is a variable in realizable k–ε model; kP is turbulence kinetic energy at P; yP is distance from point P to the wall; is kinematic viscosity of the fluid. Among these functions, standard wall functions are popular for its good performances in a wide range of wall-bounded flows. But in unsteady flows, in which the mass flow rate changes with time, y* will change. Especially when the mass flow rate decays dramatically, y* may become smaller than 11.225, which is the intersection of the logarithmic and the linear near-wall profile, and this will lead to deterioration of standard wall functions. When flow blockage occurs, the mass flow rate of the channel will decrease significantly, resulting from the increase in flow resistance. Standard wall functions therefore cannot be used to simulate the channel flow blockage transient. On the basic assumption that the viscous sub-layer is so thin that it has the same boundary with the wall, scalable wall functions (ANSYS, Inc., 2012) overcome the drawback by introducing a limiter in the y* calculations such that y˜ ∗ = max(y∗ , 11.225), and y* used for any formula in standard wall functions is replaced by y˜ ∗ . This feature of scalable wall functions makes them suitable for the present study.
3.2. CFD techniques and turbulence models
∂t
179
ε
(1)
k
interpolations, which also require extra computing time. Therefore, the size of moving zone should be minimized to lessen interface size, which can reduce numerical errors and accelerate the simulation. In this study, to guarantee the accuracy of heat transfer calculation between claddings and coolants, all interfaces are set between fluid zones and fluid zones only. Obviously, the zone between the initial plane location and the inlet of the blocked channel is being compressed during the calculation. However, this zone becomes zero-volume when the plane finally arrives at the inlet, leading to change in topology structure of domains, which is not allowed by the solver. To avoid the alteration in topology structure, an L-grid division is conducted in the blocked channel block, as illustrated in Fig. 4. The short branch near the inlet and the compressed zone in the upper plenum are combined to a new compressed zone, and the long branch beside the fuel plate is set to boundary layer mesh.
180
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
Fig. 3. Modeling of the investigated domain. (a) Front view. (b) Oblique view.
According to Best Practice Guidelines for the use of CFD in NRA (Mahaffy and Farm, 2010), grid sensitivity analysis should be conducted to make sure that the simulation is independent of mesh size. Velocity magnitude and relative static pressure of the centroid of channel-1 are compared, as shown in Table 3. Several grids have been examined; when the node number is increased to 1797,576 and 2205,604, both velocity magnitude and relative static pressure differ by less than 0.8%. When a more dense mesh with 3861,684 nodes is taken into consideration, the differences in the two finer meshes are less than 0.05%. Therefore, the 2205,604 nodes mesh is chosen for further calculations, a cross-sectional view of the mesh is depicted in Fig. 5. As advised by ANSYS FLUENT, the time step size of transient calculations is chosen from 5 × 10−6 s to 4 × 10−4 s during different simulation periods to satisfy the convergence requirements. Fig. 4. The L-grid division strategy.
Therefore, even though the plane arrives at the inlet, the new compressed zone still has a size of the short branch, which is accepted by the solver. 4. Results and discussions Because of its good performance in unsteady simulations, the PISO method is employed for pressure–velocity coupling. The gradient term and the pressure term are calculated by interpolation using the least squares cell based scheme and the PRESTO! scheme, respectively. The second order upwind scheme is used for all the variables. The first order implicit scheme is applied for the discretization of the transient formulation. According to Salama and El-Morshedy (2012), the average inlet velocity of the upper plenum is set to 2.3 m/s. In order to define an accurate flow boundary for the inlet, fully developed velocity and turbulence profiles are added to the inlet, and these quantities are fixed during transient calculations. For the lower plenum outlet, the boundary type is set to pressure-outlet. To summarize, main boundary and initial conditions are listed in Table 2.
4.1. Validation consideration According to Best Practice Guidelines for the use of CFD in NRA (Mahaffy and Farm, 2010), model validation should be performed before application to real case studies. However, due to the compact structure and extremely high heat flux, it is difficult to conduct accurate experiment on flow blockage accidents. Besides, transient experiment data are needed to validate results of the current study. Therefore, it is hard to validate results of the present study. In order to compare with the CFD results, system analysis code Relap5(Mod3.4) is used to calculate the steady results in both
Table 2 Boundary and initial conditions. Inlet of the upper plenum Outlet of the lower plenum Heat generation rate of fuels Moving speed of the blockage plane Duration of the plane moving transient
Fully developed velocity-inlet, average velocity is 2.3 m/s Pressure-outlet Chopped cosine function 2.3 m/s 0.085 s
Fig. 5. Cross-sectional view of the mesh.
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
181
Table 3 Grid sensitivity analysis. Mesh
Cells
Nodes
Velocity
1 2 3 4
1019,208 1470,198 1868,337 3246,042
1316,064 1797,576 2205,604 3861,684
4.126431 4.137421 4.140962 4.139675
non-blocked and blocked states of the domain, the corresponding model is shown in Fig. 6. System codes are widely used for their simplicity and relatively low requirements on computing resources, however there are some disadvantages in the use of such system analysis codes. As mentioned before, empirical formulas are used in one-dimensional codes. For instance, during the calculating process, Tw (z), the coolant temperature, is calculated according to the law of energy conservation first, and then Tcs (z), the cladding surface temperature, is calculated as follows: Tcs (z) = Tw (z) +
q (z) h (z)
(3)
where q (z) is the heat flux of the cladding surface, and h(z) is the convection heat transfer coefficient, which is derived from empirical formulas. However, there are many different empirical formulas that can be used for h(z) calculation, therefore different Tcs (z) are predicted. For instance, in Salama (2011), the Dittus–Bolter equation (Winterton, 1998) provides consistent results with the CFD simulation on one channel, while in Salama et al. (2015), the Sieder and Tate (1936) equation gives consistent predictions with the CFD calculation on more channels. Therefore, it is the choice of empirical equations that determines the Tcs (z) calculation. On the other hand, different modeling strategies toward the assembly structure may also influence the result. As depicted in Fig. 2, along the width direction, there are active zones and inactive zones in fuel plates. In CFD approaches, models are precisely built according to the assembly structure, while during Relap5 modeling, due to the limitation of the code, fuels and claddings are assumed
Fig. 6. The Relap5 model.
Error (%) 0.266 0.086 −0.031
Relative static pressure 10,241.457 10,298.699 10,377.161 10,382.158
Error (%) 0.559 0.762 0.048
to have the same width. This simplification changes the plate size, and the temperature field of the plate is changed consequently. Due to the differences in modeling and solving methods, like other system codes, Relap5 results are only used to compare with CFD predictions. For the non-blocked scenario, the coolant, cladding surface and fuel temperature profiles along the channel-1 in both simulations are depicted in Fig. 7. For the CFD approach, temperature is calculated for each volume. In order to compare with the onedimensional results, temperatures should be averaged in every cross section along the channel. However, proper weight should be selected when calculating the average temperature of the coolant. The area-weighted temperature and the z-velocity-weighted temperature represent the average temperature of the cross section statically and dynamically, respectively. As shown in Fig. 7, for the non-blocked scenario, these two kinds of average values are quite close. When the comparison is done between results calculated by CFD and Relap5, good closeness is also obtained for coolant temperatures. Though the cladding surface and fuel temperatures calculated by CFD method are higher than those obtained by Relap5, differences in the models can explain differences in results. 4.2. Velocity transients Like system analysis, mass flow rate of a channel can be obtained by CFD method. Besides, there many more details of the simulated domain can be obtained through CFD method. 4.2.1. Mass flow rate transients It is obvious that when the plane stops at the inlet of the blocked channel, flow resistance increases significantly, resulting in fast changes in mass flow rate of the investigated channels, especially in the blocked one. In this paper, for the sake of comparisons, all mass flow rates are converted to those in the half-channel correspondingly. Transient mass flow rates in each channel are shown in Fig. 8. When the plane is moving with the coolant, mass flow rate varies little. While when the plane arrives at the inlet, mass flow rate in
Fig. 7. Temperature profiles along channel-1(non-blocked scenario).
182
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
Fig. 8. Transitional mass flow rate in different channels.
the blocked channel drops dramatically, with those in the other two channels both increasing. The plane stops at 0.085 s after the simulation, when flow time comes to 0.1 s, the rapid change in mass flow rate of all the channels converts to a quite slow one. Therefore, the redistribution of the mass flow rate can be considered to complete in only 0.015 s. However, the result of the redistribution is different from that acquired by Salama and El-Morshedy (2012), which suggests that flow in channel-3 almost remains unchanged in different blockage ratio scenarios. The difference results from the different blockage modeling methods. In their investigation, the blockage is caused by buckling of the plate. The sizes of all the inlets and outlets are consistent with those in the non-blocked scenario, so the mass flow rate redistribution is mainly driven by flow resistances along different channels. Considering there is a flow area expansion in channel-2 corresponding to the flow area contraction in channel-1, and the flow area in channel-3 is unchanged, therefore, the flow resistance of channel-3 is larger than that in channel-2, but lower than that in channel-1. Being the channel with the smallest flow resistance, channel-2 has the highest mass flow rate in blocked scenarios. In the present study, the blockage occurs at the inlet of the flow channel, and it is the local flow resistances of different inlets that determine the mass flow rate redistribution. Compared to channel3, channel-2 is closer to the blockage, and is more affected by the blockage. Therefore the mass flow rate of channel-2 is lower than that of channel-3. The dissimilarity in mass flow rate redistribution prediction shows that the blockage form affects the flow field in the assembly strongly. For blockages caused by buckling, more coolant flows through channel-2 than that in inlet flow blockage cases, so the fuel plate adjacent to the channel-1 is cooled better. Therefore, the inlet flow blockage may lead to a severer consequence than the buckling blockage does. 4.2.2. Reversed flow around the outlet of the blocked channel A small reversed flow region occurs around the outlet of the blocked channel at about 0.1 s. Fig. 9 shows the z-velocity distribution at the outlet of the blocked channel. This tiny region then becomes larger, after that reaches a relatively steady size. However, it is hard to explain the reversed flow around the outlet, especially for channels through which coolant flows downward. As will be seen later, due to the presence of the blockage, quite a large area of coolant in the blocked channel are influenced, and flow field is not fully developed in channel-1. In addition, when coolant from channels flows to the lower plenum, the flow field is quite
Fig. 9. Profiles of z-velocity on the outlet of the blocked channel.
complex, even in normal scenario. For the blocked condition, mass flow rate in the blocked channel is much lower than other channels, which makes the flow field more complicated. Therefore, the large mass flow rate difference between the blocked and unblocked channels should be another possible reason. 4.2.3. Flow field in the blocked channel The flow field in channel-1 changes rapidly during the calculation, especially when the plane is about to arrive at the inlet and a very short time after its arrival. As shown in Fig. 10, because of the shrink of the flow area and the inertia of the upstream coolant, a high speed jet-flow occurs when the plane stops at the inlet, which leads to backflow behind the plane at the same time. The maximum velocity in the blocked channel rockets to as high as 50 m/s, and the corresponding maximum backflow velocity also reaches its peak, which is almost 13 m/s, and maximum and minimum values of them during the simulation are depicted in Fig. 11. Again, the fast change takes place in about 0.015 s. A problem therefore arouses, that the high speed jet-flow may threaten the assembly structure. As depicted in Fig. 12, the high speed region spreads instantly from the inlet to the downstream of the blocked channel after the formation of the blockage. Due to the viscosity of the coolant, a larger region is affected by the blockage during the spread, with the maximum velocity decaying sharply both in flow direction and backflow direction. Under the influence of both the high speed region and the small reversed flow region around the outlet, a new backflow region occurs around the middle part of the channel, in about 1.5 s after the calculation. After that, the new region grows larger, and the velocity field structure becomes stable. 4.3. Temperature transients As a consequence of the redistribution of the flow field, the temperature field of the domain changes with time. The maximum temperatures of the cladding and the fuel during calculation are depicted in Fig. 13. Since the melting point of the cladding is 855.15 K, the peak temperatures during calculation are far from their corresponding melting points, therefore the integrity of the fuel plate is maintained in this simulation. The boiling point
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
Fig. 10. Profiles of flow direction velocity in the blocked channel (flow time form 0.085 s to 0.0999 s).
of the coolant is 388.32 K under the operating pressure, therefore no boiling of the coolant occurs during the transient. The peak temperature around 0.8 s can be explained as follow. As the power distribution is a chopped cosine function along the channel, the peak always appears around the middle of the plate. The mass flow rate of the blocked channel decreases sharply after the formation of the blockage. Temperatures, except those in high speed regions, rise dramatically due to the bad heat transfer condition. As mentioned before, the high speed regions will spread to the
Fig. 11. The maximum and minimum z-velocity in the blocked channel.
183
Fig. 12. Profiles of flow direction velocity in the blocked channel (flow time form 0.299 s to 1.998 s).
down stream; at about 0.8 s, the middle part is involved, as shown in Fig. 12. High speed coolant with lower temperature enhances heat transfer in this region. As a consequence, the maximum temperature starts to drop at about 0.8 s. The previous discussion also indicates that, the high speed regions will slow down and then a backflow region will occur around the middle part of the channel. This leads to the rebound of the highest temperatures of the cladding and fuel. This explains the closeness between the time of the formation of the backflow region and that of the start of the bouncing.
Fig. 13. Maximum temperatures during calculation.
184
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
Fig. 14. Temperature and velocity field in the blocked channel.
4.4. Steady results of blocked scenario 4.4.1. Details calculated by CFD approach As the calculation continues, both flow field and temperature field of the investigated domain gradually transform to a new steady distribution. Fig. 14 shows the temperature and velocity fields in the blocked channel. A big vortex occurs behind the blockage, extending to the middle part. There is another vortex below the big one, but with small scale. As mentioned before, there is a reversed flow region around the outlet. Actually, there are two tiny vortex regions below the outlet. The highest coolant temperature does not occur in the middle of the channel where is the hottest area of the plate, as shown in Fig. 15, because coolant around the middle part is of high speed. It is the region where the largest vortex occurs that the maximum coolant temperature appears, because of the relatively slow circle motion of the coolant. For the cladding, besides the highest temperature in the middle, another small area is of a relatively high temperature a little past the middle, because of the corresponding vortex. 4.4.2. Comparisons with Relap5 To simulate the blocked scenario in Relap5, the flow area of the single junction, which is used to connect the upper plenum and the blocked channel, is reduced to 10% of its initial value. The mass flow rates in all channels are depicted in Fig. 16, in which mass flow rate in channel-1 calculated by FLUENT is 19.8% higher than that simulated by Relap5. This large difference results from the quite small mass flow rates in channel-1; actually, the difference between two predictions is about 0.00812 kg/s, which is only 1.2% of the total mass flow rate. Close predictions are obtained for mass flow rates in channel-2 and channel-3, and the dissimilarity is 1.8% and 0.5%, respectively. However, different profiles for channel-2 and channel-3 in two treatments appear.
Fig. 15. Cladding temperature field of the blocked channel.
In CFD approach, as the adjacent channel of the blocked channel, channel-2 is more affected by the blockage than channel-3. More coolants, therefore, pass through channel-3 than channel-2. While in Relap5, channel-2 and channel-3 are not recognizable in geometric positions. Therefore, compared to channel-3, a larger mass flow
Fig. 16. Mass flow rate comparisons.
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
185
Fig. 17. Power redistribution comparisons. Fig. 19. Temperature profiles along channel-1 (blocked scenario).
rate for channel-2 is predicted, as a sequence of the higher heating power in channel-2, which will be discussed in the following. Because of the big mass flow rate differences of the blocked channel and its adjacent one, power generated in the fuel plate between them redistributes, as shown in Fig. 17, and quite similar profiles are obtained. Compared to coolants in channel-2, those in channel-1 are of very low Re number. Therefore, the heat transfer coefficients on both sides of fuel Plate 1are quite different. For the wall side toward channel-1, the heat transfer coefficient is much lower than the coefficient on the wall side toward channel-2. Consequently, huge asymmetric heat transfer occurs in fuel Plate 1. While for fuel Plate 2, due to the little difference in mass flow rates of channel-2 and channel-3, the distributions are almost symmetric for both simulations. Fig. 18 shows temperature profiles of different outlet in both normal and blocked cases. The outlet temperature of the blocked channel dramatically rises due to the huge drop in mass flow rate of the channel. As mentioned above, channel-1 are almost equally heated in both calculations, however a higher mass flow rate of channe-1 is predicted by CFD approach. The outlet temperature, therefore, is a little lower than that calculated by Relap5 codes. Due to the total mass flow rate is fixed in both simulations, the mass flow rates in channel-2 and channel-3 both become larger. Channel-2 gets more energy form fuel Plate 1, therefore, the outlet temperature rises a little. For channel-3, the energy source is almost
same with the normal case, increase in mass flow rate leads to a slight drop in outlet temperature. Both the comparisons of heat transfer ratios of fuel plates and outlet temperature profiles indicate that, when considering heat transfer effects of the flow blockage, only the blocked channel and plates that constitute the blocked channel are mostly affected. In another word, the blocked channel and its adjacent ones are significant for investigation. The coolant, cladding surface and fuel temperature profiles along the flow direction are compared in Fig. 19. For Relap5 treatment, except for the detailed values, tendency of profiles are quite similar to those in the normal scenario. This indicates that flow details in the blocked channel are totally ignored. For CFD approach, three-dimensional results are converted to one-dimensional profiles, and the shapes of the cladding surface temperature and fuel temperature profiles are different from those in the non-blocked case, especially around the middle of the plate. This difference indicates that power redistribution occurs in the flow direction, which is not considered in Relap5 treatment. Since a relatively higher temperature is predicted in the CFD approach, the z-velocity-weighted temperature is lower than that calculated by Relap5 with a similar profile shape. However, the area-weighted temperature is quite irregular due to the complex flow and temperature fields. In Relap5, outlet temperature is the highest temperature of the coolant, as shown in Fig. 18, and it is higher than the outlet temperature that is predicted by CFD method. However, the value is still lower than the maximum area-weighted temperature in CFD approach. Therefore, due to the neglect of flow details, conservative perdition in the mass flow rate cannot guarantee a conservative temperature result. 5. Conclusions and outlooks
Fig. 18. Temperature profiles of different outlet in both normal and blocked cases.
Using dynamic mesh technique, the 90% inlet flow blockage of an MTR plate-type assembly is investigated by three dimensional CFD methods. Relap5 simulation is conducted for comparison. No boiling of the coolant is predicted, and both CFD and Relap5 simulations indicate that, when concerning heat transfer effects, the blocked channel and its adjacent ones are mostly affected. Besides, quite similar asymmetry power redistributions are predicted by these two methods. The transient CFD simulation shows that the development of flow blockage transients is so fast that the rapid redistribution of mass flow rates occurs in only 0.015 s after the formation of the blockage. Besides, more significant details of the domain are captured by the CFD method, such as the jet-flow and vortexes in the blocked channel, and the irregular temperature profiles caused by the blockage.
186
W. Fan et al. / Nuclear Engineering and Design 292 (2015) 177–186
The present study shows advantages of using transient CFD method to investigate flow blockages which are caused by debris. However, relevant experiments are needed for validations. As mentioned before, it is quite difficult to conduct experiment. Using current results, targeted measurements may make the experiment more practical. Besides, further simulations can be conducted on larger domain and toward different forms of debris. Acknowledgements The authors would like to thank the Science and Technology on Reactor System Design Technology Laboratory of Nuclear Power Institute of China for financial and experiential support of this work. This study is also supported by National Natural Science Foundation of China (Grant no. 11305169). References Adorni, M., Bousbia-Salah, A., Hamidouche, T., et al., 2005. Analysis of partial and total blockage of a single fuel assembly of an MTR research reactor core. Ann. Nucl. Energy 32, 1679–1692. ANSYS, Inc., 2009. ANASYS FLUENT 12.0 Theory Guide. ANSYS, Inc., Canonsburg, PA, USA. ANSYS, Inc., 2012. ANASYS 14.5 Help. ANSYS, Inc., Canonsburg, PA, USA. Bestion, D., 2004. Recommendation on use of CFD codes for nuclear reactor safety analysis. In: EVOL-ECORA-D14. Bousbia-salah, A., Hamidouche, T., 2005. Analysis of the IAEA research reactor benchmark problem by the RETRAC-PC code. Nucl. Eng. Des. 235, 661–674. Chatzidakisa, S., Ikonomopoulosa, A., Ridikasb, D., 2013. Evaluation of RELAP5/MOD3 behavior against loss of flow experimental results from two research reactor facilities. Nucl. Eng. Des. 255, 321–329. El-Messiry, A.M., 2000. Cobalt irradiation box ejection accident of ETRR-2. Nucl. Eng. Des. 198, 287–293. El-Morshedy, S.E., 2011. Prediction, analysis and solution of the flow inversion phenomenon in a typical MTR-reactor with upward core cooling. Nucl. Eng. Des. 241, 226–235.
Hainouna, A., Ghazia, N., Mansour Abdul-Moaizb, B., 2010. Safety analysis of the IAEA reference research reactor during loss of flow accident using the code MERSAT. Nucl. Eng. Des. 240, 1132–1138. Hamidouche, T., Bousbia-salah, A., 2010. Assessment of RELAP5 point kinetic model against reactivity insertion transient in the IAEA 10MW MTR research reactor. Nucl. Eng. Des. 240, 672–677. Höhne, T., Krepper, E., Rohde, U., 2010. Application of CFD codes in nuclear reactor safety analysis. Sci. Technol. Nucl. Instal. 2010, 1–8. IAEA, 1980. IAEA Research Reactor Core Conversion from the use of high-enriched uranium to the use of low enriched uranium fuels Guidebook. In: IAEA-TECDOC233. IAEA. IAEA, 1992. Research reactor core conversion guidebook, volume 3: analytical verification (Appendices G and H). In: IAEA-TECDOC-643. IAEA. Lu, Q., Qiu, S., Su, G.H., 2009. Flow blockage analysis of a channel in a typical material test reactor core. Nucl. Eng. Des. 239, 45–50. Mahaffy, J., Farm, W.S., 2010. Development of best practice guidelines for CFD in nuclear reactor safety. Nucl. Eng. Des. 42, 377–381. Salama, A., 2011. CFD investigation of flow inversion in typical MTR research reactor undergoing thermal–hydraulic transients. Ann. Nucl. Energy 38, 1578–1592. Salama, A., 2012. CFD analysis of fast loss of flow accident in typical MTR reactor under going partial and full blockage: the average channel scenario. Prog. Nucl. Energy 60, 1–13. Salama, A., El-Amin, M.F., Sun, S., 2015. Three-dimensional, numerical investigation of flow and heat transfer in rectangular channels subject to partial blockage. Heat Transfer Eng. 36, 152–165. Salama, A., El-Morshedy, S.E., 2011. CFD simulation of the IAEA 10 MW generic MTR reactor under loss of flow transient. Ann. Nucl. Energy 38, 564–577. Salama, A., El-Morshedy, S.E., 2012. CFD simulation of flow blockage through a coolant channel of a typical material testing reactor core. Ann. Nucl. Energy 41, 26–39. Sieder, E.N., Tate, G.E., 1936. Heat transfer and pressure drop of liquids in tubes. Ind. Eng. Chem. 28, 1429–1435. Sim, T.M., Tabor, W.H., 1964. Report on Fuel-plate Melting at the Oak Ridge Research Reactor: ORNL-TM-627. ORNL, USA, Oak Ridge. Varvayanni, M., Stakakis, E., Catsaros, N., Antonopoulos-Domis, M., 2005. Void induced reactivity in a mixed MTR core. Nucl. Eng. Des. 235, 855–865. Winterton, R.H.S., 1998. Where did the Dittus and Boelter equation come from? Int. J. Heat Mass Transfer 41, 809–810. Yadigaroglu, G., Andreani, M., Dreier, J., Coddington, P., 2003. Trends and needs in experimentation and numerical simulation for LWR safety. Nucl. Eng. Des. 221, 205–223.