Development of an Analytical Method on Water-Vapor Boiling Two-Phase Flow Characteristics in BWR Fuel Assemblies Under Earthquake Condition

an analytical method for boiling two-phase flow in a fuel assembly under earthquake conditions by improving ACE-3D and shows how the three-dimensional behavior of boiling two-phase flow under these conditions is evaluated by the improved ACE-3D. fuel assembly under oscillation conditions was evaluated using a simulated fuel assembly excited by oscillation acceleration. On the basis of this evaluation, it was confirmed that void fraction fluctuation This book presents a comprehensive review of studies in nuclear reactors technology from authors across the globe. Topics discussed in this compilation include: thermal hydraulic investigation of TRIGA type research reactor, materials testing reactor and high temperature gas-cooled reactor; the use of radiogenic lead recovered from ores as a coolant for fast reactors; decay heat in reactors and spent-fuel pools; present status of two-phase flow studies in reactor components; thermal aspects of conventional and alternative fuels in supercritical water?cooled reactor; two-phase flow coolant behavior in boiling water reactors under earthquake condition; simulation of nuclear reactors core; fuel life control in light-water reactors; methods for monitoring and controlling power in nuclear reactors; structural materials modeling for the next generation of nuclear reactors; application of the results of finite group theory in reactor physics; and the usability of vermiculite as a shield for nuclear reactor.


Introduction
Safe operation of nuclear reactors under earthquake conditions cannot be guaranteed because the behavior of thermal fluids under such conditions is not yet known. For instance, the behavior of gas-liquid two-phase flow during earthquakes is unknown. In particular, fluctuation in the void fraction is an important consideration for the safe operation of a nuclear reactor, especially for a boiling water reactor (BWR). The void fraction in the coolant is one of the physical parameters important in determining the thermal power of the reactor core, and fluctuations in the void fraction are expected to affect the power of the plant.
To evaluate fluctuation in the void fraction, numerical simulation is the most effective and realistic approach. In this study, we have developed a numerical simulation technique to predict boiling two-phase flow behavior, including fluctuation in the void fraction, in a fuel assembly under earthquake conditions.
In developing this simulation technique, we selected a three-dimensional two-fluid model as an analytical method to simulate boiling two-phase flow in a fuel assembly because this model can calculate the three-dimensional time variation in boiling two-phase flow in a large-scale channel such as a fuel assembly while incurring only a realistic computational cost. In addition, this model has been used to successfully predict the void fraction for a steady-state boiling two-phase flow simulation (Misawa, et al., 2008). We expect that the development of the boiling two-phase flow analysis method for a fuel assembly under earthquake conditions can be achieved by improving the three-dimensional two-fluid model analysis code ACE-3D (Ohnuki, et al., 2001;Misawa, et al., 2008), which has been developed by the Japan Atomic Energy Agency. This paper describes an analytical method for boiling two-phase flow in a fuel assembly under earthquake conditions by improving ACE-3D and shows how the three-dimensional behavior of boiling two-phase flow under these conditions is evaluated by the improved ACE-3D. www.intechopen.com

Development of an analytical method of boiling two-phase flow in a fuel assembly under earthquake conditions 2.1 Overview of ACE-3D
ACE-3D has been developed by the Japan Atomic Energy Agency (JAEA) to simulate watervapor or water-air two-phase flows at subcritical pressures. The basic equations of ACE-3D are shown below.
Mass conservation for vapor and liquid phases: Here, t represents time; x, a spatial coordinate; U, velocity; P, pressure; g, acceleration due to gravity; e, internal energy; and ρ, density. Subscripts g and l indicate vapor and liquid phases, respectively. Subscripts i and j indicate spatial coordinate components. If the subscripts of the spatial coordinate components are duplicative, the summation convention is applied to the term. The terms q w and q int indicate wall heat flux and interfacial heat flux, respectively, and h sat represents saturation enthalpy. Summation of volume ratios α g and α l is equal to one. Γ + , Γ -, and Γ represent vapor generation rates. Γ + in Eq.
(3) is equal to Γ if Γ is positive and is zero if Γ is negative. On the other hand, Γis equal to Γ if Γ is negative and is zero if Γ is positive. τ g,ij and τ l,ij in Eqs.
Summation of interface stress, M int , of the vapor and liquid phases is equal to zero. The interface stress is determined by correlations between turbulent diffusion force (Lehey, 1991), lift force (Tomiyama, 1995), interface friction force, and bubble diameter (Lilies, 1988). The lift force is calculated as follows.
In Eq. (7), σ is the surface tension coefficient; E t , Eotvos number; Re bubble , the bubble Reynolds number; and D b , bubble diameter. c lift0 in Eq. (7) describes the effect of shear flow and is a positive value. c wake describes the effect of bubble deformation and significantly depends on Eotvos number of bubble diameter. The coefficient of lift force c lift is estimated by the summation of c lift0 and c wake . If the Eotvos number of bubble diameter is small, c lift is positive. However, if the Eotvos number of bubble diameter increases, c lift is negative. Therefore, the lift force corresponding to a small bubble diameter acts in a direction opposite to that corresponding to a large bubble diameter.
Bubble diameter is evaluated by the following equations.
where σ is the surface tension coefficient; We, a critical Weber number of 5.0; and D bmax , the maximum bubble diameter, which is an input parameter. Bubble diameter, D b , is estimated by linear interpolation of small bubble diameter, D bubble , and slug diameter, D slug , with a coefficient, X slug , which is dependent on the void fraction. Therefore, D b increases with an increase in the void fraction.
In ACE-3D, the two-phase flow turbulent model based on the standard k-ε model (Bertodano, 1994) is introduced below. Turbulent energy conservation for vapor and liquid phases: Turbulent energy dissipation rate conservation for vapor and liquid phases: The basic equations represented in Eqs.
(1) to (12) are expanded to a boundary fitted coordinate system (Yang, 1994). ACE-3D adopts the finite difference method using constructed grids, although it is difficult to construct fuel assembly geometry by using only constructed grids. Therefore, a computational domain, which consists of constructed grids, is regarded as one block, and complex geometry such as that of a fuel assembly is divided into more than one block. An analysis of complex geometry can be performed by a calculation that takes the interaction between blocks into consideration. Parallelization based on the Message Passing Interface (MPI) was also introduced to ACE-3D to enable the analysis of large-scale domains.
From a previous experiment in which the three-dimensional distribution of the void fraction (vapor volume ratio) in a tight-lattice fuel assembly was measured (Kureta, 1994), it is known that the vapor void (void fraction) is concentrated in the narrowest region between adjacent fuel rods near the starting point of boiling and as the elevation of the flow channel increases, the vapor void spreads over a wide region surrounding the fuel rods. This tendency of the vapor void to redistribute has been described by a past analysis using ACE-3D (Misawa, et al., 2008). Therefore, it is confirmed that the boiling two-phase flow analysis can be carried out using ACE-3D under steady-state conditions and with no oscillation.

Improvement of three-dimensional two-fluid model for earthquake conditions
In order to simulate the boiling two-phase flow in a fuel assembly under earthquake conditions, it is necessary to consider the influence of structural oscillation of reactor equipment on boiling two-phase flow. If the coordinate system for an analysis is fixed to an oscillating fuel assembly under earthquake conditions, it can be seen that a fictitious force acts on the boiling two-phase flow in the fuel assembly. Therefore, a new external force term, f i , which simulates the acceleration of oscillation, was added to the momentum conservation equations (Eqs. (3) and (4)).
We assume that the analysis of boiling two-phase flow in a fuel assembly under earthquake conditions can be performed by using time-series data as an input if the time-series data of oscillation acceleration can be obtained from structural analysis results for a reactor (Yoshimura, et al., 2002) or if the measurement data of actual earthquakes can be obtained by seismographic observation.
In order to apply this improved method to the analysis of boiling two-phase flow in a fuel assembly under earthquake conditions, it is necessary to confirm that the simulation of boiling two-phase flow under oscillation conditions can be performed using the interface stress models shown in the preceding section; these stress models are empirical correlations and are based on experimental results under steady-state conditions. In the case of boiling two-phase flow analysis under oscillation conditions, these interface stress models may cause instability in simulation results.
In addition, it is necessary that large-scale analysis be performed within limited computable physical time and that it be consistent with the time-series data of oscillation acceleration obtained from the results of structural analysis in a reactor or with the measurement data from actual earthquakes. In structural analysis in a reactor (Yoshimura, et al., 2002), the minimum time interval of the analysis is limited to 0.01 s (100 Hz). Seismographic observation is also frequently performed with a sampling period of 100 Hz. If a highfrequency oscillation acceleration of over 100 Hz influences boiling two-phase flow in a fuel assembly, boiling two-phase flow analysis, which is consistent with the structural analysis in a reactor, cannot be performed. Therefore, it is necessary to evaluate the highest frequency necessary for this improved method to be consistent with the time-series data of oscillation acceleration.
A computable physical time of about 1 s is preferred for the boiling two-phase flow analysis in a fuel assembly because this analysis requires a large number of computational grids in order to simulate a large-scale domain such as a fuel assembly. If the results of the boiling two-phase flow analysis show quasi-steady time variation for long-period oscillation acceleration, it is not efficient to perform the analysis with a computable physical time span longer than the long period. Effective analysis can be performed if the analysis with a time span subequal to the shortest period of oscillation acceleration, for which the boiling twophase flow shows quasi-steady time variation, by extracting earthquake motion at any time during the earthquake. Therefore, it is necessary to evaluate the shortest period of oscillation acceleration for which the boiling two-phase flow shows quasi-steady time variation.
The boiling two-phase flow was simulated in a heated parallel-plate channel, which is a simplification of a single subchannel in a fuel assembly. The channel was excited by vertical www.intechopen.com and horizontal oscillation to simulate an earthquake in order to confirm that the boiling twophase flow simulation can be performed under oscillation conditions.
In addition, the influence of the oscillation period on the boiling two-phase flow behavior in a fuel assembly was investigated in order to evaluate the highest frequency necessary for the improved method to be consistent with the time-series data of oscillation acceleration and the shortest period of oscillation acceleration for which the boiling two-phase flow shows quasi-steady time variation.

Confirmation of stability of the boiling two-phase flow simulation under oscillation conditions
The parallel-plate channel, which simulates a single subchannel in a fuel assembly, was adopted as the computational domain, as shown in Fig. 1. Both plates were heated with a uniform heat flux of 270 kW/m 2 . The single-phase water flows into the parallel-plate channel vertically from the inlet. The hydraulic diameter and heated length of the computational domain are equal to those of the single subchannel in the fuel assembly of a current BWR. The outlet pressure of 7.1 MPa, the inlet velocity of 2.2 m/s, and the inlet temperature of 549.15 K also reflect the operating conditions in the current BWR. The adiabatic wall region is set up on the top of the heated region in order to eliminate the influence of the outlet boundary condition. In this analysis, the maximum bubble diameter in Eq. (8) is set to the channel width of 8.2 mm.
First, an analysis was performed without applying oscillation acceleration. After a steady boiling flow was attained, oscillation acceleration was applied. The time when the oscillation acceleration was applied is regarded as t = 0 s.

Fig. 1. Computational domain
Two cases of oscillation acceleration, in the vertical direction (Z axis) and in the horizontal direction (X axis), were applied. In both cases, the oscillation acceleration was a sine wave with a magnitude of 400 Gal and a period of 0.3 s, as shown in Fig. 2. The magnitude and period of the oscillation accelerations were taken from actual earthquake acceleration data.
www.intechopen.com   the oscillation acceleration in boiling flow. In this case, the liquid phase is driven by the oscillation acceleration because the influence of the oscillation acceleration is relatively large owing to a high liquid density; on the other hand, the vapor phase is driven by the horizontal pressure gradient because the influence of the oscillation acceleration is less than that of the horizontal pressure gradient owing to the low vapor density. This explains why the vapor velocity and the void fraction moved in a direction opposite to that of the oscillation acceleration. A comparison between Fig. 4 and Fig. 6 indicates that the magnitude of the void fraction fluctuation for the horizontal oscillation acceleration case was greater than that for the vertical oscillation acceleration case at any vertical position.
It can therefore be confirmed that the fluctuation of the void fraction with the same period as the oscillation acceleration can be calculated in the case of both horizontal and vertical oscillation acceleration.

Investigation of the effect of oscillation period on boiling two-phase flow behavior
The computational domain and thermal hydraulic conditions are the same as those for boiling two-phase flow in the parallel-plate channel, as described in the preceding section. The oscillation acceleration was applied at t = 0 s, after steady boiling flow was obtained.
Nine cases of oscillation acceleration, as shown in Table 1, were applied in order to investigate the influence of the oscillation period of the oscillation acceleration upon the boiling two-phase flow behavior. As shown in the preceding section, the influence of the horizontal oscillation acceleration upon boiling flow was greater than the influence of the vertical oscillation acceleration. Therefore, only the horizontal oscillation acceleration was investigated in these analyses. The minimum oscillation period of 0.005 s, as listed in Table  1, is equal to half of the minimum time interval of structural analysis in a reactor. The maximum oscillation period of 1.2 s is almost equal to the computable physical time of about 1 s. In all cases, magnitude of the oscillation acceleration was set to 400 Gal. Case G in Table 1 Figure 8 shows the standard deviation distribution of void fraction fluctuation. In cases where the oscillation period is less than 0.01 s, the influence of the oscillation acceleration is small because the magnitude of the void fraction fluctuation is very small compared to that in the cases where the oscillation period is greater than 0.02 s. When the oscillation period is greater than 0.02 s, although the magnitude of the void fraction fluctuation increases with elevation, it decreases near the top of the heated region.
In cases where the oscillation period is between 0.02 s and 0.30 s, the standard deviation distributions varied significantly with the variation in the oscillation period. In Case F, the magnitude of the void fraction fluctuation was highest locally. Therefore, the distribution of void fraction fluctuation was significantly dependent on the oscillation period in this range.
Case A Case B Case C Case D Case E Case F Case G Case H Case I

Fig. 8. Standard deviation distribution of void fraction distribution
On the other hand, in cases where the oscillation period was greater than 0.30 s, the standard deviation distributions hardly varied with the variation in the oscillation period. Therefore, the influence of the oscillation acceleration is small in this range.
From the information above, it can be confirmed that the boiling two-phase flow analysis, which is consistent with the time-series data of oscillation acceleration and has a time period greater than 0.01 s, can be performed. This is because oscillation acceleration with an oscillation period of less than 0.01 s has very little influence on the boiling two-phase flow. In addition, the time variations in the void fraction in cases where the oscillation period is greater than 0.30 s are close to quasi-steady variation. This means that the computable physical time of about 1 s is enough to evaluate the response of the boiling two-phase flow to the oscillation acceleration. Therefore, it can be confirmed that effective analysis can be performed by extracting an earthquake motion of about 1 s at any time during an earthquake.

Application to the boiling two-phase flow analysis in a simulated fuel assembly excited by oscillation acceleration
Boiling two-phase flow in a simulated fuel assembly excited by oscillation acceleration was performed by the improved ACE-3D in order to investigate how the three-dimensional behavior of boiling two-phase flow in a fuel assembly under oscillation conditions is evaluated by the improved ACE-3D.

Computational condition
In this analysis, a 7 × 7 fuel assembly in a current BWR core is simulated, as shown in Fig. 9. Fuel rod diameter is 10.8 mm; the narrowest gap between fuel rods is 4.4 mm, and the axial heat length is 3.66 m.
Four subchannels surrounded by nine fuel rods without channel boxes are adopted as the computational domain shown in Fig. 9; this is the smallest domain that can describe the www.intechopen.com three-dimensional behavior of boiling two-phase flow. This computational domain was determined to reflect the basic thermal-hydraulic characteristics in fuel assemblies under earthquake conditions.
In this domain, single-phase water flows in from the bottom of the channel with a mass velocity of 1673 kg/m 2 s and inlet temperature of 549.15 K. At the exit of the computational domain, pressure was fixed at 7.1 MPa. The mass velocity, inlet temperature, and exit pressure reflect the operating conditions in a current BWR core. The core thermal power is 351.9 W. The axial power distribution of the fuel-rod surfaces is shown in Fig. 9 and it simulates the power distribution in a current BWR core. Figure 10 shows the boundary conditions and the computational block divisions. Here, the non-slip condition is set for each fuel-rod surface, and the slip condition is set for each symmetric boundary. In this analysis, the computational domain was divided into 9 blocks. The computational grids in each block have 10 and 256 grids in the radial and axial directions, respectively. The number of grids in the peripheral direction is as follows: 30 in block 1, block 3, block 7, and block 9; 60 in block 2, block 4, block 6, and block 8; and 120 in block 5.
In this study, the boiling two-phase flow analysis was performed under steady-state conditions to obtain a steady boiling two-phase flow. Subsequently, oscillation acceleration was applied. The time when the oscillation acceleration was applied is regarded as t = 0 s.  In this analysis, in-phase sine wave acceleration was applied in the X and Y directions as shown by the black arrow in Fig. 10. The magnitude and oscillation period of the oscillation acceleration in the X and Y directions were 400 Gal and 0.2 s, respectively, as shown in Fig.  11; these values are based on actual earthquake data measured in the Kashiwazaki-Kariwa nuclear power plant. The computable physical time in this analysis after applying the oscillation acceleration was 1 s based on the results described in section 2.

Results and discussions
After applying the oscillation acceleration, the void fraction distribution fluctuated with the same period as the oscillation acceleration. Figure 12 shows the isosurfaces of the void fraction at t = 0.8 s. In the whole area where boiling occurs, the void fraction in the center of the subchannel was relatively low, and the void fraction concentrated in the positive directions along the X and Y axes was high.  Figure 13 shows the time variation in the void fraction at Z = 2.3 m in the upstream region of Fig. 12. The oscillation acceleration did not act at t = 0.8 s and t = 0.9 s and acted in the direction of the black arrow shown in Fig. 13(b). At t = 0.8 s, a high void fraction could be seen near the fuel-rod surface in the narrowest region between the fuel rods, as indicated by red circles in Fig. 13(a). At t = 0.85 s, a high void fraction moved in a direction opposite to the oscillation acceleration as shown in Fig. 13(b). At t = 0.9 s, a high void fraction could be seen near the fuel-rod surface in the regions marked by red circles in Fig. 13(c). This indicates that the magnitude of void fraction fluctuation at Z = 2.3 m is particularly large near the fuel-rod surface. This tendency of void fraction fluctuation is the same between t = 0.9 s and 1.0 s, when the oscillation acceleration acted in a direction opposite to that of the black arrow. Figure 14 shows the time variation in the void fraction at Z = 3.4 m in the downstream region of Fig. 12. The black arrow shows the direction in which the oscillation acceleration acts. At t = 0.78 s, the vapor phase moved in a direction opposite to the black arrow and concentrated in the regions marked by red circles, shown in Fig. 14(a). At t = 0.8 s, the void fraction in the region marked by the red circles in Fig. 14(b) increased. In addition, a high void fraction could also be seen away from the fuel rod surface, as shown by the blue circles in Fig. 14(b). The high void fraction in the regions marked by the blue circles in Fig. 14(b) split, and the high void fraction in the blue circles in Fig. 14(c) was formed at t = 0.82 s. The high void fraction regions represented by red and blue circles in Fig. 14(c) moved in a direction opposite to the black arrow as shown Fig. 14(d). While the void fraction regions indicated by the red and blue circles in Fig. 14(d) decreased as shown in Fig. 14(e), high void fraction was concentrated in the regions marked by red circles; high void fraction could also be seen in the regions away from the fuel rod surface, such as the regions indicated by the blue circles at t = 0.9 s, as shown in Fig. 14(f). Near the fuel rod surface, void fluctuation with a different period to that of the oscillation acceleration was seen while the magnitude of the void fraction was relatively small. The black arrow shows the direction in which the oscillation acceleration acts. The lift force in the red circles in Fig. 16(a) to Fig. 16(c) acted in a direction facing away from the fuel-rod surface. Hence, the vapor velocity was directed along the red arrow, as shown in Fig. 15; this is a direction opposite, but not parallel, to the black arrow. A high void fraction could be seen away from the fuel-rod surface, as shown by blue circles in Fig. 14(b). The lift force in the regions of the blue circles in Fig. 16(b) and Fig. 16(c) also acts in a direction facing away from the fuel rod surface. Hence, a high void fraction could be seen away from the fuel rod surface, as shown by red and blue circles in Fig. 14(c) and Fig. 14(d). Near the fuel rod surface, the magnitude and direction of the lift force were not uniform along the fuel rod and fluctuated with a period different from that of the oscillation acceleration. Consequently, void fraction fluctuation at Z = 3.4 m was significantly dependent on lift force fluctuation. Figure 17 shows the time variation in the Eotvos number at Z = 3.4 m and also shows a range of Eotvos number from 4 to 10 for which the effect of bubble deformation upon the lift force is dependent upon Eotvos number, as shown in Eq. (7). The black arrow shows the direction in which the oscillation acceleration acts. The red and blue circles in Fig. 17 correspond to regions where the magnitude of the lift force was large; the lift force acted in a direction facing away from the fuel rod surface, as shown in Fig. 16. In these regions, the effect of bubble deformation on the lift force was dominant because the Eotvos number exhibited high values. Near the fuel rod surface, the Eotvos numbers less than 4 and greater than 10 were mixed, indicating that the magnitude and direction of the lift force were not uniform near the fuel rod surface.  Figure 18 shows the variation in bubble diameter with time at Z = 3.4 m. The black arrow shows the direction in which the oscillation acceleration acts. Bubble diameters greater than 7 mm are distributed in the region where the Eotvos number is greater than 10, as shown in Fig. 17. The bubble diameter distribution shown in Fig. 18 is strongly inhomogeneous and physically invalid because large bubble diameters are mainly observed in small regions in the subchannel, while small bubble diameters of less than 3 mm are observed in the center of the subchannel. This strongly inhomogeneous bubble diameter distribution resulted in locally high Eotvos numbers and fluctuation in the direction of the lift force vectors.
The region where large bubble diameters are seen corresponds to the region of high void fraction, as shown in Fig. 14. According to Eq. (8), the bubble diameter is significantly dependent on the void fraction, and a local high void fraction results in a local large bubble diameter. Thus, a strongly inhomogeneous bubble diameter distribution results from void fraction fluctuation.
It is necessary to adequately evaluate the influence of the void fraction upon bubble diameter in order to avoid a strongly inhomogeneous bubble diameter distribution under oscillation conditions.
According to our results, void fraction fluctuation in the downstream region is significantly dependent on the lift force caused by a strongly inhomogeneous bubble diameter distribution.

Conclusion
A new external force term, which can simulate the oscillation acceleration, was added to the momentum conservation equations in order to apply the three-dimensional two-fluid model analysis code ACE-3D under earthquake conditions.
A boiling two-phase flow excited by applying vertical and horizontal oscillation acceleration was simulated in order to confirm that the simulation can be performed under oscillation conditions. It was confirmed that the void fraction fluctuation with the same period as that of the oscillation acceleration could be calculated in the case of both horizontal and vertical oscillation acceleration.
The influence of the oscillation period of the oscillation acceleration on the boiling twophase flow behavior in a fuel assembly was investigated in order to evaluate the highest frequency necessary for the improved method to be consistent with the time-series data of oscillation acceleration and the shortest period of oscillation acceleration for which the boiling two-phase flow shows quasi-steady time variation. It was confirmed that a boiling two-phase flow analysis consistent with the time-series data of oscillation acceleration and with a time interval greater than 0.01 s, can be performed. It was also shown that an effective analysis can be performed by extracting an earthquake motion of about 1 s at any time during the earthquake.
The three-dimensional behavior of boiling two-phase flow in a fuel assembly under oscillation conditions was evaluated using a simulated fuel assembly excited by oscillation acceleration. On the basis of this evaluation, it was confirmed that void fraction fluctuation www.intechopen.com in the downstream region is significantly dependent on the lift force caused by a strongly inhomogeneous bubble diameter distribution and that it is necessary to adequately evaluate the influence of void fraction on bubble diameter in order to avoid strongly inhomogeneous bubble diameter distribution under oscillation conditions.

Acknowledgment
The present study includes the result of "Research of simulation technology for estimation of quake-proof strength of nuclear power plant" conducted by the University of Tokyo as Core Research for Evolutional Science and Technology (CREST). This research was conducted using a supercomputer of the Japan Atomic Energy Agency.