- Email: [email protected]
A coupling model for the two-stage core calculation method with subchannel analysis for boiling water reactors
Annals of Nuclear Energy 102 (2017) 77–84
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
A coupling model for the two-stage core calculation method with subchannel analysis for boiling water reactors Takeshi Mitsuyasu a,b,⇑, Motoo Aoyama a, Akio Yamamoto b a b
Hitachi, Ltd., 7-1-1, Omika, Hitachi-shi, Ibaraki 319-1292, Japan Department of Materials, Physics and Energy Engineering, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan
a r t i c l e
i n f o
Article history: Received 14 June 2016 Received in revised form 21 November 2016 Accepted 27 November 2016
Keywords: BWR Control rod Subchannel analysis Void fraction
a b s t r a c t The two-stage core analysis method is widely used for BWR core analysis. The purpose of this study is to develop a core analysis model coupled with subchannel analysis within the two-stage calculation scheme using an assembly-based thermal-hydraulics calculation in the core analysis. The model changes the 2D lattice physics scheme, and couples with 3D subchannel analysis which evaluates the thermal-hydraulics characteristics within the coolant flow area divided as some subchannel regions. In order to couple with these two analyses, some BWR fuel assembly parameters are assumed and verified. The developed model is evaluated for the heterogeneous problem with and without a control rod. The present model is especially effective for the control rod inserted condition. The present model can incorporate the subchannel effect into the current two-stage core calculation method. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction The neutronics analysis for LWR cores is conventionally divided into two stages, the first one is lattice physics and the second one is core analysis. The lattice physics basically covers a LWR fuel assembly, where neutronics analysis is carried out without coupling to the thermal-hydraulics calculation. The geometry of an assembly is accurately modeled, and finer or continuous energy groups are adopted than those used for core analysis. Since core analysis needs cross sections that vary for the state parameters in the core nodes, the lattice physics tabulates homogeneous cross section constants as the results of parametric calculations. Core analysis covers a whole LWR core where the neutronics analysis is coupled with the thermal-hydraulics calculation. Especially in BWR core analysis, some of the assembly cross sections (constants) are prepared in terms of the void fraction because the nodal-wise void fractions depend on the thermal-hydraulics calculation in core analysis. Conventional lattice physics uses a uniform void fraction inside an assembly because the calculated cross section is functionalized by one value equivalent to the nodal average void fraction. As fuel design conditions become more complicated, e.g., use of a large water rod, however, it is known that the radial void distri-
⇑ Corresponding author at: Hitachi, Ltd., 7-1-1, Omika, Hitachi-shi, Ibaraki 3191292, Japan. E-mail address: [email protected] (T. Mitsuyasu). http://dx.doi.org/10.1016/j.anucene.2016.11.045 0306-4549/Ó 2016 Elsevier Ltd. All rights reserved.
bution in a BWR fuel assembly, which is called the subchannel void distribution, becomes heterogeneous and has an impact on the neutronics core characteristics. Katono et al. (2015) experimentally observed the subchannel void distribution using X-ray CT. A single BWR assembly study of neutronics coupled with a subchannel code revealed a difference in the multiplication factor due to the uniform and subchannel void distributions (Jatuff et al., 2006). Furthermore, a model coupled with neutronics and subchannel analysis was proposed for the single-assembly geometry (Ama et al., 2002). This method adopted an axially stacked fine mesh 2D neutronics analysis and assembly-based subchannel analysis. Ikehara et al. (2008) improved the method so that it could carry out the two-step neutronics analysis with subchannel analysis which consisted of the 2D fine mesh neutronics and axially stacked 1D problem with the nodal homogenized cross section obtained from the 2D neutronics analysis. On the other hand, full core analysis methods including subchannel effects have also been presented. Tada et al. (2011) directly combined the fine mesh SP3 method and subchannel analysis. In this paper, the authors extend their previous work (Mitsuyasu et al., 2014) to include more comprehensive results and they also demonstrate the heterogeneous problem with a control rod. In the present study, a calculation scheme is developed and evaluated which incorporates the subchannel analysis method into the two-stage BWR core analysis while maintaining the conventional neutronics and thermal-hydraulics calculation scheme. Utilization of the present core analysis is useful for back-fitting
78
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
to the conventional full core BWR analysis method. The conventional BWR core analysis in this study means the core analysis coupling with the assembly-based thermal-hydraulics calculation. In the paper, the calculation model is described in Section 2, numerical results are described in Section 3, and conclusions are given in Section 4.
to the one required by the lattice physics calculation. They are normalized to the average void fraction in the lattice physics because, in most cases, the nodal average void fraction obtained from subchannel analysis does not completely satisfy the average void fraction required by the lattice physics calculation. The normalized subchannel void fractions are given as:
2. Calculation model
v~ v~ s ¼ v s;n ~ L vn
Typically, lattice physics is performed for 2D radial cross sections of a fuel assembly for which the flow area has a fixed uniform average void fraction. On the other hand, subchannel analysis is performed for a 3D fuel assembly because the thermal-hydraulics model in the subchannel code needs to trace the changes of water and the void condition from the bottom to the top of the assembly. To be consistent with lattice physics and subchannel analysis, changes should be essentially calculated by a 3D fuel assembly. In order to incorporate the effect of the subchannel void fraction into the 2D lattice physics calculation, a new calculation model should be established. The proposed calculation scheme is illustrated in Fig. 1. The 2D lattice physics calculation normally has a fixed average void fraction over all the subchannel regions, for example, 0, 40 or 70%. In the present model, the 2D lattice physics calculation must be done with subchannel-wise void fractions. For subchannel analysis, the pin-by-pin radial power distribution can be obtained from the lattice physics. The axial power profiles however cannot be obtained from the lattice physics. Then, axial power profiles should be assumed. In the model, the axial cosine profile is used for the validation, but any axial power profile is allowed since the dependency of axial power profiles is negligible. The verification of the dependency on axial power profiles is described in Section 3.2. The model assumed the radial pin power distribution of the infinite lattice even though the radial pin power distribution is affected by the neighbor assemblies in an actual BWR core. The lowest flow and highest power of the assembly during operation is necessary as the thermal hydraulic condition in subchannel analysis because the average void fraction at the top of the assembly at the rated power and flow would not reach the void fraction which is necessary for the lattice physics calculation. The subchannel void fractions in subchannel analysis are obtained from the axial node for which the nodal average void fraction is closest
ð1Þ
~ s is a normalized subchannel void fraction at the region of where v subchannel number s, v s;n is the void fraction at the subchannel region number s and axial node n obtained from subchannel analy~ n is the average void fraction in the axial node n of the subsis, v channel analysis, and v~ L is the required average void fraction in lattice physics. Lattice physics and subchannel analysis are iterated until their results convergence. The convergence criteria are given as:
k1;iþ1 k1;i < 3rk k
ð2Þ
1;i
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uXN 2 u t j ðPiþ1;j Pi;j Þ N
< 3rp
ð3Þ
where k1 is k-infinity, i is the number of iteration and rk is the standard deviation of k-infinity of lattice physics using the Monte Carlo method. P j is fission rate of the jth fuel rod, N is number of fuel rods and rp is the standard deviation of pin-wise fission rate. Eq. (3) means the root mean square of the pin-wise fission rate difference between the current iteration and previous iteration. rk and rp are determined by a Monte Carlo calculation during the lattice physics calculation. Those iterations are carried out for all cases of voided conditions. In the depletion calculation, the subchannel analysis is coupled after the depletion calculation over all burnup points, instead of iteration at each point, in order to simplify the calculation scheme. In other words, subchannel void fraction is fixed during the depletion calculation. The core analysis is performed in the same manner as the conventional methods because the nuclear constants obtained by this model are tabulated just as those for the conventional methods are. 3. Numerical results and discussion 3.1. Calculation codes and conditions
Start Lattice physics Converged?
No
Pin power
Axial power T-H condition (Assumed)
Yes End
Subchannel analysis Select the node & Normalize the void distribution
Subchannel void distribution
Depletion calculation Burnup Branch calculation Point of calculation coupling with subchannel analysis Fig. 1. The calculation scheme of subchannel coupling analysis.
For evaluation of the developed model, VMONT (Morimoto et al., 1989) is used for the lattice physics calculation, SILFEED (Tomiyama et al., 1988) for the subchannel analysis and 3D-DRM (Hino et al., 2010) for the core analysis. JENDL-4.0 is used as the nuclear data (Shibata et al., 2011). The VMONT code is based on the Monte Carlo neutron transport method with a pseudoscattering method. This code uses a multi-group model for the neutron spectrum calculation, and the total number of energy groups is 190. The SILFEED code is a subchannel analysis code which consists of a subchannel analysis part and a film flow analysis part. To evaluate the developed model, the two parts are combined using results of the subchannel analysis to correct the vapor and droplet mass fluxes for the film flow analysis. SILFEED can predict the critical power within ±7% error (Tomiyama et al., 1988). 3D-DRM is introduced in order to minimize the error of neutronics characteristics and to clarify the effect of the developed model. The VMONT with multiple assemblies is used for the reference. The statistical uncertainty of the reference neutron infinite multiplication factor is 0.02%dk/k and that of the fuel rod-wise fission rate is 0.2%.
79
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
2.0
Relative power
In VMONT, the number of tracked neutrons is 1.0 107 at each state point. The statistical uncertainty of the neutron infinite multiplication factor is approximately 0.02%dk/k and that of the fuel rod-wise fission rate is approximately 0.2%. These values are set to rk and rp in Eqs. (2) and (3). As a result, convergence criteria of 0.06%dk/k and 0.6% are obtained for k-infinity and pin-wise fission rate, respectively. These conditions are the same as the reference with multiple assemblies. In SILFEED, the cross-sectional geometry is the same as in VMONT. The axial power profile and thermal-hydraulics condition, however, are assumed. Fuel active length is 3.7 m, mass flow is 45 t/h and assembly power is 7 MW. The number of axial nodes is 72. The subchannel domain decomposition of the 8 8 BWR fuel assembly with a large water rod is shown in Fig. 2. The subchannel region is defined as the gap between fuel rods. The surfaces of the fuel rods are treated as heated walls and the surfaces of the channel wall and the large water rod are treated as unheated walls. The fuel rod power is radially distributed according to fuel rod power obtained by VMONT. In 3D-DRM, the number of energy groups is 3: Group 1 is from 5.53 keV to 10.0 MeV, Group 2 is from 0.625 eV to 5.53 keV, and Group 3 is from 0.0 to 0.625 eV. The surface of the core node is divided into 4 sub-surfaces in each direction and 4 sub-angles. Fuel rods in a node are axially divided by 4. These conditions are the same as for the authors’ previous core analysis (Mitsuyasu et al., 2014). The sub-response matrices are calculated in VMONT according to the above surface, fuel rods and energy decomposition. This number of surfaces is indicated to be appropriate for core analysis (Ishii et al., 2009). For burnup analysis, the spectral history method using boundary current is used for this evaluation (Mitsuyasu et al., 2010).
1.5 1.0 cosine peaked lower peaked upper peaked
0.5 0.0 1
Bottom
11
21 31 41 51 61 Sub channel analysis axial node
71
Top
Fig. 3. Axial power profiles.
Channel box Fuel rod Water rod Water gap
Low enrichment fuel
Gd-containing fuel rod
3.2. Dependency of axial power profiles In order to show the dependency of axial power profiles in a subchannel analysis, three typical power profiles are used. Fig. 3 shows their profiles. The test 2D assemblies are the two types shown in Fig. 4. Each assembly consists of 8 8 fuel rods and a large water rod. The assemblies have two values of average enrichment. The average enrichment of low enrichment fuel is 1.3 wt%. The high enrichment fuel has 3.8 wt% enrichment and this assembly includes eight Gd-containing fuel rods. The boundary condition is mirror reflective for both radial and axial boundary. The average void fractions in the low and high enrichment fuel are 40%. The calculations for two assemblies are carried out using each of the three axial power profiles. Table 1 shows the results of each case. The kinfinity differences are less than 0.05%dk and root mean squares (RMSs) of pin-wise fission rate difference are less than 0.45%. All values are less than the convergence criteria. As a result of those comparisons, the effect of axial power profile is concluded to be negligible in the present analysis.
Channel box Fuel rod Water rod Subchannel region Boundary of subchannel Fig. 2. Example of subchannel domain decomposition.
High enrichment fuel Fig. 4. Test assemblies.
3.3. Dependency of the number of the iteration In order to show the dependency of the number of the iteration between the lattice physics and subchannel analysis calculations, the convergences of k-infinity and pin-wise fission rate are evaluated. The test assemblies and the average void fraction are the same as those in Section 3.2. Fig. 5 shows the convergence of kinfinity. The difference of k-infinity indicates the difference between each iteration and the 10th iteration. The k-infinity converged just after the 2nd iteration. In terms of pin-wise fission rate convergence, Fig. 6 shows the RMSs of pin-wise power difference between each iteration and the 10th iteration. The RMS is 0.4% after the 2nd iteration. The differences between the 2nd and 3rd iterations are below the convergence criteria which are shown in Eqs. (2) and (3), consequently, at least 3 iterations are necessary in the model. 3.4. Heterogeneous problem without a control rod For the comparison, some definitions are clarified. The ‘‘reference subchannel effect” is defined as the difference of reference results with and without subchannel coupling analysis. It is defined as:
kR ¼ kRðuniformÞ kRðsubchannelÞ
ð4Þ
80
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
Table 1 k-Infinity and pin-wise fission rate differences for three axial power profiles. k-infinity difference (%dk)
Pin-wise fission rate difference (RMS%)
Profile
High
Low
High
Low
Cosine peaked Lower peaked Upper peaked
Base 0.05 0.05
Base 0.05 0.01
Base 0.45 0.36
Base 0.30 0.22
For the assembly power, the ‘‘reference subchannel effect” and the ‘‘model difference” are defined as follows:
k-infinity difference (%dk)
0.40
AR ¼ ðARðuniformÞ ARðsubchannelÞ Þ=ARðsubchannelÞ
0.30 High enrichment fuel Low enrichment fuel
0.20
AD ¼ ðARðuniformÞ ARðsubchannelÞ AMðuniformÞ þ AMðsubchannelÞ Þ=ARðsubchannelÞ
0.10
-0.10 1
2
3
4 5 6 Number of iteration
7
8
ð8Þ
where AR is the reference subchannel effect on the assembly power evaluated by the reference method, AM is the model subchannel effect on the assembly power evaluated by the two-stage core calculation method, and AD is the model difference of the assembly power. The subscripts uniform and subchannel are the same as above. For the pin power, the differences are defined as follows:
0.00
9
PR ¼ ðPRðuniformÞ PRðsubchannelÞ Þ=P RðsubchannelÞ
Fig. 5. Convergence of k-infinity.
ð9Þ
PD ¼ ðPRðuniformÞ PRðsubchannelÞ PMðuniformÞ
RMS of pin-wise fission rate difference (%)
1.50
þ PMðsubchannelÞ Þ=PRðsubchannelÞ
0.50
0.00 1
2
3
4 5 6 Number of iteration
7
8
9
Fig. 6. Convergence of pin-wise fission rate.
where kRðuniformÞ is the result of k-infinity using the reference method without subchannel coupling analysis, kRðsubchannelÞ is the result of k-infinity using the reference method with subchannel coupling analysis. Then, kR shows the subchannel effect by the reference method. On the other hand, the ‘‘model subchannel effect” estimated by the present method is defined as:
kM ¼ kMðuniformÞ kMðsubchannelÞ
ð10Þ
where PR is the reference subchannel effect on the pin power evaluated by the reference method, PM is the model subchannel effect on the pin power evaluated by the two-stage core calculation method, and PD is the model difference of the pin power. The subscripts uniform and subchannel are the same as above. All results are compared using kR ; kD ; AR ; AD ; P R and P D . The model is appropriate to incorporate the subchannnel effect if each model difference (kD ; AD and P D ) is lower than the reference subchannel effect (kR ; AR and P R ). The model is evaluated in a heterogeneous problem with multiple enrichments without a control rod. The problem includes four assemblies in the 2D geometry which is shown in Fig. 7. Nodal average void fractions in each assembly are the same as those used in Section 3.2. The cosine power shape is used for the axial power profile in subchannel analysis. The calculation conditions about
High enrichment fuel Low enrichment fuel
1.00
ð5Þ
where kMðuniformÞ is the result of k-infinity using the two-stage core calculation method without subchannel coupling analysis. It is the conventional analysis using macroscopic cross section or subresponse matrices. kMðsubchannelÞ is the result of k-infinity using the two-stage core calculation method with subchannel coupling analysis. In general, kM differs from kR because the developed method includes some assumptions. If the difference between kR and kM , which is defined as ‘‘model difference”, is almost zero throughout the burnup, the effect of the assumptions adopted in the present study should be negligible. The ‘‘model difference” is defined as:
kD ¼ kR kM
ð7Þ
ð6Þ
Low enrichment fuel Channel box Fuel rod Water rod Water gap Gd-containing fuel rod
High enrichment fuel Boundary condition: mirror reflective in axial and radial boundaries Fig. 7. The 2D geometry of the heterogeneous problem without a control rod.
81
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
Table 2 Calculation conditions for burnup.
*
Reference (MultiVMONT)
Two-stage calculation method with the proposed method VMONT
3D-DRM
Direct calculation
- Burnup 0,40,70%* without the control rod - Branch 0,40,70%* with/without the control rod Power density 25,50,75 W/cc with/without the control rod (for Xenon calculation)
Spectral history method
Nodal average void fraction.
Difference of k-infinity (%dk)
0.50 0.40 0.30 0.20 0.10 0.00 -0.10 -0.20 -0.30
(kD) Model difference (D)
-0.40
Reference subchannel effect (k (R) R)
-0.50 0.0
5.0 10.0 Average exposure (GWd/t)
15.0
Fig. 8. Difference of k-infinity.
Bundle power difference (%)
2.00 1.50 1.00 0.50 0.00 -0.50
(AD) Model difference (D)
-1.00
Reference subchannel effect (A (R)R)
-1.50 -2.00 0.0
5.0 10.0 Exposure (GWd/t)
15.0
Fig. 9. Assembly power difference (low enrichment fuel).
2.00
Bundle power difference (%)
burnup are summarized in Table 2. The burnup cases in VMONT are 0, 40, 70% nodal average void fractions. The branching cases are 12 which are three void fractions and three power densities with/without the control rod. The k-infinity difference during burnup is shown in Fig. 8. The model difference is about 0.1%dk at the later stage of burnup, while the difference due to the subchannel effect is 0.25%dk at 0 GWd/t, the beginning of burnup. The reference subchannel effect changes from 0.3%dk at 0 GWd/t to 0.2%dk at 13 GWd/t. In the fourassembly geometry, estimated pin power differs from that in the single-assembly geometry with the reflective mirror boundary condition. Interference from a neighbor fuel assembly is largest at the beginning of burnup due to the ‘‘negative feedback” of burnup. Namely, if power of a fuel pin is higher at the beginning of burnup, its burnup becomes higher than that of other fuel pins, resulting in lower power at the later stage of burnup. In this context, it is natural that the reference and the model subchannel effect on k-infinity are largest at the beginning of the cycle and this effect generally decreases during burnup. The results of assembly power differences in low and high enrichment fuels are shown in Figs. 9 and 10, respectively. The reference subchannel effect is large at the beginning of burnup, and decreases with the burnup in both fuels. Both model differences are lower than the reference subchannel effect. The result is due to the Gd-containing rods in the high enrichment fuel. Fig. 11 shows the subchannel void fraction distribution in the high enrichment fuel obtained by the single assembly calculation at 0 GWd/t. Subchannel void fraction around the Gd-containing fuel rods is lower due to their lower power density. The lower void fraction promotes the neutron slowing down to produce thermal neutrons. Higher thermal neutron density around the Gd-containing fuel rods causes higher neutron absorption by Gd. As a result, the high enrichment fuel gives lower assembly power and thus the negative assembly power difference (Ama et al., 2002). Fig. 12 shows the pin-wise fission rate differences at 0 GWd/t. In the low enrichment fuel, the RMSs of pin-wise fission rate difference due to the reference subchannel effect and the model difference are 0.93% and 0.35%, respectively. In the high enrichment fuel, the respective RMSs are 1.36% and 0.93%. The model differences at both fuels are lower than the reference subchannel effect. In the low enrichment fuel, the higher reference subchannel effect is seen at the corner fuel rods, and the model differences at the corner fuel rods are lower than the reference subchannel effect. In the high enrichment fuel, the higher reference subchannel effect is seen at the fuel rods surrounding the water rod, and the model differences at the rods are lower than the reference subchannel effect. That is why the large subchannel effect is found along unheated surfaces such as the channel box walls and water rod surface because the subchannel void fractions decrease near the unheated surfaces. As a result, the developed model can decrease the reference subchannel effect, especially the difference of the pin power
1.50
Model difference (D) (AD)
1.00
(AR) Reference subchannel effect (R)
0.50 0.00 -0.50 -1.00 -1.50 -2.00 0.0
5.0 10.0 Exposure (GWd/t)
15.0
Fig. 10. Assembly power difference (high enrichment fuel).
37 42 42 40 38 38 40 41 35
42 48 46 44 40 40 45 47 41
42 46 42 38 37 33 39 45 40
40 44 38 33 33 30 33 40 39
38 40 37 33 33 37 40 38
38 40 33 30 33 33 38 44 40
40 45 39 33 37 38 42 46 42
41 47 45 40 40 44 46 48 42
35 41 40 39 38 40 42 42 37
Fig. 11. Subchannel void distribution of the high enrichment fuel (red values indicate the subchannels surrounding the Gd-containing rod).
82
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
0.83 -1.9 0.7
0.74 -1.2 0.7 1.10 -0.3 0.0
0.68 -0.6 0.4 1.27 0.2 0.4 1.19 0.7 0.4
1.09 -0.2 0.2 1.21 0.6 0.3 1.20 0.9 -0.3
Low
High
High
Low
1.07 0.2 0.4 1.18 0.8 0.2 1.18 0.9 0.0
1.09 -0.4 -0.2 1.19 0.4 0.2 1.13 0.7 -0.7 1.15 1.2 0.0 1.14 1.1 0.2 1.10 0.9 0.4
0.68 -1.3 -0.2 1.00 -0.3 -0.4 1.15 0.3 -0.4 1.11 0.4 -0.2 1.11 0.5 0.1 1.14 0.1 0.2 0.97 -0.5 0.0
0.73 -2.3 -0.1 0.65 -1.1 -0.2 1.00 -0.6 -0.4 0.95 -0.3 -0.3 0.96 -0.4 0.0 0.60 -0.8 0.1 0.64 -1.6 -0.6 0.71 -2.3 -0.3
Reference pin-wise fission rate Reference subchannel effect (PR) [%] Model difference (PD) [%]
1.26 1.3 1.2 1.36 1.0 0.7 1.33 1.1 1.4 1.24 0.6 0.9 1.28 1.1 2.0 1.38 1.4 0.6 1.42 1.5 0.5 1.34 1.0 0.3
1.28 0.9 -0.3 1.10 0.2 -0.3 0.98 -0.5 -0.2 0.24 0.0 -1.0 0.90 -0.4 1.0 1.08 0.3 0.6 0.00 0.5 -0.4
1.19 0.0 -0.5 0.94 -1.2 -0.3 0.24 -0.4 -0.4 0.88 -3.3 -0.2 0.94 -2.3 0.6 0.24 -0.4 -0.3
1.08 -0.6 -0.3 0.23 -0.4 -1.0 0.87 -3.1 0.3
1.10 -0.3 -1.0 0.82 -1.5 -0.6 0.91 -3.4 -1.8
1.18 0.3 -1.1 0.97 -0.4 -1.4 0.24 -0.4 0.1
1.21 0.9 -1.1 1.05 0.4 -1.2
1.13 1.1 -2.2
RMS difference [%] Low
High
Reference subchannel effect
0.93
1.36
Model difference
0.35
0.93
Fig. 12. Difference of pin-wise fission rate distribution.
Low enrichment fuel Channel box Fuel rod Water rod Water gap Gd-containing fuel rod
0.50 Difference of k-infinity (%dk)
Control rod
0.40 0.30 0.20 0.10 0.00 -0.10 -0.20
(kD) Model difference (D)
-0.30
(kR) Reference subchannel effect (R)
-0.40 -0.50 0.0
5.0 10.0 Average exposure (GWd/t)
15.0
Fig. 14. Difference of k-infinity at control rod inserted condition.
High enrichment fuel Boundary condition: mirror reflective in axial and radial boundaries Fig. 13. The 2D geometry of the heterogeneous problem with a control rod.
Table 3 Difference of k-infinity and assembly power of the heterogeneous problem with a control rod at 0 GWd/t. k-infinity (%dk) Case
Reference subchannel effect Model difference
0.37 0.21
Assembly power (%) Low with CR
High
1.61 0.39
0.88 0.46
Low without CR 1.13 0.94
at the fuel rods with large subchannel void differences, and the model has no negative influence. 3.5. Heterogeneous problem with a control rod Another evaluation of the proposed model is done with the control rod. The control rod is inserted during burnup. The calculated geometry is shown in Fig. 13. The average void fraction and axial power profile are the same as those used in Section 3.4. In the lattice physics calculation of the present evaluation, the burnup cal-
Bundle power difference (%)
2.00 1.50 1.00 0.50 0.00 -0.50 -1.00
(AD) Model difference (D)
-1.50
Reference subchannel effect (A (R)R)
-2.00 0.0
5.0 10.0 Exposure (GWd/t)
15.0
Fig. 15. Difference of assembly power (low enrichment fuel with a control rod at the north-west position).
culation without the control rod is performed in order to avoid a computational burden. The branch calculation is used to simulate the control rod inserted condition. These conditions are same as Table 2. Difference of k-infinity and assembly power at 0 GWd/t are summarized in Table 3. All differences of the reference subchannel effect are larger than the model difference at the beginning of burnup since the large power tilt is caused by the control rod and it makes lower void fraction around the control rod. The trend of kinfinity difference is shown in Fig. 14. The reference subchannel
83
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
2.00
1.50
(AD) Model difference (D)
1.00
Reference subchannel effect (R) (AR)
Bundle power difference (%)
Bundle power difference (%)
2.00
0.50 0.00 -0.50 -1.00 -1.50
1.50 1.00 0.50 0.00 -0.50 -1.00
(AD) Model difference (D)
-1.50
(AR) Reference subchannel effect (R)
-2.00
-2.00 0.0
5.0 10.0 Exposure (GWd/t)
0.0
15.0
effect on k-infinity is large at 0 GWd/t, and it decreases with burnup due to the negative feedback described in Section 3.4. But there is an undershoot until 0.12%dk at 12 GWd/t. The model difference for k-infinity is lower than the reference subchannel effect at 0 GWd/t and decreases by 5 GWd/t. This means that the model has good agreement with the reference analysis when burnup exceeds 5 GWd/t. 0.33 -12.9 0.4 0.61 -13.1 -0.8
0.35 -10.9 -1.4 0.83 -9.1 -0.9 0.93 -4.5 0.3
0.59 -7.6 -0.4 0.89 -5.5 -1.0 1.07 -2.0 0.6
Low High with CR High
0.65 -5.2 -1.1 0.98 -3.5 -0.8 1.16 0.0 -0.1
0.76 -2.8 0.4 1.09 -1.1 0.2 1.20 1.3 0.1 1.35 3.3 0.8 1.43 4.2 0.5 1.45 4.3 0.3
0.59 -2.7 0.2 1.03 -0.4 0.3 1.32 0.9 -0.7 1.38 3.0 0.2 1.45 4.2 0.3 1.56 4.4 0.3 1.37 4.3 0.7
Low w/o CR
15.0
Fig. 17. Difference of assembly power (low enrichment fuel at the south-east position).
Fig. 16. Difference of assembly power (high enrichment fuel).
0.32 -14.0 0.0
5.0 10.0 Exposure (GWd/t)
The respective changes with burnup of assembly power differences of low enrichment fuel with the control rod at the northwest position, high enrichment fuel, and low enrichment fuel at the south-east position are shown in Figs. 15–17. In low enrichment fuel with the control rod at the north-west position, the model difference is lower than the reference subchannel effect
0.80 -3.4 -1.3 0.76 -1.6 -1.3 1.24 0.5 -0.6 1.25 2.6 0.4 1.31 3.0 -0.2 0.84 3.4 0.7 0.93 3.1 0.9 1.06 2.2 0.9
1.09 -0.8 -1.0 1.19 -0.1 -0.8 1.17 0.6 0.4 1.12 1.7 1.1 1.19 2.1 1.6 1.31 2.8 0.6 1.38 3.1 0.8 1.33 2.9 1.2
1.16 -1.7 -2.8 1.00 -1.1 -1.7 0.90 -0.6 -0.6 0.23 0.0 -0.8 0.86 0.3 1.1 1.06 0.6 0.1 1.18 2.2 0.5 1.44 2.6 1.1
1.14 -2.9 -3.3 0.89 -2.8 -1.9 0.23 -0.9 -0.7 0.88 -3.3 -0.2 0.94 -2.7 0.3 0.25 -0.4 0.1 1.09 1.2 1.3 1.42 2.0 1.0
1.07 -3.1 -2.5 0.23 -0.9 -1.1 0.88 -3.7 0.1
1.11 -1.5 -1.8 0.84 -2.5 -1.3 0.93 -3.7 -1.8
0.97 -2.4 0.6 0.93 0.6 2.2 1.33 2.1 3.1
0.91 -2.2 1.2 0.25 0.4 -0.2 1.31 1.8 2.2
0.67 -2.5 -1.3
0.62 -2.4 -1.9 0.94 -1.2 -0.8
0.58 -1.2 -0.4 1.12 -0.4 -0.4 1.09 0.2 -0.4
0.94 -0.9 -0.6 1.10 -0.2 -0.6 1.13 0.5 -0.4
0.95 -1.0 -1.0 1.11 0.2 -0.3 1.15 0.6 -0.7
Reference pin-wise fission rate Reference subchannel effect (PR ) [%] Model difference (PD ) [%]
RMS difference [%]
Low with CR
High
Low w/o CR
Reference subchannel effect
5.35
1.84
0.99
Model difference
0.70
1.48
0.67
1.21 -1.1 -2.0 1.00 -1.5 -2.1 0.25 -0.4 0.7 0.94 -3.2 -1.2 0.90 -2.9 1.0 0.25 -0.4 0.2 1.03 0.5 1.1 1.41 1.7 2.4
1.24 -0.2 -1.9 1.08 -0.7 -1.9 1.00 -0.5 -1.0 0.85 -1.5 -0.3 0.24 0.0 0.0 0.98 -0.1 1.3 1.16 0.9 0.9 1.44 2.7 2.7
1.16 0.5 -2.4 1.25 0.2 -1.5 1.22 0.2 -0.7 1.15 -1.3 -1.6 1.13 -0.9 -0.2 1.25 -0.1 -0.2 1.35 1.7 0.9 1.34 2.2 2.4
0.99 -0.5 -0.3 1.14 1.0 0.3 1.13 0.8 -0.6 1.19 1.5 0.7 1.21 1.4 0.3 1.21 1.2 1.0
0.65 -0.9 0.0 1.00 0.2 0.2 1.20 0.4 0.3 1.20 0.5 -0.1 1.22 0.8 0.6 1.29 0.1 0.4 1.12 0.4 0.9
0.73 -1.8 0.6 0.68 -0.9 0.3 1.10 -0.4 -0.1 1.08 0.1 0.4 1.11 0.3 0.7 0.70 0.0 1.0 0.76 -1.2 0.8 0.85 -1.8 1.1
Fig. 18. Difference of pin-wise fission rate distribution of 0 GWd/t for control rod inserted condition.
84
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
during burnup, but it fluctuates. The pin power tilts due to the neighbor fuel assemblies change the subchannel void distribution from that in the infinite lattice. The model however calculates the lattice properties based on the subchannel void distribution made by pin powers in the infinite lattice. In this case, the model overestimates the assembly power in low enrichment fuel with the control rod, and causes undershoot later in the burnup due to the negative feedback. In high enrichment fuel, the model difference during burnup is lower than the reference subchannel effect. In low enrichment fuel at the south-east position, the model difference is almost the same as the reference subchannel effect except from 5 GWd/t to 12 GWd/t burnups. The duration of large difference in low enrichment at the south-east position agrees with the duration of the undershoot in low enrichment fuel with the control rod at the north-west position and the model difference is also estimated smaller than the difference of assembly power in high enrichment fuel. The difference is considered as a secondary effect because the total assembly power which is the summation of assembly powers is fixed and the underestimations should cause the overestimation. In the context of core analysis, the assembly with higher assembly power is important for evaluating the thermal margin such as the critical power ratio, and the present model can lower the difference in high enrichment fuel. Fig. 18 shows differences of the pin-wise fission rate at 0 GWd/t. In the low enrichment fuel with the control rod at the north-west position, the high enrichment fuel, and the low enrichment fuel at the south-east position, the respective RMSs of the reference subchannel effect are 5.35, 1.84 and 0.99%. And the respective RMSs of the model difference are 0.70, 1.48 and 0.67. The model can decrease the pin-wise fission rate difference at all assemblies, especially for the corner fuel rod neighboring the control rod. As a result, the model is effective for calculating the pin-wise fission rate with the control rod. 4. Conclusions The model of the two-stage core calculation method coupled with subchannel analysis for BWRs has been developed. The present model, which consists of the lattice physics analysis and the subchannel analysis parts, can be directly incorporated into the current core analysis model. In the present model, the lattice physics couples with the subchannel analysis which treats the 3D full assembly and uses the flow and power conditions to cover the maximum void fraction in the operation core. In order to simulate the power distribution inside a fuel assembly for the subchannel analysis, the cosine axial power profile was assumed. Through sensitivity calculations, the effect of the axial power profile was found to be negligible. The number of coupling iterations required for the model was three, which is sufficiently small. The present model
was validated in the heterogeneous problem with and without a control rod. The reference subchannel effect, which indicates the difference due to void distribution inside a fuel assembly, was always largest at the beginning of the burnup because the pin power tilt due to the neutronics heterogeneity was largest at the beginning of cycle and decreased as the burnup increased. The model difference, which is the difference between the reference subchannel effect and the model subchannel effect, was lower than the reference subchannel effect in terms of the k-infinity, assembly power, pin-wise fission rate and RMS. The present model was judged to be especially effective for the control rod inserted condition. The overall results indicated that the present model can incorporate the subchannnel effect into the current two-stage core calculation method. References Ama, T., Ikeda, H., Kosaka, S., Hyoudou, H., Takeda, T., 2002. Effect of moderator density distribution of annular flow on fuel assembly neutronic characteristics in boiling water reactor cores. J. Nucl. Sci. Technol. 39 (5), 487–498. Ama, T., Hyoudou, H., Takeda, T., 2002. Effect of radial void distribution within fuel assembly on assembly neutronic characteristics. J. Nucl. Sci. Technol. 39 (1), 90– 100. Hino, T., Ishii, K., Mitsuyasu, T., Aoyama, M., 2010. BWR core simulator using threedimensional direct response matrix and analysis of cold critical experiments. J. Nucl. Sci. Technol. 47 (5), 482–491. Ikehara, T., Kudou, Y., Tamitani, M., Yamamoto, M., 2008. Effect of Subchannel void fraction distribution on lattice physics parameters for boiling water reactor fuel bundles. J. Nucl. Sci. Technol. 45 (12), 1237–1251. Ishii, K., Hino, T., Mitsuyasu, T., Aoyama, M., 2009. Three-dimensional direct response matrix method using a Monte Carlo calculation. J. Nucl. Sci. Technol. 46 (3), 259–267. Jatuff, F., Giust, F., Krouthen, J., Helmersson, S., Chawla, R., 2006. Effects of void uncertainties on the void reactivity coefficient and pin power distributions for a 10 10 BWR assembly. Ann. Nucl. Energy 33 (2), 119–125. Katono, K., Nukaga, J., Fujimoto, K., Nagayoshi, T., Yasuda, K., 2015. Threedimensional time-averaged void fraction distribution measurement technique for BWR thermal hydraulic conditions using an X-ray CT system. J. Nucl. Sci. Technol. 52 (3), 388–395. Mitsuyasu, T., Ishii, K., Hino, T., Aoyama, M., 2010. Development of spectral history and pin power correction methods for pin-by-pin core analysis method using three-dimensional direct response matrix. J. Nucl. Sci. Technol. 47 (12), 1124– 1130. Mitsuyasu, T., Ishii, K., Nakadozono, N., Aoyama, M., A model of two-stage core calculation method coupled with subchannel analysis for boiling water reactors. In: Proc. PHYSOR-2014; 2014 Sep 28-Oct 4; Kyoto (Japan). [CD-ROM]. Morimoto, Y., Maruyama, H., Ishii, K., Aoyama, M., 1989. Neutronic analysis code for fuel assembly using a vectorized Monte Carlo method. Nucl. Sci. Eng. 103 (4), 351–358. Shibata, K., Iwamoto, O., Nakagawa, T., Iwamoto, N., Ichihara, A., Kunieda, S., Chiba, S., Furukawa, K., Otuka, N., Ohsawa, T., Murata, T., Matsunobu, H., Zukaran, A., Kameda, S., Katakura, J., 2011. JENDL-4.0: a new library for nuclear science and technology. J. Nucl. Sci. Technol. 48 (1), 1–30. Tada, K., Fujita, T., Endo, T., Yamamoto, A., Kosaka, S., Hirano, G., Nozaki, K., 2011. Application of quick subchannel analysis method for three-dimensional pin-bypin BWR core calculations. J. Nucl. Sci. Technol. 48 (12), 1437–1452. Tomiyama, A., Yokomizo, O., Yoshimoto, Y., Sugawara, S., 1988. Method of critical power prediction based on film flow model coupled with subchannel analysis. J. Nucl. Sci. Technol. 25 (12), 914–928.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
A coupling model for the two-stage core calculation method with subchannel analysis for boiling water reactors Takeshi Mitsuyasu a,b,⇑, Motoo Aoyama a, Akio Yamamoto b a b
Hitachi, Ltd., 7-1-1, Omika, Hitachi-shi, Ibaraki 319-1292, Japan Department of Materials, Physics and Energy Engineering, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan
a r t i c l e
i n f o
Article history: Received 14 June 2016 Received in revised form 21 November 2016 Accepted 27 November 2016
Keywords: BWR Control rod Subchannel analysis Void fraction
a b s t r a c t The two-stage core analysis method is widely used for BWR core analysis. The purpose of this study is to develop a core analysis model coupled with subchannel analysis within the two-stage calculation scheme using an assembly-based thermal-hydraulics calculation in the core analysis. The model changes the 2D lattice physics scheme, and couples with 3D subchannel analysis which evaluates the thermal-hydraulics characteristics within the coolant flow area divided as some subchannel regions. In order to couple with these two analyses, some BWR fuel assembly parameters are assumed and verified. The developed model is evaluated for the heterogeneous problem with and without a control rod. The present model is especially effective for the control rod inserted condition. The present model can incorporate the subchannel effect into the current two-stage core calculation method. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction The neutronics analysis for LWR cores is conventionally divided into two stages, the first one is lattice physics and the second one is core analysis. The lattice physics basically covers a LWR fuel assembly, where neutronics analysis is carried out without coupling to the thermal-hydraulics calculation. The geometry of an assembly is accurately modeled, and finer or continuous energy groups are adopted than those used for core analysis. Since core analysis needs cross sections that vary for the state parameters in the core nodes, the lattice physics tabulates homogeneous cross section constants as the results of parametric calculations. Core analysis covers a whole LWR core where the neutronics analysis is coupled with the thermal-hydraulics calculation. Especially in BWR core analysis, some of the assembly cross sections (constants) are prepared in terms of the void fraction because the nodal-wise void fractions depend on the thermal-hydraulics calculation in core analysis. Conventional lattice physics uses a uniform void fraction inside an assembly because the calculated cross section is functionalized by one value equivalent to the nodal average void fraction. As fuel design conditions become more complicated, e.g., use of a large water rod, however, it is known that the radial void distri-
⇑ Corresponding author at: Hitachi, Ltd., 7-1-1, Omika, Hitachi-shi, Ibaraki 3191292, Japan. E-mail address: [email protected] (T. Mitsuyasu). http://dx.doi.org/10.1016/j.anucene.2016.11.045 0306-4549/Ó 2016 Elsevier Ltd. All rights reserved.
bution in a BWR fuel assembly, which is called the subchannel void distribution, becomes heterogeneous and has an impact on the neutronics core characteristics. Katono et al. (2015) experimentally observed the subchannel void distribution using X-ray CT. A single BWR assembly study of neutronics coupled with a subchannel code revealed a difference in the multiplication factor due to the uniform and subchannel void distributions (Jatuff et al., 2006). Furthermore, a model coupled with neutronics and subchannel analysis was proposed for the single-assembly geometry (Ama et al., 2002). This method adopted an axially stacked fine mesh 2D neutronics analysis and assembly-based subchannel analysis. Ikehara et al. (2008) improved the method so that it could carry out the two-step neutronics analysis with subchannel analysis which consisted of the 2D fine mesh neutronics and axially stacked 1D problem with the nodal homogenized cross section obtained from the 2D neutronics analysis. On the other hand, full core analysis methods including subchannel effects have also been presented. Tada et al. (2011) directly combined the fine mesh SP3 method and subchannel analysis. In this paper, the authors extend their previous work (Mitsuyasu et al., 2014) to include more comprehensive results and they also demonstrate the heterogeneous problem with a control rod. In the present study, a calculation scheme is developed and evaluated which incorporates the subchannel analysis method into the two-stage BWR core analysis while maintaining the conventional neutronics and thermal-hydraulics calculation scheme. Utilization of the present core analysis is useful for back-fitting
78
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
to the conventional full core BWR analysis method. The conventional BWR core analysis in this study means the core analysis coupling with the assembly-based thermal-hydraulics calculation. In the paper, the calculation model is described in Section 2, numerical results are described in Section 3, and conclusions are given in Section 4.
to the one required by the lattice physics calculation. They are normalized to the average void fraction in the lattice physics because, in most cases, the nodal average void fraction obtained from subchannel analysis does not completely satisfy the average void fraction required by the lattice physics calculation. The normalized subchannel void fractions are given as:
2. Calculation model
v~ v~ s ¼ v s;n ~ L vn
Typically, lattice physics is performed for 2D radial cross sections of a fuel assembly for which the flow area has a fixed uniform average void fraction. On the other hand, subchannel analysis is performed for a 3D fuel assembly because the thermal-hydraulics model in the subchannel code needs to trace the changes of water and the void condition from the bottom to the top of the assembly. To be consistent with lattice physics and subchannel analysis, changes should be essentially calculated by a 3D fuel assembly. In order to incorporate the effect of the subchannel void fraction into the 2D lattice physics calculation, a new calculation model should be established. The proposed calculation scheme is illustrated in Fig. 1. The 2D lattice physics calculation normally has a fixed average void fraction over all the subchannel regions, for example, 0, 40 or 70%. In the present model, the 2D lattice physics calculation must be done with subchannel-wise void fractions. For subchannel analysis, the pin-by-pin radial power distribution can be obtained from the lattice physics. The axial power profiles however cannot be obtained from the lattice physics. Then, axial power profiles should be assumed. In the model, the axial cosine profile is used for the validation, but any axial power profile is allowed since the dependency of axial power profiles is negligible. The verification of the dependency on axial power profiles is described in Section 3.2. The model assumed the radial pin power distribution of the infinite lattice even though the radial pin power distribution is affected by the neighbor assemblies in an actual BWR core. The lowest flow and highest power of the assembly during operation is necessary as the thermal hydraulic condition in subchannel analysis because the average void fraction at the top of the assembly at the rated power and flow would not reach the void fraction which is necessary for the lattice physics calculation. The subchannel void fractions in subchannel analysis are obtained from the axial node for which the nodal average void fraction is closest
ð1Þ
~ s is a normalized subchannel void fraction at the region of where v subchannel number s, v s;n is the void fraction at the subchannel region number s and axial node n obtained from subchannel analy~ n is the average void fraction in the axial node n of the subsis, v channel analysis, and v~ L is the required average void fraction in lattice physics. Lattice physics and subchannel analysis are iterated until their results convergence. The convergence criteria are given as:
k1;iþ1 k1;i < 3rk k
ð2Þ
1;i
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uXN 2 u t j ðPiþ1;j Pi;j Þ N
< 3rp
ð3Þ
where k1 is k-infinity, i is the number of iteration and rk is the standard deviation of k-infinity of lattice physics using the Monte Carlo method. P j is fission rate of the jth fuel rod, N is number of fuel rods and rp is the standard deviation of pin-wise fission rate. Eq. (3) means the root mean square of the pin-wise fission rate difference between the current iteration and previous iteration. rk and rp are determined by a Monte Carlo calculation during the lattice physics calculation. Those iterations are carried out for all cases of voided conditions. In the depletion calculation, the subchannel analysis is coupled after the depletion calculation over all burnup points, instead of iteration at each point, in order to simplify the calculation scheme. In other words, subchannel void fraction is fixed during the depletion calculation. The core analysis is performed in the same manner as the conventional methods because the nuclear constants obtained by this model are tabulated just as those for the conventional methods are. 3. Numerical results and discussion 3.1. Calculation codes and conditions
Start Lattice physics Converged?
No
Pin power
Axial power T-H condition (Assumed)
Yes End
Subchannel analysis Select the node & Normalize the void distribution
Subchannel void distribution
Depletion calculation Burnup Branch calculation Point of calculation coupling with subchannel analysis Fig. 1. The calculation scheme of subchannel coupling analysis.
For evaluation of the developed model, VMONT (Morimoto et al., 1989) is used for the lattice physics calculation, SILFEED (Tomiyama et al., 1988) for the subchannel analysis and 3D-DRM (Hino et al., 2010) for the core analysis. JENDL-4.0 is used as the nuclear data (Shibata et al., 2011). The VMONT code is based on the Monte Carlo neutron transport method with a pseudoscattering method. This code uses a multi-group model for the neutron spectrum calculation, and the total number of energy groups is 190. The SILFEED code is a subchannel analysis code which consists of a subchannel analysis part and a film flow analysis part. To evaluate the developed model, the two parts are combined using results of the subchannel analysis to correct the vapor and droplet mass fluxes for the film flow analysis. SILFEED can predict the critical power within ±7% error (Tomiyama et al., 1988). 3D-DRM is introduced in order to minimize the error of neutronics characteristics and to clarify the effect of the developed model. The VMONT with multiple assemblies is used for the reference. The statistical uncertainty of the reference neutron infinite multiplication factor is 0.02%dk/k and that of the fuel rod-wise fission rate is 0.2%.
79
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
2.0
Relative power
In VMONT, the number of tracked neutrons is 1.0 107 at each state point. The statistical uncertainty of the neutron infinite multiplication factor is approximately 0.02%dk/k and that of the fuel rod-wise fission rate is approximately 0.2%. These values are set to rk and rp in Eqs. (2) and (3). As a result, convergence criteria of 0.06%dk/k and 0.6% are obtained for k-infinity and pin-wise fission rate, respectively. These conditions are the same as the reference with multiple assemblies. In SILFEED, the cross-sectional geometry is the same as in VMONT. The axial power profile and thermal-hydraulics condition, however, are assumed. Fuel active length is 3.7 m, mass flow is 45 t/h and assembly power is 7 MW. The number of axial nodes is 72. The subchannel domain decomposition of the 8 8 BWR fuel assembly with a large water rod is shown in Fig. 2. The subchannel region is defined as the gap between fuel rods. The surfaces of the fuel rods are treated as heated walls and the surfaces of the channel wall and the large water rod are treated as unheated walls. The fuel rod power is radially distributed according to fuel rod power obtained by VMONT. In 3D-DRM, the number of energy groups is 3: Group 1 is from 5.53 keV to 10.0 MeV, Group 2 is from 0.625 eV to 5.53 keV, and Group 3 is from 0.0 to 0.625 eV. The surface of the core node is divided into 4 sub-surfaces in each direction and 4 sub-angles. Fuel rods in a node are axially divided by 4. These conditions are the same as for the authors’ previous core analysis (Mitsuyasu et al., 2014). The sub-response matrices are calculated in VMONT according to the above surface, fuel rods and energy decomposition. This number of surfaces is indicated to be appropriate for core analysis (Ishii et al., 2009). For burnup analysis, the spectral history method using boundary current is used for this evaluation (Mitsuyasu et al., 2010).
1.5 1.0 cosine peaked lower peaked upper peaked
0.5 0.0 1
Bottom
11
21 31 41 51 61 Sub channel analysis axial node
71
Top
Fig. 3. Axial power profiles.
Channel box Fuel rod Water rod Water gap
Low enrichment fuel
Gd-containing fuel rod
3.2. Dependency of axial power profiles In order to show the dependency of axial power profiles in a subchannel analysis, three typical power profiles are used. Fig. 3 shows their profiles. The test 2D assemblies are the two types shown in Fig. 4. Each assembly consists of 8 8 fuel rods and a large water rod. The assemblies have two values of average enrichment. The average enrichment of low enrichment fuel is 1.3 wt%. The high enrichment fuel has 3.8 wt% enrichment and this assembly includes eight Gd-containing fuel rods. The boundary condition is mirror reflective for both radial and axial boundary. The average void fractions in the low and high enrichment fuel are 40%. The calculations for two assemblies are carried out using each of the three axial power profiles. Table 1 shows the results of each case. The kinfinity differences are less than 0.05%dk and root mean squares (RMSs) of pin-wise fission rate difference are less than 0.45%. All values are less than the convergence criteria. As a result of those comparisons, the effect of axial power profile is concluded to be negligible in the present analysis.
Channel box Fuel rod Water rod Subchannel region Boundary of subchannel Fig. 2. Example of subchannel domain decomposition.
High enrichment fuel Fig. 4. Test assemblies.
3.3. Dependency of the number of the iteration In order to show the dependency of the number of the iteration between the lattice physics and subchannel analysis calculations, the convergences of k-infinity and pin-wise fission rate are evaluated. The test assemblies and the average void fraction are the same as those in Section 3.2. Fig. 5 shows the convergence of kinfinity. The difference of k-infinity indicates the difference between each iteration and the 10th iteration. The k-infinity converged just after the 2nd iteration. In terms of pin-wise fission rate convergence, Fig. 6 shows the RMSs of pin-wise power difference between each iteration and the 10th iteration. The RMS is 0.4% after the 2nd iteration. The differences between the 2nd and 3rd iterations are below the convergence criteria which are shown in Eqs. (2) and (3), consequently, at least 3 iterations are necessary in the model. 3.4. Heterogeneous problem without a control rod For the comparison, some definitions are clarified. The ‘‘reference subchannel effect” is defined as the difference of reference results with and without subchannel coupling analysis. It is defined as:
kR ¼ kRðuniformÞ kRðsubchannelÞ
ð4Þ
80
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
Table 1 k-Infinity and pin-wise fission rate differences for three axial power profiles. k-infinity difference (%dk)
Pin-wise fission rate difference (RMS%)
Profile
High
Low
High
Low
Cosine peaked Lower peaked Upper peaked
Base 0.05 0.05
Base 0.05 0.01
Base 0.45 0.36
Base 0.30 0.22
For the assembly power, the ‘‘reference subchannel effect” and the ‘‘model difference” are defined as follows:
k-infinity difference (%dk)
0.40
AR ¼ ðARðuniformÞ ARðsubchannelÞ Þ=ARðsubchannelÞ
0.30 High enrichment fuel Low enrichment fuel
0.20
AD ¼ ðARðuniformÞ ARðsubchannelÞ AMðuniformÞ þ AMðsubchannelÞ Þ=ARðsubchannelÞ
0.10
-0.10 1
2
3
4 5 6 Number of iteration
7
8
ð8Þ
where AR is the reference subchannel effect on the assembly power evaluated by the reference method, AM is the model subchannel effect on the assembly power evaluated by the two-stage core calculation method, and AD is the model difference of the assembly power. The subscripts uniform and subchannel are the same as above. For the pin power, the differences are defined as follows:
0.00
9
PR ¼ ðPRðuniformÞ PRðsubchannelÞ Þ=P RðsubchannelÞ
Fig. 5. Convergence of k-infinity.
ð9Þ
PD ¼ ðPRðuniformÞ PRðsubchannelÞ PMðuniformÞ
RMS of pin-wise fission rate difference (%)
1.50
þ PMðsubchannelÞ Þ=PRðsubchannelÞ
0.50
0.00 1
2
3
4 5 6 Number of iteration
7
8
9
Fig. 6. Convergence of pin-wise fission rate.
where kRðuniformÞ is the result of k-infinity using the reference method without subchannel coupling analysis, kRðsubchannelÞ is the result of k-infinity using the reference method with subchannel coupling analysis. Then, kR shows the subchannel effect by the reference method. On the other hand, the ‘‘model subchannel effect” estimated by the present method is defined as:
kM ¼ kMðuniformÞ kMðsubchannelÞ
ð10Þ
where PR is the reference subchannel effect on the pin power evaluated by the reference method, PM is the model subchannel effect on the pin power evaluated by the two-stage core calculation method, and PD is the model difference of the pin power. The subscripts uniform and subchannel are the same as above. All results are compared using kR ; kD ; AR ; AD ; P R and P D . The model is appropriate to incorporate the subchannnel effect if each model difference (kD ; AD and P D ) is lower than the reference subchannel effect (kR ; AR and P R ). The model is evaluated in a heterogeneous problem with multiple enrichments without a control rod. The problem includes four assemblies in the 2D geometry which is shown in Fig. 7. Nodal average void fractions in each assembly are the same as those used in Section 3.2. The cosine power shape is used for the axial power profile in subchannel analysis. The calculation conditions about
High enrichment fuel Low enrichment fuel
1.00
ð5Þ
where kMðuniformÞ is the result of k-infinity using the two-stage core calculation method without subchannel coupling analysis. It is the conventional analysis using macroscopic cross section or subresponse matrices. kMðsubchannelÞ is the result of k-infinity using the two-stage core calculation method with subchannel coupling analysis. In general, kM differs from kR because the developed method includes some assumptions. If the difference between kR and kM , which is defined as ‘‘model difference”, is almost zero throughout the burnup, the effect of the assumptions adopted in the present study should be negligible. The ‘‘model difference” is defined as:
kD ¼ kR kM
ð7Þ
ð6Þ
Low enrichment fuel Channel box Fuel rod Water rod Water gap Gd-containing fuel rod
High enrichment fuel Boundary condition: mirror reflective in axial and radial boundaries Fig. 7. The 2D geometry of the heterogeneous problem without a control rod.
81
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
Table 2 Calculation conditions for burnup.
*
Reference (MultiVMONT)
Two-stage calculation method with the proposed method VMONT
3D-DRM
Direct calculation
- Burnup 0,40,70%* without the control rod - Branch 0,40,70%* with/without the control rod Power density 25,50,75 W/cc with/without the control rod (for Xenon calculation)
Spectral history method
Nodal average void fraction.
Difference of k-infinity (%dk)
0.50 0.40 0.30 0.20 0.10 0.00 -0.10 -0.20 -0.30
(kD) Model difference (D)
-0.40
Reference subchannel effect (k (R) R)
-0.50 0.0
5.0 10.0 Average exposure (GWd/t)
15.0
Fig. 8. Difference of k-infinity.
Bundle power difference (%)
2.00 1.50 1.00 0.50 0.00 -0.50
(AD) Model difference (D)
-1.00
Reference subchannel effect (A (R)R)
-1.50 -2.00 0.0
5.0 10.0 Exposure (GWd/t)
15.0
Fig. 9. Assembly power difference (low enrichment fuel).
2.00
Bundle power difference (%)
burnup are summarized in Table 2. The burnup cases in VMONT are 0, 40, 70% nodal average void fractions. The branching cases are 12 which are three void fractions and three power densities with/without the control rod. The k-infinity difference during burnup is shown in Fig. 8. The model difference is about 0.1%dk at the later stage of burnup, while the difference due to the subchannel effect is 0.25%dk at 0 GWd/t, the beginning of burnup. The reference subchannel effect changes from 0.3%dk at 0 GWd/t to 0.2%dk at 13 GWd/t. In the fourassembly geometry, estimated pin power differs from that in the single-assembly geometry with the reflective mirror boundary condition. Interference from a neighbor fuel assembly is largest at the beginning of burnup due to the ‘‘negative feedback” of burnup. Namely, if power of a fuel pin is higher at the beginning of burnup, its burnup becomes higher than that of other fuel pins, resulting in lower power at the later stage of burnup. In this context, it is natural that the reference and the model subchannel effect on k-infinity are largest at the beginning of the cycle and this effect generally decreases during burnup. The results of assembly power differences in low and high enrichment fuels are shown in Figs. 9 and 10, respectively. The reference subchannel effect is large at the beginning of burnup, and decreases with the burnup in both fuels. Both model differences are lower than the reference subchannel effect. The result is due to the Gd-containing rods in the high enrichment fuel. Fig. 11 shows the subchannel void fraction distribution in the high enrichment fuel obtained by the single assembly calculation at 0 GWd/t. Subchannel void fraction around the Gd-containing fuel rods is lower due to their lower power density. The lower void fraction promotes the neutron slowing down to produce thermal neutrons. Higher thermal neutron density around the Gd-containing fuel rods causes higher neutron absorption by Gd. As a result, the high enrichment fuel gives lower assembly power and thus the negative assembly power difference (Ama et al., 2002). Fig. 12 shows the pin-wise fission rate differences at 0 GWd/t. In the low enrichment fuel, the RMSs of pin-wise fission rate difference due to the reference subchannel effect and the model difference are 0.93% and 0.35%, respectively. In the high enrichment fuel, the respective RMSs are 1.36% and 0.93%. The model differences at both fuels are lower than the reference subchannel effect. In the low enrichment fuel, the higher reference subchannel effect is seen at the corner fuel rods, and the model differences at the corner fuel rods are lower than the reference subchannel effect. In the high enrichment fuel, the higher reference subchannel effect is seen at the fuel rods surrounding the water rod, and the model differences at the rods are lower than the reference subchannel effect. That is why the large subchannel effect is found along unheated surfaces such as the channel box walls and water rod surface because the subchannel void fractions decrease near the unheated surfaces. As a result, the developed model can decrease the reference subchannel effect, especially the difference of the pin power
1.50
Model difference (D) (AD)
1.00
(AR) Reference subchannel effect (R)
0.50 0.00 -0.50 -1.00 -1.50 -2.00 0.0
5.0 10.0 Exposure (GWd/t)
15.0
Fig. 10. Assembly power difference (high enrichment fuel).
37 42 42 40 38 38 40 41 35
42 48 46 44 40 40 45 47 41
42 46 42 38 37 33 39 45 40
40 44 38 33 33 30 33 40 39
38 40 37 33 33 37 40 38
38 40 33 30 33 33 38 44 40
40 45 39 33 37 38 42 46 42
41 47 45 40 40 44 46 48 42
35 41 40 39 38 40 42 42 37
Fig. 11. Subchannel void distribution of the high enrichment fuel (red values indicate the subchannels surrounding the Gd-containing rod).
82
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
0.83 -1.9 0.7
0.74 -1.2 0.7 1.10 -0.3 0.0
0.68 -0.6 0.4 1.27 0.2 0.4 1.19 0.7 0.4
1.09 -0.2 0.2 1.21 0.6 0.3 1.20 0.9 -0.3
Low
High
High
Low
1.07 0.2 0.4 1.18 0.8 0.2 1.18 0.9 0.0
1.09 -0.4 -0.2 1.19 0.4 0.2 1.13 0.7 -0.7 1.15 1.2 0.0 1.14 1.1 0.2 1.10 0.9 0.4
0.68 -1.3 -0.2 1.00 -0.3 -0.4 1.15 0.3 -0.4 1.11 0.4 -0.2 1.11 0.5 0.1 1.14 0.1 0.2 0.97 -0.5 0.0
0.73 -2.3 -0.1 0.65 -1.1 -0.2 1.00 -0.6 -0.4 0.95 -0.3 -0.3 0.96 -0.4 0.0 0.60 -0.8 0.1 0.64 -1.6 -0.6 0.71 -2.3 -0.3
Reference pin-wise fission rate Reference subchannel effect (PR) [%] Model difference (PD) [%]
1.26 1.3 1.2 1.36 1.0 0.7 1.33 1.1 1.4 1.24 0.6 0.9 1.28 1.1 2.0 1.38 1.4 0.6 1.42 1.5 0.5 1.34 1.0 0.3
1.28 0.9 -0.3 1.10 0.2 -0.3 0.98 -0.5 -0.2 0.24 0.0 -1.0 0.90 -0.4 1.0 1.08 0.3 0.6 0.00 0.5 -0.4
1.19 0.0 -0.5 0.94 -1.2 -0.3 0.24 -0.4 -0.4 0.88 -3.3 -0.2 0.94 -2.3 0.6 0.24 -0.4 -0.3
1.08 -0.6 -0.3 0.23 -0.4 -1.0 0.87 -3.1 0.3
1.10 -0.3 -1.0 0.82 -1.5 -0.6 0.91 -3.4 -1.8
1.18 0.3 -1.1 0.97 -0.4 -1.4 0.24 -0.4 0.1
1.21 0.9 -1.1 1.05 0.4 -1.2
1.13 1.1 -2.2
RMS difference [%] Low
High
Reference subchannel effect
0.93
1.36
Model difference
0.35
0.93
Fig. 12. Difference of pin-wise fission rate distribution.
Low enrichment fuel Channel box Fuel rod Water rod Water gap Gd-containing fuel rod
0.50 Difference of k-infinity (%dk)
Control rod
0.40 0.30 0.20 0.10 0.00 -0.10 -0.20
(kD) Model difference (D)
-0.30
(kR) Reference subchannel effect (R)
-0.40 -0.50 0.0
5.0 10.0 Average exposure (GWd/t)
15.0
Fig. 14. Difference of k-infinity at control rod inserted condition.
High enrichment fuel Boundary condition: mirror reflective in axial and radial boundaries Fig. 13. The 2D geometry of the heterogeneous problem with a control rod.
Table 3 Difference of k-infinity and assembly power of the heterogeneous problem with a control rod at 0 GWd/t. k-infinity (%dk) Case
Reference subchannel effect Model difference
0.37 0.21
Assembly power (%) Low with CR
High
1.61 0.39
0.88 0.46
Low without CR 1.13 0.94
at the fuel rods with large subchannel void differences, and the model has no negative influence. 3.5. Heterogeneous problem with a control rod Another evaluation of the proposed model is done with the control rod. The control rod is inserted during burnup. The calculated geometry is shown in Fig. 13. The average void fraction and axial power profile are the same as those used in Section 3.4. In the lattice physics calculation of the present evaluation, the burnup cal-
Bundle power difference (%)
2.00 1.50 1.00 0.50 0.00 -0.50 -1.00
(AD) Model difference (D)
-1.50
Reference subchannel effect (A (R)R)
-2.00 0.0
5.0 10.0 Exposure (GWd/t)
15.0
Fig. 15. Difference of assembly power (low enrichment fuel with a control rod at the north-west position).
culation without the control rod is performed in order to avoid a computational burden. The branch calculation is used to simulate the control rod inserted condition. These conditions are same as Table 2. Difference of k-infinity and assembly power at 0 GWd/t are summarized in Table 3. All differences of the reference subchannel effect are larger than the model difference at the beginning of burnup since the large power tilt is caused by the control rod and it makes lower void fraction around the control rod. The trend of kinfinity difference is shown in Fig. 14. The reference subchannel
83
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
2.00
1.50
(AD) Model difference (D)
1.00
Reference subchannel effect (R) (AR)
Bundle power difference (%)
Bundle power difference (%)
2.00
0.50 0.00 -0.50 -1.00 -1.50
1.50 1.00 0.50 0.00 -0.50 -1.00
(AD) Model difference (D)
-1.50
(AR) Reference subchannel effect (R)
-2.00
-2.00 0.0
5.0 10.0 Exposure (GWd/t)
0.0
15.0
effect on k-infinity is large at 0 GWd/t, and it decreases with burnup due to the negative feedback described in Section 3.4. But there is an undershoot until 0.12%dk at 12 GWd/t. The model difference for k-infinity is lower than the reference subchannel effect at 0 GWd/t and decreases by 5 GWd/t. This means that the model has good agreement with the reference analysis when burnup exceeds 5 GWd/t. 0.33 -12.9 0.4 0.61 -13.1 -0.8
0.35 -10.9 -1.4 0.83 -9.1 -0.9 0.93 -4.5 0.3
0.59 -7.6 -0.4 0.89 -5.5 -1.0 1.07 -2.0 0.6
Low High with CR High
0.65 -5.2 -1.1 0.98 -3.5 -0.8 1.16 0.0 -0.1
0.76 -2.8 0.4 1.09 -1.1 0.2 1.20 1.3 0.1 1.35 3.3 0.8 1.43 4.2 0.5 1.45 4.3 0.3
0.59 -2.7 0.2 1.03 -0.4 0.3 1.32 0.9 -0.7 1.38 3.0 0.2 1.45 4.2 0.3 1.56 4.4 0.3 1.37 4.3 0.7
Low w/o CR
15.0
Fig. 17. Difference of assembly power (low enrichment fuel at the south-east position).
Fig. 16. Difference of assembly power (high enrichment fuel).
0.32 -14.0 0.0
5.0 10.0 Exposure (GWd/t)
The respective changes with burnup of assembly power differences of low enrichment fuel with the control rod at the northwest position, high enrichment fuel, and low enrichment fuel at the south-east position are shown in Figs. 15–17. In low enrichment fuel with the control rod at the north-west position, the model difference is lower than the reference subchannel effect
0.80 -3.4 -1.3 0.76 -1.6 -1.3 1.24 0.5 -0.6 1.25 2.6 0.4 1.31 3.0 -0.2 0.84 3.4 0.7 0.93 3.1 0.9 1.06 2.2 0.9
1.09 -0.8 -1.0 1.19 -0.1 -0.8 1.17 0.6 0.4 1.12 1.7 1.1 1.19 2.1 1.6 1.31 2.8 0.6 1.38 3.1 0.8 1.33 2.9 1.2
1.16 -1.7 -2.8 1.00 -1.1 -1.7 0.90 -0.6 -0.6 0.23 0.0 -0.8 0.86 0.3 1.1 1.06 0.6 0.1 1.18 2.2 0.5 1.44 2.6 1.1
1.14 -2.9 -3.3 0.89 -2.8 -1.9 0.23 -0.9 -0.7 0.88 -3.3 -0.2 0.94 -2.7 0.3 0.25 -0.4 0.1 1.09 1.2 1.3 1.42 2.0 1.0
1.07 -3.1 -2.5 0.23 -0.9 -1.1 0.88 -3.7 0.1
1.11 -1.5 -1.8 0.84 -2.5 -1.3 0.93 -3.7 -1.8
0.97 -2.4 0.6 0.93 0.6 2.2 1.33 2.1 3.1
0.91 -2.2 1.2 0.25 0.4 -0.2 1.31 1.8 2.2
0.67 -2.5 -1.3
0.62 -2.4 -1.9 0.94 -1.2 -0.8
0.58 -1.2 -0.4 1.12 -0.4 -0.4 1.09 0.2 -0.4
0.94 -0.9 -0.6 1.10 -0.2 -0.6 1.13 0.5 -0.4
0.95 -1.0 -1.0 1.11 0.2 -0.3 1.15 0.6 -0.7
Reference pin-wise fission rate Reference subchannel effect (PR ) [%] Model difference (PD ) [%]
RMS difference [%]
Low with CR
High
Low w/o CR
Reference subchannel effect
5.35
1.84
0.99
Model difference
0.70
1.48
0.67
1.21 -1.1 -2.0 1.00 -1.5 -2.1 0.25 -0.4 0.7 0.94 -3.2 -1.2 0.90 -2.9 1.0 0.25 -0.4 0.2 1.03 0.5 1.1 1.41 1.7 2.4
1.24 -0.2 -1.9 1.08 -0.7 -1.9 1.00 -0.5 -1.0 0.85 -1.5 -0.3 0.24 0.0 0.0 0.98 -0.1 1.3 1.16 0.9 0.9 1.44 2.7 2.7
1.16 0.5 -2.4 1.25 0.2 -1.5 1.22 0.2 -0.7 1.15 -1.3 -1.6 1.13 -0.9 -0.2 1.25 -0.1 -0.2 1.35 1.7 0.9 1.34 2.2 2.4
0.99 -0.5 -0.3 1.14 1.0 0.3 1.13 0.8 -0.6 1.19 1.5 0.7 1.21 1.4 0.3 1.21 1.2 1.0
0.65 -0.9 0.0 1.00 0.2 0.2 1.20 0.4 0.3 1.20 0.5 -0.1 1.22 0.8 0.6 1.29 0.1 0.4 1.12 0.4 0.9
0.73 -1.8 0.6 0.68 -0.9 0.3 1.10 -0.4 -0.1 1.08 0.1 0.4 1.11 0.3 0.7 0.70 0.0 1.0 0.76 -1.2 0.8 0.85 -1.8 1.1
Fig. 18. Difference of pin-wise fission rate distribution of 0 GWd/t for control rod inserted condition.
84
T. Mitsuyasu et al. / Annals of Nuclear Energy 102 (2017) 77–84
during burnup, but it fluctuates. The pin power tilts due to the neighbor fuel assemblies change the subchannel void distribution from that in the infinite lattice. The model however calculates the lattice properties based on the subchannel void distribution made by pin powers in the infinite lattice. In this case, the model overestimates the assembly power in low enrichment fuel with the control rod, and causes undershoot later in the burnup due to the negative feedback. In high enrichment fuel, the model difference during burnup is lower than the reference subchannel effect. In low enrichment fuel at the south-east position, the model difference is almost the same as the reference subchannel effect except from 5 GWd/t to 12 GWd/t burnups. The duration of large difference in low enrichment at the south-east position agrees with the duration of the undershoot in low enrichment fuel with the control rod at the north-west position and the model difference is also estimated smaller than the difference of assembly power in high enrichment fuel. The difference is considered as a secondary effect because the total assembly power which is the summation of assembly powers is fixed and the underestimations should cause the overestimation. In the context of core analysis, the assembly with higher assembly power is important for evaluating the thermal margin such as the critical power ratio, and the present model can lower the difference in high enrichment fuel. Fig. 18 shows differences of the pin-wise fission rate at 0 GWd/t. In the low enrichment fuel with the control rod at the north-west position, the high enrichment fuel, and the low enrichment fuel at the south-east position, the respective RMSs of the reference subchannel effect are 5.35, 1.84 and 0.99%. And the respective RMSs of the model difference are 0.70, 1.48 and 0.67. The model can decrease the pin-wise fission rate difference at all assemblies, especially for the corner fuel rod neighboring the control rod. As a result, the model is effective for calculating the pin-wise fission rate with the control rod. 4. Conclusions The model of the two-stage core calculation method coupled with subchannel analysis for BWRs has been developed. The present model, which consists of the lattice physics analysis and the subchannel analysis parts, can be directly incorporated into the current core analysis model. In the present model, the lattice physics couples with the subchannel analysis which treats the 3D full assembly and uses the flow and power conditions to cover the maximum void fraction in the operation core. In order to simulate the power distribution inside a fuel assembly for the subchannel analysis, the cosine axial power profile was assumed. Through sensitivity calculations, the effect of the axial power profile was found to be negligible. The number of coupling iterations required for the model was three, which is sufficiently small. The present model
was validated in the heterogeneous problem with and without a control rod. The reference subchannel effect, which indicates the difference due to void distribution inside a fuel assembly, was always largest at the beginning of the burnup because the pin power tilt due to the neutronics heterogeneity was largest at the beginning of cycle and decreased as the burnup increased. The model difference, which is the difference between the reference subchannel effect and the model subchannel effect, was lower than the reference subchannel effect in terms of the k-infinity, assembly power, pin-wise fission rate and RMS. The present model was judged to be especially effective for the control rod inserted condition. The overall results indicated that the present model can incorporate the subchannnel effect into the current two-stage core calculation method. References Ama, T., Ikeda, H., Kosaka, S., Hyoudou, H., Takeda, T., 2002. Effect of moderator density distribution of annular flow on fuel assembly neutronic characteristics in boiling water reactor cores. J. Nucl. Sci. Technol. 39 (5), 487–498. Ama, T., Hyoudou, H., Takeda, T., 2002. Effect of radial void distribution within fuel assembly on assembly neutronic characteristics. J. Nucl. Sci. Technol. 39 (1), 90– 100. Hino, T., Ishii, K., Mitsuyasu, T., Aoyama, M., 2010. BWR core simulator using threedimensional direct response matrix and analysis of cold critical experiments. J. Nucl. Sci. Technol. 47 (5), 482–491. Ikehara, T., Kudou, Y., Tamitani, M., Yamamoto, M., 2008. Effect of Subchannel void fraction distribution on lattice physics parameters for boiling water reactor fuel bundles. J. Nucl. Sci. Technol. 45 (12), 1237–1251. Ishii, K., Hino, T., Mitsuyasu, T., Aoyama, M., 2009. Three-dimensional direct response matrix method using a Monte Carlo calculation. J. Nucl. Sci. Technol. 46 (3), 259–267. Jatuff, F., Giust, F., Krouthen, J., Helmersson, S., Chawla, R., 2006. Effects of void uncertainties on the void reactivity coefficient and pin power distributions for a 10 10 BWR assembly. Ann. Nucl. Energy 33 (2), 119–125. Katono, K., Nukaga, J., Fujimoto, K., Nagayoshi, T., Yasuda, K., 2015. Threedimensional time-averaged void fraction distribution measurement technique for BWR thermal hydraulic conditions using an X-ray CT system. J. Nucl. Sci. Technol. 52 (3), 388–395. Mitsuyasu, T., Ishii, K., Hino, T., Aoyama, M., 2010. Development of spectral history and pin power correction methods for pin-by-pin core analysis method using three-dimensional direct response matrix. J. Nucl. Sci. Technol. 47 (12), 1124– 1130. Mitsuyasu, T., Ishii, K., Nakadozono, N., Aoyama, M., A model of two-stage core calculation method coupled with subchannel analysis for boiling water reactors. In: Proc. PHYSOR-2014; 2014 Sep 28-Oct 4; Kyoto (Japan). [CD-ROM]. Morimoto, Y., Maruyama, H., Ishii, K., Aoyama, M., 1989. Neutronic analysis code for fuel assembly using a vectorized Monte Carlo method. Nucl. Sci. Eng. 103 (4), 351–358. Shibata, K., Iwamoto, O., Nakagawa, T., Iwamoto, N., Ichihara, A., Kunieda, S., Chiba, S., Furukawa, K., Otuka, N., Ohsawa, T., Murata, T., Matsunobu, H., Zukaran, A., Kameda, S., Katakura, J., 2011. JENDL-4.0: a new library for nuclear science and technology. J. Nucl. Sci. Technol. 48 (1), 1–30. Tada, K., Fujita, T., Endo, T., Yamamoto, A., Kosaka, S., Hirano, G., Nozaki, K., 2011. Application of quick subchannel analysis method for three-dimensional pin-bypin BWR core calculations. J. Nucl. Sci. Technol. 48 (12), 1437–1452. Tomiyama, A., Yokomizo, O., Yoshimoto, Y., Sugawara, S., 1988. Method of critical power prediction based on film flow model coupled with subchannel analysis. J. Nucl. Sci. Technol. 25 (12), 914–928.