**1. Introduction**

156 Nuclear Reactors

Reisch, F. (2009). High Pressure Boiling Water Reactor, HP-BWR. Royal Institute of

Ross, S. B., El-Genk, M. S., and Matthews, R. B. (1988). Thermal Conductivity Correlation for

Ross, S. B., El-Genk, M. S., and Matthews, R. B. (1990). Uranium nitride fuel swelling

Routbort, J. L. (1972). High-Temperature Deformation of Polycrystalline Uranium Carbide. J.

Routbort, J. L., and Singh, R. N. (1975). Carbide and Nitride Nuclear Fuels. J. Nuclear

Saltanov, E., Monichan, R., Tchernyavskaya, E., and Pioro, I. (2009). Steam-Reheat Option

Schlichting, K. W., Padture, N. P., Klemens, P. G. (2001). Thermal conductivity of dense and porous yttria-stabilized zirconia.J. of Materials Science 31, 3003-3010. Seltzer, M. S., Wright, T. R., and Moak, D. P. (1975). Creep Behavior of Uranium Carbide-

Shoup, R. D., and Grace, W. R., (1977). Process Variables in the Preparation of UN

Stellrecht, D. E., Farkas, M. S., and Moak, D. P. (1968). Compressive Creep of Uranium

Tokar, M., Nutt, A. W., and Leary, J. A., (1970). Mechanical Properties of Carbide and

U.S. DOE Nuclear Energy Research Advisory Committee (2002). A Technology Roadmap

Villamere, B., Allison, L., Grande, L., Mikhael, S., Rodriguez-Prado, A., and Pioro, I. (2009).

Wang, S., Yuan, L. Q., and Leung, L. K. (2010). Assessment of Supercritical Heat-Transfer

Wang, Z., Naterer, G. F., and Gabriel, K. S. (2010). Thermal Integration of SCWR Nuclear

Watts, M. J., and Chou, C-T. (1982), Mixed Convection Heat Transfer to Supercritical

Wheeler, M. J. (1965). Thermal Diffusivity at Incandescent Temperatures by a Modulated

Zahlan, H., Groeneveld, D., and Tavoularis, S., Mokry, S., Pioro, I. (2011). Assessment of

for Generation IV Nuclear Energy Systems. Retrieved July 12, 2010, from The

Thermal Aspects for Uranium Carbide and Uranium Dicarbide Fuels in Supercritical Water-cooled Nuclear Reactors. Proc. ICONE-17, Brussels, Belgium,

Correlations against AECL Database for Tubes. Proc. 2nd Canada-China Joint Workshop on Supercritical Water-Cooled Reactors (CCSC-2010), Toronto, Ontario,

and Thermochemical Hydrogen Plants. Proc. 2nd Canada-China Joint Workshop on Supercritical Water-Cooled Reactors (CCSC-2010), Toronto, Ontario, Canada:

Supercritical Heat Transfer Prediction Methods. Proc. of the 5th International Symposium on Supercritical Water-Cooled Reactors (ISSCWR-5), Vancouver,

for SCWRs. Proc. ICONE-17, July 12-16, Brussels, Belgium, Paper #76061, 10 pages.

Technology, Nuclear Power Safety, Stockholm, Sweden.

Uranium Nitride Fuel. J. Nuclear Materials 151, 313-317.

Based Alloys. J. American Ceramic Society 58, 138-142.

Special Metals.(n.d.). Inconel alloy 600. Retrieved December 2, 2008, from http://www.specialmetals.com/products/inconelalloy600.php.

Carbide. J. American Ceramic Society 51, No. 8, 455-458.

http://www.osti.gov/bridge/purl.cover.jsp?purl=/4100394-tE8dXk/.

http://www.ne.doe.gov/genIV/documents/gen\_iv\_roadmap.pdf.

Canada: Canadian Nuclear Society, April 25-28, 12 pages.

Pressure Water, Proc. 7th, IHTC, Munich, Germany, 495–500.

Electron Beam Technique. Brit. J. Appl. Phys., Vol. 16, 365-376.

Canadian Nuclear Society, April 25-28, 12 pages.

British Colombia, Canada, March 13-16.

Nitride Reactor Fuels. J. Reactor Technology.

Generation IV International Forum:

July 12-16, Paper #75990, 12 pages.

Microspheres. J. American Ceramic Society 60, No. 7-8, 332-335.

correlation. J. Nuclear Materials 170, 169-177.

Nuclear Materials 44, 24-30.

Materials 58, 78-114.

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.

Development of an Analytical Method on Water-Vapor Boiling Two-Phase

= −

*lift*

( ) ρ ρ

*E*

to that corresponding to a large bubble diameter.

Bubble diameter is evaluated by the following equations.

σ

ρ

increase in the void fraction.

In Eq. (7),

where σ σ

σ

*t bubble*

*c*

α ρ ( ) − × int *gl l M c U U rotU lift lift g g*

− − <sup>=</sup> + − <sup>0</sup>

*lift lift wake* = +<sup>0</sup> *cc c*

1 exp 0.242 Re 0.288 1 exp 0.242 Re

<sup>&</sup>lt; =− + ≤ ≤

− >

where

*l g b l g l b*

Reynolds number; and *Db*, bubble diameter. *clift*0 in Eq. (7) describes the effect of shear flow and is a positive value. *cwake* describes the effect of bubble deformation and significantly depends on Eotvos number of bubble diameter. The coefficient of lift force *clift* is estimated by the summation of *clift*0 and *cwake*. If the Eotvos number of bubble diameter is small, *clift* is positive. However, if the Eotvos number of bubble diameter increases, *clift* is negative. Therefore, the lift force corresponding to a small bubble diameter acts in a direction opposite

*g D U UD*

*wake t t*

− − = = 2

, Re

*D X D XD <sup>b</sup>* =− + ( ) 1 *slug bubble slug slug*

2 3 *X X XX slug s s s g* 3 2 , 4 0.25

*U U g*

<sup>2</sup> max , min ,30 *bubble slug <sup>b</sup>*

*We D DD*

=− = − ( )

= =

*g l lg l*

maximum bubble diameter, which is an input parameter. Bubble diameter, *Db*, is estimated by linear interpolation of small bubble diameter, *Dbubble*, and slug diameter, *Dslug*, with a coefficient, *Xslug*, which is dependent on the void fraction. Therefore, *Db* increases with an

α

<sup>−</sup> <sup>−</sup>

is the surface tension coefficient; *We*, a critical Weber number of 5.0; and *Dbmax*, the

( )

 σ

ρ ρ

*cE E*

The lift force is calculated as follows.

Flow Characteristics in BWR Fuel Assemblies Under Earthquake Condition 159

1991), lift force (Tomiyama, 1995), interface friction force, and bubble diameter (Lilies, 1988).

( ) ( )

0 ( 4) 0.096 0.384 (4 10) 0.576 ( 10)

*bubble*

*bubble*

(7)

(8)

*t*

*E*

*t*

*E*

ρ

is the surface tension coefficient; *Et*, Eotvos number; Re*bubble*, the bubble

 ν

#### **2. 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 water– vapor 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:

$$\frac{\partial}{\partial t}(\alpha\_{\mathcal{g}}\rho\_{\mathcal{g}}) + \frac{\partial}{\partial \alpha\_{\mathcal{j}}}(\alpha\_{\mathcal{g}}\rho\_{\mathcal{g}}\mathcal{U}\_{\mathcal{g},\mathcal{j}}) = \Gamma \tag{1}$$

$$\frac{\partial}{\partial t}(\alpha\_l \rho\_l) + \frac{\partial}{\partial x\_j}(\alpha\_l \rho\_l \iota I\_{l,j}) = -\Gamma \tag{2}$$

Momentum conservation for vapor and liquid phases:

$$\frac{\partial \mathcal{U}\_{\mathcal{g},i}}{\partial t} + \mathcal{U}\_{\mathcal{g},j} \frac{\partial \mathcal{U}\_{\mathcal{g},i}}{\partial \mathbf{x}\_{j}} = -\frac{1}{\rho\_{\mathcal{g}}} \frac{\partial P}{\partial \mathbf{x}\_{i}} - \frac{M\_{\mathcal{g},i}^{\text{int}}}{\alpha\_{\mathcal{g}} \rho\_{\mathcal{g}}} - \frac{\Gamma^{\star}}{\alpha\_{\mathcal{g}} \rho\_{\mathcal{g}}} (\mathcal{U}\_{\mathcal{g},i} - \mathcal{U}\_{l,i}) + \frac{1}{\alpha\_{\mathcal{g}} \rho\_{\mathcal{g}}} \frac{\partial \mathcal{T}\_{\mathcal{g},ij}}{\partial \mathbf{x}\_{j}} + g\_{i} \tag{3}$$

$$\frac{\partial \mathcal{U}\_{l,i}}{\partial t} + \mathcal{U}\_{l,j} \frac{\partial \mathcal{U}\_{l,i}}{\partial \mathbf{x}\_{j}} = -\frac{1}{\rho\_{l}} \frac{\partial P}{\partial \mathbf{x}\_{i}} - \frac{M\_{l,i}^{\text{int}}}{\alpha\_{l} \rho\_{l}} - \frac{\Gamma^{-}}{\alpha\_{l} \rho\_{l}} (\mathcal{U}\_{g,i} - \mathcal{U}\_{l,i}) + \frac{1}{\alpha\_{l} \rho\_{l}} \frac{\partial \tau\_{l,ij}}{\partial \mathbf{x}\_{j}} + \mathcal{g}\_{i} \tag{4}$$

Internal energy conservation for vapor and liquid phases:

$$\frac{\partial}{\partial t} \left( \alpha\_{\mathcal{g}} \rho\_{\mathcal{g}} e\_{\mathcal{g}} \right) + \frac{\partial \left( \alpha\_{\mathcal{g}} \rho\_{\mathcal{g}} e\_{\mathcal{g}} \mathcal{U}\_{\mathcal{g},i} \right)}{\partial \mathbf{x}\_{j}} = -P \left| \frac{\partial \alpha\_{\mathcal{g}}}{\partial t} + \frac{\partial \left( \alpha\_{\mathcal{g}} \mathcal{U}\_{\mathcal{g},i} \right)}{\partial \mathbf{x}\_{j}} \right| + q\_{\mathcal{g}}^{w} + q\_{\mathcal{g}}^{\mathrm{int}} + \Gamma \cdot h\_{\mathcal{g}}^{\mathrm{sat}} \tag{5}$$

$$\frac{\partial}{\partial t}(\alpha\_l \rho\_l e\_l) + \frac{\partial \left(\alpha\_l \rho\_l e\_l \mathbf{U}\_{l,j}\right)}{\partial \mathbf{x}\_j} = -P \left[\frac{\partial \alpha\_l}{\partial t} + \frac{\partial \left(\alpha\_l \mathbf{U}\_{l,j}\right)}{\partial \mathbf{x}\_j}\right] + q\_l^w + q\_l^{\text{int}} - \Gamma \cdot h\_l^{\text{sat}} \tag{6}$$

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 *qw* and *qint* indicate wall heat flux and interfacial heat flux, respectively, and *hsat* 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. (3) and (4), respectively, are shear stress tensors. The liquid shear stress tensor, τ *l,ij*, is a summation of shear induced turbulence obtained and bubble induced turbulence (Bertodano, 1994).

Summation of interface stress, *Mint*, of the vapor and liquid phases is equal to zero. The interface stress is determined by correlations between turbulent diffusion force (Lehey,

ACE-3D has been developed by the Japan Atomic Energy Agency (JAEA) to simulate water– vapor or water–air two-phase flows at subcritical pressures. The basic equations of ACE-3D

> αρ

 αρ

> αρ

 αρ

α

*g gg gg g j j e U U*

α

Here, *t* represents time; *x*, a spatial coordinate; *U*, velocity; *P*, pressure; *g*, acceleration due to

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 *qw* and *qint* indicate wall heat flux and interfacial heat flux,

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

Summation of interface stress, *Mint*, of the vapor and liquid phases is equal to zero. The interface stress is determined by correlations between turbulent diffusion force (Lehey,

respectively, and *hsat* represents saturation enthalpy. Summation of volume ratios

*e U U*

*l ll ll l j j*

+ =− + + + −Γ⋅

+ =− + + + +Γ⋅

+ =− − − − + + ∂∂∂ ∂ int , ,, , , , , *l i l i* 1 1 *l i l ij*

*UU M <sup>P</sup> <sup>U</sup> U U <sup>g</sup> txx <sup>x</sup>*

+ =− − − − + +

*UU M <sup>P</sup> <sup>U</sup> U U <sup>g</sup> t xx <sup>x</sup>*

*g j gi li i j g i gg gg gg j*

*l j gi li i j l i ll ll ll j*

> ( ) α

*<sup>e</sup> <sup>P</sup> qq h t x tx* (5)

 ( ) α

, , int *l l l lj l l lj w sat*

*<sup>e</sup> <sup>P</sup> qq h t x tx* (6)

, , int *g g g gj g g gj w sat*

+ =Γ

+ =−Γ

*U*

( )

( )

α ρ

α ρ

, density. Subscripts *g* and *l* indicate vapor and liquid

*l,ij* in Eqs. (3) and (4), respectively, are shear stress tensors.

*l,ij*, is a summation of shear induced turbulence obtained and

*U*

(1)

(2)

τ

τ

(3)

(4)

α*g* and α*l*

( ) ∂ ∂

∂ ∂ *g g g g gj*, *j*

( ) ∂ ∂

∂ ∂ *l l l l lj*, *j*

( ) αρ

( ) αρ

*t x*

 αρ

> αρ

∂ ∂ <sup>∂</sup> <sup>∂</sup>

∂ ∂ ∂∂

∂ ∂ <sup>∂</sup> <sup>∂</sup>

∂ ∂ ∂∂

ρ

τ

τ

∂ ∂ ∂ Γ<sup>−</sup> ∂

∂ ∂ ∂ Γ<sup>+</sup> ∂

∂∂∂ ∂ int ,, , , , , , *g i g i* 1 1 *g i g ij*

*t x*

**2. Development of an analytical method of boiling two-phase flow in a fuel** 

**assembly under earthquake conditions** 

Mass conservation for vapor and liquid phases:

Momentum conservation for vapor and liquid phases:

ρ

ρ

Internal energy conservation for vapor and liquid phases:

( ) ( ) α ρ

( ) ( ) α ρ

> τ*g,ij* and

bubble induced turbulence (Bertodano, 1994).

α ρ

> α ρ

gravity; *e*, internal energy; and

is zero if Γ is positive.

The liquid shear stress tensor,

**2.1 Overview of ACE-3D** 

are shown below.

1991), lift force (Tomiyama, 1995), interface friction force, and bubble diameter (Lilies, 1988). The lift force is calculated as follows.

$$\begin{aligned} \label{eq:1} \mathcal{M}\_{lift}^{\text{int}} &= -c\_{lift} \sigma\_{\mathcal{S}} \mathcal{P}\_{\mathcal{S}} \left( \overline{U}\_{\mathcal{S}} - \overline{U}\_{l} \right) \times rot \overline{U}\_{l} \\\\ \begin{aligned} c\_{lift} &= c\_{lift0} + c\_{unck} \\\\ c\_{lift0} &= 0.288 \frac{1 - \exp\left(-0.242 \operatorname{Re}\_{luibtle}\right)}{1 + \exp\left(-0.242 \operatorname{Re}\_{luibtle}\right)} \end{aligned} \end{aligned} \tag{7}$$
 
$$\mathcal{L}\_{unckc} = \begin{cases} 0 & \text{( $E\_t < 4$ )} \\ -0.096 E\_t + 0.384 & \text{( $4 \le E\_t \le 10$ )} \\ -0.576 & \text{( $E\_t > 10$ )} \end{cases} \tag{7}$$
 
$$\text{where}$$

where

$$E\_t = \frac{g\left(\rho\_l - \rho\_g\right)D\_b^2}{\sigma}, \quad \mathrm{Re}\_{bubble} = \frac{\rho\_l \left|\mathcal{U}\_g - \mathcal{U}\_l\right| D\_b}{\nu}.$$

In Eq. (7), σ is the surface tension coefficient; *Et*, Eotvos number; Re*bubble*, the bubble Reynolds number; and *Db*, bubble diameter. *clift*0 in Eq. (7) describes the effect of shear flow and is a positive value. *cwake* describes the effect of bubble deformation and significantly depends on Eotvos number of bubble diameter. The coefficient of lift force *clift* is estimated by the summation of *clift*0 and *cwake*. If the Eotvos number of bubble diameter is small, *clift* is positive. However, if the Eotvos number of bubble diameter increases, *clift* 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.

$$\begin{aligned} D\_b &= \left( 1 - X\_{slug} \right) D\_{bubble} + X\_{slug} D\_{slug} \\\\ X\_{slug} &= 3X\_s^2 - 2X\_s^3, \quad X\_s = 4 \left( \alpha\_\mathcal{g} - 0.25 \right) \\\\ D\_{bubble} &= \frac{\sigma \mathcal{W} \mathcal{e}}{\rho\_l \left| \mathcal{U}\_\mathcal{g} - \mathcal{U}\_l \right|^2}, \quad D\_{slug} = \min \left( D\_{b \text{max}}, 30 \sqrt{\frac{\sigma}{\mathcal{g} \left( \rho\_\mathcal{g} - \rho\_l \right)}} \right) \end{aligned} \tag{8}$$

where σ is the surface tension coefficient; *We*, a critical Weber number of 5.0; and *Dbmax*, the maximum bubble diameter, which is an input parameter. Bubble diameter, *Db*, is estimated by linear interpolation of small bubble diameter, *Dbubble*, and slug diameter, *Dslug*, with a coefficient, *Xslug*, which is dependent on the void fraction. Therefore, *Db* increases with an increase in the void fraction.

Development of an Analytical Method on Water-Vapor Boiling Two-Phase

ρ

ρ

 αρ

> αρ

∂ ∂ ∂ Γ<sup>−</sup> ∂

∂∂∂ ∂ int , ,, , , , , *l i l i* 1 1 *l i l ij*

∂ ∂ ∂ Γ<sup>+</sup> ∂

∂∂∂ ∂ int ,, , , , , , *g i g i* 1 1 *g i g ij*

conservation equations (Eqs. (3) and (4)).

by seismographic observation.

cause instability in simulation results.

oscillation acceleration.

Flow Characteristics in BWR Fuel Assemblies Under Earthquake Condition 161

term, *fi*, which simulates the acceleration of oscillation, was added to the momentum

 αρ

 αρ

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

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

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

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

+ =− − − − + + +

+ =− − − − + + +

*g j gi li i i j g i gg gg gg j UU M <sup>P</sup> <sup>U</sup> U U g f t xx <sup>x</sup>*

*l j gi li i i j l i ll ll ll j UU M <sup>P</sup> <sup>U</sup> U U g f txx <sup>x</sup>*

( )

( )

α ρ

α ρ τ

τ

(13)

(14)

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:

$$
\alpha\_g \frac{\partial \mathbf{k}\_g}{\partial t} + \alpha\_g \mathbf{U}\_{g,j} \frac{\partial \mathbf{k}\_g}{\partial \mathbf{x}\_j} = \frac{\partial}{\partial \mathbf{x}\_j} \left( \alpha\_g \frac{\nu\_g^t}{\sigma\_k} \frac{\partial \mathbf{k}\_g}{\partial \mathbf{x}\_j} \right) + \alpha\_g \left( \Phi\_g - \varepsilon\_g \right) \tag{9}
$$

$$
\alpha\_l \frac{\partial \mathbf{k}\_l}{\partial t} + \alpha\_l \mathbf{U}\_{l,j} \frac{\partial \mathbf{k}\_l}{\partial \mathbf{x}\_j} = \frac{\partial}{\partial \mathbf{x}\_j} \left( \alpha\_l \frac{\nu\_l^t}{\sigma\_k} \frac{\partial \mathbf{k}\_l}{\partial \mathbf{x}\_j} \right) + \alpha\_l (\Phi\_l - \mathbf{c}\_l) \tag{10}
$$

Turbulent energy dissipation rate conservation for vapor and liquid phases:

$$
\alpha\_g \frac{\partial \varepsilon\_g}{\partial t} + \alpha\_g \mathcal{U}\_{\mathcal{g},j} \frac{\partial \varepsilon\_g}{\partial \mathbf{x}\_j} = \frac{\partial}{\partial \mathbf{x}\_j} \left( \alpha\_g \frac{\nu\_g^t}{\sigma\_e} \frac{\partial \varepsilon\_g}{\partial \mathbf{x}\_j} \right) + \alpha\_g \left( \mathcal{C}\_{\varepsilon 1} \frac{\Phi\_g \varepsilon\_g}{k\_\mathcal{g}} - \mathcal{C}\_{\varepsilon 2} \frac{\varepsilon\_g^2}{k\_\mathcal{g}} \right) \tag{11}
$$

$$\alpha\_l \frac{\partial \varepsilon\_l}{\partial t} + \alpha\_l \mathcal{U}\_{l,j} \frac{\partial \varepsilon\_l}{\partial \mathbf{x}\_j} = \frac{\partial}{\partial \mathbf{x}\_j} \left( \alpha\_l \frac{\nu\_l^t}{\sigma\_\varepsilon} \frac{\partial \varepsilon\_l}{\partial \mathbf{x}\_j} \right) + \alpha\_l \left( \mathcal{C}\_{\varepsilon 1} \frac{\Phi\_l \varepsilon\_l}{k\_l} - \mathcal{C}\_{\varepsilon 2} \frac{\varepsilon\_l^2}{k\_l} \right) \tag{12}$$

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.

#### **2.2 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

In ACE-3D, the two-phase flow turbulent model based on the standard k-ε model

 α

 α

*l l lj l lll j j kj*

ε

ε

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

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.

**2.2 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

*t l l l l ll l*

σ

 νε

*t g g g g gg g*

∂∂ ∂ Φ ∂ += + − ∂ ∂∂ ∂

∂∂ ∂ Φ ∂ += + − ∂ ∂∂ ∂

*U C C*

, 1 2

*U C C*

σ

 νε

+ = + Φ−

*g g gj g ggg j j kj*

+ = + Φ−

∂∂ ∂ ∂

*kk k*

*g g g g*

∂ ∂∂ ∂

∂∂ ∂ ∂

*l l ll*

*kk k <sup>U</sup>*

Turbulent energy dissipation rate conservation for vapor and liquid phases:

*g g gj g g*

*l l lj l l*

∂ ∂∂ ∂ ,

> α

 α

,

*U*

 ε

> ε

ν

*t*

σ

ν

*t*

σ

, 1 2

( )

αε

( )

ε

 ε  ε

> ε

2

 ε  ε

2

 αε

*t xx x* (10)

ε

ε

*t xx x k k* (12)

 α

*t xx x k k* (11)

*j j j gg*

 α

*j j j ll*

*t xx x* (9)

(Bertodano, 1994) is introduced below.

ε

αα

ε

analysis of large-scale domains.

αα

Turbulent energy conservation for vapor and liquid phases:

αα

αα

term, *fi*, which simulates the acceleration of oscillation, was added to the momentum conservation equations (Eqs. (3) and (4)).

$$\frac{\partial \mathcal{U}\_{\mathbf{g},i}}{\partial t} + \mathcal{U}\_{\mathbf{g},j} \frac{\partial \mathcal{U}\_{\mathbf{g},i}}{\partial \mathbf{x}\_{j}} = -\frac{1}{\rho\_{\mathcal{g}}} \frac{\partial P}{\partial \mathbf{x}\_{i}} - \frac{M\_{\mathbf{g},i}^{\text{int}}}{a\_{\mathcal{g}} \rho\_{\mathcal{g}}} - \frac{\Gamma^{+}}{a\_{\mathcal{g}} \rho\_{\mathcal{g}}} \left(\mathcal{U}\_{\mathbf{g},i} - \mathcal{U}\_{l,i}\right) + \frac{1}{a\_{\mathcal{g}} \rho\_{\mathcal{g}}} \frac{\partial \mathsf{\mathcal{T}}\_{\mathbf{g},\vec{\eta}}}{\partial \mathbf{x}\_{j}} + g\_{i} + f\_{i} \tag{13}$$

$$\frac{\partial \mathcal{U}\_{l,i}}{\partial t} + \mathcal{U}\_{l,j} \frac{\partial \mathcal{U}\_{l,i}}{\partial \mathbf{x}\_{j}} = -\frac{1}{\rho\_{l}} \frac{\partial P}{\partial \mathbf{x}\_{i}} - \frac{M\_{l,i}^{\text{int}}}{\alpha\_{l} \rho\_{l}} - \frac{\Gamma^{-}}{\alpha\_{l} \rho\_{l}} \left(\mathcal{U}\_{\text{g.}} - \mathcal{U}\_{l,i}\right) + \frac{1}{\alpha\_{l} \rho\_{l}} \frac{\partial \sigma\_{l,ij}}{\partial \mathbf{x}\_{j}} + \mathcal{g}\_{i} + f\_{i} \tag{14}$$

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

Development of an Analytical Method on Water-Vapor Boiling Two-Phase

Fig. 2. Time variation in oscillation acceleration

Fig. 3. Void fraction distribution at *t* = 0 s

center of the channel.

Flow Characteristics in BWR Fuel Assemblies Under Earthquake Condition 163

Figure 3 shows distribution of the void fraction at *t* = 0 s. Much of the void fraction was distributed near the wall at Z = 2.0 m, because the effect of the lift force was dominant at this time, and the lift force acted toward the wall. On the other hand, much of the void fraction was distributed in the center of the channel at Z = 3.66 m, because the effect of bubble deformation on evaluation of lift force was dominant, and the lift force acted toward the

The case of horizontal oscillation acceleration is shown in Fig. 4, which shows the time variation in the void fraction at Z = 2.0 m and Z = 3.66 m. The void fraction fluctuated in the horizontal direction with the same period as the oscillation acceleration; however, it moved in the direction opposite to the oscillation acceleration at both Z = 2.0 m and Z = 3.66 m.

Figure 5 shows the time variation in the horizontal velocity of liquid and vapor at Z = 2.0 m, where a positive value of velocity corresponds to the positive direction along the X axis, and a negative value of velocity corresponds to a negative direction along the X axis. The liquid velocity fluctuated in the same direction as the oscillation acceleration, while the vapor velocity fluctuated in the opposite direction. These tendencies in liquid and vapor velocities at Z = 2.0 m can also be seen at Z = 3.66 m. If oscillation acceleration is applied in the horizontal direction, a horizontal pressure gradient arises in a direction opposite to that of

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.

#### **2.3 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/m2. 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.

Fig. 2. Time variation in oscillation acceleration

and horizontal oscillation to simulate an earthquake in order to confirm that the boiling two-

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

**2.3 Confirmation of stability of the boiling two-phase flow simulation under oscillation** 

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/m2. 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

First, an analysis was performed without applying oscillation acceleration. After a steady boiling flow was attained, oscillation acceleration was applied. The time when the

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.

phase flow simulation can be performed under oscillation conditions.

quasi-steady time variation.

in Eq. (8) is set to the channel width of 8.2 mm.

Fig. 1. Computational domain

oscillation acceleration was applied is regarded as *t =* 0 s.

**conditions** 

Figure 3 shows distribution of the void fraction at *t* = 0 s. Much of the void fraction was distributed near the wall at Z = 2.0 m, because the effect of the lift force was dominant at this time, and the lift force acted toward the wall. On the other hand, much of the void fraction was distributed in the center of the channel at Z = 3.66 m, because the effect of bubble deformation on evaluation of lift force was dominant, and the lift force acted toward the center of the channel.

Fig. 3. Void fraction distribution at *t* = 0 s

The case of horizontal oscillation acceleration is shown in Fig. 4, which shows the time variation in the void fraction at Z = 2.0 m and Z = 3.66 m. The void fraction fluctuated in the horizontal direction with the same period as the oscillation acceleration; however, it moved in the direction opposite to the oscillation acceleration at both Z = 2.0 m and Z = 3.66 m.

Figure 5 shows the time variation in the horizontal velocity of liquid and vapor at Z = 2.0 m, where a positive value of velocity corresponds to the positive direction along the X axis, and a negative value of velocity corresponds to a negative direction along the X axis. The liquid velocity fluctuated in the same direction as the oscillation acceleration, while the vapor velocity fluctuated in the opposite direction. These tendencies in liquid and vapor velocities at Z = 2.0 m can also be seen at Z = 3.66 m. If oscillation acceleration is applied in the horizontal direction, a horizontal pressure gradient arises in a direction opposite to that of

Development of an Analytical Method on Water-Vapor Boiling Two-Phase

vertical oscillation acceleration case at any vertical position.

oscillation acceleration.

Flow Characteristics in BWR Fuel Assemblies Under Earthquake Condition 165

acceleration. The vertical oscillation acceleration caused fluctuations in the pressure in the channel, causing expansion and contraction of the vapor phase. This explains why the void

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

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

(a) Z = 2.0 m (b) Z = 3.66 m

**2.4 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 is the same as the horizontal oscillation acceleration case shown in section 2.3.

Fig. 6. Time variation in void fraction in the vertical oscillation acceleration case

fraction fluctuated with the same period as that of the oscillation acceleration.

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.

Fig. 4. Time variation in void fraction in the horizontal oscillation acceleration case

Fig. 5. Time variation in liquid and vapor velocities in the horizontal oscillation acceleration case

The case of vertical oscillation acceleration is shown in Fig. 6, which also shows the time variation in the void fraction at Z = 2.0 m and Z = 3.66 m. The distribution of the void fraction at Z = 2.0m and Z = 3.66 m fluctuated with the same period as that of the oscillation

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

(a) Z=2.0 m (b) Z=3.66 m

(a) Liquid velocity (b) Vapor velocity Fig. 5. Time variation in liquid and vapor velocities in the horizontal oscillation acceleration

The case of vertical oscillation acceleration is shown in Fig. 6, which also shows the time variation in the void fraction at Z = 2.0 m and Z = 3.66 m. The distribution of the void fraction at Z = 2.0m and Z = 3.66 m fluctuated with the same period as that of the oscillation

Fig. 4. Time variation in void fraction in the horizontal oscillation acceleration case

oscillation acceleration.

case

acceleration. The vertical oscillation acceleration caused fluctuations in the pressure in the channel, causing expansion and contraction of the vapor phase. This explains why the void fraction fluctuated with the same period as 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.

Fig. 6. Time variation in void fraction in the vertical oscillation acceleration case
