Response Surface Methodology Control Rod Position Optimization of a Pressurized Water Reactor Core Considering Both High Safety and Low Energy Dissipation

Abstract: Response Surface Methodology (RSM) is introduced to optimize the control rod positions in a pressurized water reactor (PWR) core. The widely used 3D-IAEA benchmark problem is selected as the typical PWR core and the neutron flux field is solved. Besides, some additional thermal parameters are assumed to obtain the temperature distribution. Then the total and local entropy production is calculated to evaluate the energy dissipation. Using RSM, three directions of optimization are taken, which aim to determine the minimum of power peak factor Pmax, peak temperature Tmax and total entropy production Stot. These parameters reflect the safety and energy dissipation in the core. Finally, an optimization scheme was obtained, which reduced Pmax, Tmax and Stot by 23%, 8.7% and 16%, respectively. The optimization results are satisfactory.


Introduction
In China, energy consumption has grown rapidly driven by the improvement of the economic level.In order to balance economic development and environmental protection, China's government has established a green development strategy.An important measure is to adjust the energy structure and promote the consumption of clean energy sources such as nuclear energy.Thus, the use of nuclear energy increased from 44.19 billion kWh in 2000 to 290.75 billion kWh in 2013 [1], and 26 nuclear power units (28,528 MWe) are currently under construction in China [2,3].
The design of a nuclear reactor system includes shielding design, thermodynamics, fluid flow and heat transfer, fuel element design, radioactivity releases, etc. [4].In nuclear physics calculations, the neutron flux density and its distribution are usually the primary target.The Neutron Transport Equation (NTE) is used to describe the variation of neutron flux density in a physical field.For simplification, the Neutron Diffusion Equation (NDE) is also used.These equations can be solved by the Finite Different Method (FDM) [5][6][7], the Finite Volume Method (FVM) [8][9][10][11], the Finite Element Method (FEM) [10,[12][13][14], etc. Normally, the neutron flux density distribution matches the fission heat source distribution.
Heat transfer is also an important physical process in nuclear reactor core.It is an irreversible process which increases the entropy of the whole system.The differences in heat transfer in different regions causes a distribution of local entropy production.In 1979, Bejan [15] presented the fundamental equations describing entropy production due to finite temperature gradients and fluid flow.Then Bejan [16][17][18] discussed the irreversibility in thermodynamic systems by entropy production analysis.Since then, many researchers have started to use the entropy production method to evaluate heat transfer processes in various structures.For example, Ibáñez et al. [19] investigated the entropy production minimization of a solid slab with uniform internal heating and asymmetric convective cooling.Makinde and Aziz [20] analyzed the inherent irreversibility and thermal stability in the model of a long hollow cylinder with temperature-dependent internal heating and asymmetric convective cooling boundaries.Aziz and Khan [21] analyzed the entropy production of steady conduction in the model of a hollow sphere with temperature-dependent internal heating and asymmetric convective cooling boundaries.Malvandi et al. [22] studied the entropy production of steady two-dimensional boundary layer flow of nano-fluids over a flat plate.Torabi and Zhang [23] used classical entropy production analysis to investigate heat transfer in cooled homogenous and functionally graded material slabs with variation of internal heat generation with temperature, and convective-radiative boundaries.
For coupling the neutronic and thermal-hydraulic calculation, the traditional technique is the Operator Splitting (OS) method [24,25].By this method, different physical processes are solved separately.The output of neutronic calculation is the input of thermal-hydraulic part.Some other advanced multi-physics algorithms have also been developed.One of the most famous algorithms is the Jacobian-Free Newton-Krylov (JFNK) method [26,27].By the JFNK method, several physical processes are treated together as an entire unit.The biggest advantage of the OS method is its ease of implementation, while JFNK is more suitable for nonlinear problems and has higher accuracy.
Response Surface Methodology (RSM) is a collection of mathematical and statistical techniques used to explore the relationships between several independent variables and one or more response variables.This methodology was proposed by Box and Wilson in the 1950s [28,29].Since then this methodology has been widely used in analytical chemistry [28], bioprocessing [29,30], structural reliability [31,32], food chemistry [33][34][35], etc.The biggest advantage of this methodology is that the RSM model is easy and convenient to establish, even when little information about the process can be obtained.
This paper aimed to optimize the insertion positions of the control rods in a Pressurized Water Reactor (PWR).The neutron diffusion and heat transfer processes are coupled and both high safety and low energy dissipation are taken into consideration.The chosen safety factors are the power peak factor and maximum temperature in the core, and the energy dissipation is calculated by the entropy production method.RSM is utilized to find a more improved scheme.

Description of Problem
The International Atomic Energy Agency (IAEA) has published a three-dimensional PWR benchmark problem named the 3D-IAEA [36].It is a typical PWR problem which is widely used in verification of neutron diffusion calculation codes.This PWR is taken as a study object in this paper to demonstrate the RSM optimization of control rod insertion positions.The nuclear data given in [36] are used to solve the neutron flux distribution.Some essential thermal parameters are also assumed to calculate the temperature field and local and total entropy production.
There are 177 groups of Fuel Assemblies (FAs) arranged in the reactor core, and 13 groups of these FAs are have inserted control rods.Five kinds of materials are applied in this core.A quarter of horizontal and one half of vertical cross section of the 3D-IAEA problem are illustrated in Figures 1 and 2, respectively.The multi-group neutron diffusion approximation that most widely used for commercial reactors is the two-energy groups approximation [37,38], which is also applied in 3D-IAEA benchmark problem.The neutron energy spectrum is divided into two groups, which are called fast and thermal neutrons, respectively.The two-group nuclear data are given in Table 1.Some additional thermal parameters are given in Table 2.

Neutron Diffusion
The NTE based on transport theory can describe the neutron behaviors in more details, but to balance of amount of calculation and perduring accuracy, in the engineering area the diffusion theory and NDE are always used to solve the neutron flux, which can be written as: where g is the index of energy group, g = 1,2, . . .,G; G is the total number which the energy spectrum is divided into; v is the velocity of neutrons, cm•s −1 ; t is the time, s; D is the diffusion constant, cm; Σ t is the total cross section, cm −1 ; Σ s,g'→g is the scattering cross section from g' energy group to g energy group, cm −1 ; Σ f is the fission cross section, cm −1 ; → r is the position, cm; ϕ is the neutron scalar flux density, neutrons•cm −2 •s −1 ; υ is the average number of neutrons that are emitted from each fission process; χ is the fission spectrum function; k eff is the multiplication factor; S e is the extra neutron source, neutrons•cm −3 •s −1 .
When using in time independent problem with no extra neutron source and no upscattering neutrons, Equation (1) can be simplified as: Thus, for the two-group diffusion problem, the governing equations can be written as: where Σ a is the absorb cross section, cm −1 , and ∑ t,1 ).The local fission reaction rate per unit volume is: where R f is the fission reaction rate, reaction times•cm −3 •s −1 .
It can be assumed that the total heat is generated at the position where the fission reaction occurs, therefore the power of local fission reaction can be defined as [39]: where E f is the total heat generation in one fission reaction.
According to the thermal parameters shown in Table 2, the reactor core is designed at the normal thermal power P core , which means that under the normal conditions, the summary of local powers is P core : Average subassembly power P k is defined as: where k is the designated number of FA which can be seen in Figure 1.
In order to compare our results with Ref. [36], a type of normalized neutron fluxes ϕ g ( → r ) is introduced.In the normalized form, the volumetrically weighted average of released fission neutron fluxes in the reactor core is set as 1, which can be written as: The aim of normalization is to eliminate the impact of initial conditions on the neutron flux distribution.

Heat Transfer and Energy Dissipation
The governing equation of conductive heat transfer can be written as: where ρ is the density, kg•m −3 ; c p is the specific heat capacity, J•kg T is the temperature, K; P core is the total thermal power of the reactor core, W. For the time independent problem, Equation ( 7) can be simplified as: The irreversible energy dissipation can be calculated by the entropy production.Normally, the local entropy production of conductive heat transfer can be written as [17,18]: The total entropy production is:

Response Surface Methodology
By RSM, if all the independent variables can be measured, the response surface can be expressed as: where X i represents the independent variable, Y is the response variable.
It is assumed that the independent variables are continuous.The goal is to find a suitable approximate relationship between independent variables and the response variable.Usually, a second-order model is utilized [40,41], which can be written as: Entropy 2017, 19, 63 where β is the undetermined coefficient, ε is a random error.Equation ( 14) can be written in matrix form: The solution of Equation ( 15) can be obtained by the matrix approach: In Equation ( 16), the superscript T represents matrix transposition, the superscript −1 represents matrix inversion.
There are several two-order designs for RSM.The Central Composite Design (CCD) and Box-Behnken Design (BBD) are two typical ones.In this paper, BBD is adopted to get the coefficients of Equation ( 14).Based on the symmetry of the reactor core along x and y axis, only one-eighth of the whole core need to be designated, which can be seen in Figure 1.Thus, the insertion positions of four control rods can be chosen as independent variables, namely Z 1 , Z 2 , Z 3 and Z 4 , which represent the positions of FA.1, FA.18, FA.5 and FA.31, respectively.

Standard Problem Solution
The control rod insertion positions shown in Figure 2 were chosen as the standard and reference solution.This example solution is the basis of the optimization.The neutron diffusion equations and heat transfer equations are all solved by FVM.In each calculation, the considered domain is one quarter of the whole core with reflective and non-return external boundaries.The convergence criterion is set as maximum relative flux change on each inner iteration = 10 −10 , maximum k eff change on outer iteration = 10 −6 .Four groups of regular hexahedral meshes are taken into consideration.The mesh quantities and calculation results are given in Table 3.The error in Table 3 represents the relative error between using a mesh group and its further refinement.Considering both of the calculation accuracy and time cost, mesh No. 2 was selected to be applied for the next optimization scheme.This quantity of mesh is also used in [36], thus our neutron diffusion calculation results can be checked.The calculated k eff is 1.02855 and the reference at the same mesh quantity is 1.02864, with the difference of 0.0087%.The results of fast and thermal neutron flux at the diagonal line on the x-y plane at the level of z = 195 cm are plotted in Figures 3 and 4. A good match between calculated results and references can be seen.Considering both of the calculation accuracy and time cost, mesh No. 2 was selected to be applied for the next optimization scheme.This quantity of mesh is also used in [36], thus our neutron diffusion calculation results can be checked.The calculated keff is 1.02855 and the reference at the same mesh quantity is 1.02864, with the difference of 0.0087%.The results of fast and thermal neutron flux at the diagonal line on the x-y plane at the level of z = 195 cm are plotted in Figures 3  and 4. A good match between calculated results and references can be seen.The calculation local power, temperature and entropy production results are plotted in Figures 5-7.It can be seen in Figure 5 that the control rods have an obvious effect on the neutron fission power.Both the vertical and horizontal cross sections show significant decreases of local power depending on the control rod positions.The comparison of horizontal cross section at different z levels shows that the local power near the center of the core is higher than at the upper level, as expected.Figure 6 shows that the regular temperature distribution is not as complex as the fission power distribution.It shows a continuous reduce from the core center to the edge.The control rod positions cannot be seen clearly in these temperature cloud pictures.The regular local entropy production distribution is far from the temperature distribution.The calculation local power, temperature and entropy production results are plotted in Figures 5-7.It can be seen in Figure 5 that the control rods have an obvious effect on the neutron fission power.Both the vertical and horizontal cross sections show significant decreases of local power depending on the control rod positions.The comparison of horizontal cross section at different z levels shows that the local power near the center of the core is higher than at the upper level, as expected.Figure 6 shows that the regular temperature distribution is not as complex as the fission power distribution.It shows a continuous reduce from the core center to the edge.The control rod positions cannot be seen clearly in these temperature cloud pictures.The regular local entropy production distribution is far from the temperature distribution.
Figure 7 shows that the boundary centers of the reactor core are the peaks of the local entropy production, and the volume center of the reactor core is the valley of the local entropy production.The reason is that this core is symmetric at the vertical middle cross section, thus the middle cross section is heat insulated and the temperature gradient there is zero.The rapid decrease of temperature near the boundary center causes the entropy production peak.
depending on the control rod positions.The comparison of horizontal cross section at different z levels shows that the local power near the center of the core is higher than at the upper level, as expected.Figure 6 shows that the regular temperature distribution is not as complex as the fission power distribution.It shows a continuous reduce from the core center to the edge.The control rod positions cannot be seen clearly in these temperature cloud pictures.The regular local entropy production distribution is far from the temperature distribution.In total, some lumped and feature parameters are selected to describe the thermal situation of the reactor core.These parameters are the power peak factor P max , the maximum temperature T max and the total entropy production S tot .These parameters' values are shown in Table 4.

Response Surface Design
The BBD is introduced to apply the RSM.The insert positions of control rods Z 1 , Z 2 , Z 3 , Z 4 are the input independent variables, and the feature parameters P max , T max , S tot are chosen as the response variables.The RSM can be written as follows: where: In order to establish the coefficient matrix b, several experiments should be performed.Section 4.1 presents an example of a numerical experiment.The input variables are [Z 1 Z 2 Z 3 Z 4 ] = [202802020], and the response result Y is given in Table 4. Using BBD to design the RSM for this problem, 29 groups of numerical experiments must be executed.The input and response variables are shown in Table 5.The matrix b is obtained by Equation ( 17) and the fitting results are shown in Table 6.The fitting results are usually evaluated by R-square.For this evaluation index, the closer to 1, the better result gained.The R-square value of the three response surface fitting is shown in Table 7.It can be seen in Table 7 that, although T max is a local feature, depending on its simple regular distribution, the R-square of its response surface is the highest.The second is S tot since it is a total summary of local entropy production.It has the both global and local parameter features, so its sensitivity to the input variables is weaker than that of the local parameter.The R-square of the P max response surface is not as fine as the others, due to the complex local power distribution.As these response surfaces are introduced to find out the direction of optimization, rather than obtain the exact solution, this level of fitting can also be applied.

Rod Position Optimization
In this section, a simple code is programmed to traverse all four input variable ranges in the definition domain, which is set as [30 cm, 300 cm].The coefficients in Table 6 and Equation (18) are used in this code to obtain the approximate response results.The interval step of the input variables is set as 10 cm.In every traversal calculation, one evaluation index (P max , T max or S tot ) is selected as the optimization direction.The aim of optimization is to promote safety or reduce entropy dissipation, that is to say, P max , T max or S tot should be reduced.
Although these three indexes are key metrics for a reactor core, the importance rankings are not at the same level for each of them.For every nuclear energy system, safety is always a seriously concerned for the public [42], and its status is higher than economical efficiency.T max is the most obvious safety index, as the reactor materials may melt above a certain temperature limit, which would cause a serious accident.P max is also a type of safety index, which is associated with the level of power flattening.A more flattened power distribution means a smoother reactor operation is expected.The last is the economical index S tot , which represents the energy dissipation.
Therefore, our optimization procedure is: (1) use RSM to traverse all the control rod positions to find N 1 positions approaching the lowest T max .The selected positions are put in set Ω 1 ; (2) Use RSM again to traverse all the positions in the definition domain to find N 2 positions approaching the lowest P max and name the selected positions set as Ω 2 ; (3) Use FVM to obtain accurate indexes of position elements in Ω 3 , Ω 3 = Ω 1 ∩ Ω 2 , and sort them by S tot .The best solution is chosen as the final optimization scheme.The flow diagram of this optimization procedure is shown in Figure 8.In this scheme, a smaller N1 or N2 means a narrower search range, but the indexes may become closer to the limit.In this case, they are set as N1 = N2 = 10.Finally, eight groups of control rod positions are selected, which are listed in Table 8.It can be seen in Table 8 that the safety and entropy index will reach optimization when the central control rods are inserted deeply and the peripheral control rods are inserted shallowly.In the calculation series, it can be found that when the input rod positions are [30 30 300 300], the Pmax, Tmax and Stot reach the minimum, which means that both the power and temperature are flatter, and the energy dissipation reaches the minimum.In order to observe the physical field in depth, the In this scheme, a smaller N 1 or N 2 means a narrower search range, but the indexes may become closer to the limit.In this case, they are set as N 1 = N 2 = 10.Finally, eight groups of control rod positions are selected, which are listed in Table 8.It can be seen in Table 8 that the safety and entropy index will reach optimization when the central control rods are inserted deeply and the peripheral control rods are inserted shallowly.In the calculation series, it can be found that when the input rod positions are [30 30 300 300], the P max , T max and S tot reach the minimum, which means that both the power and temperature are flatter, and the energy dissipation reaches the minimum.In order to observe the physical field in depth, the distributions of local power, temperature and local entropy production are illustrated in Figures 9-11.In order to show some details of the optimization results, the local power, temperature and local entropy production at the diagonal line on midplane are shown in Figure 12.It can be seen in Figure 12a that the local power decreases near the core center and a power valley forms between the distance of 40 cm and 50 cm, as the control rods in FA.1 and FA.18 are inserted deeply.At the same time, with the control rods in FA.5 and FA.31 rise, the local power of the optimization scheme is higher than that of the standard problem in the distance range of [65 cm, 110 cm].As a result, local power is flatter than in the standard problem and the power peak factor is reduced.For the temperature distribution, as shown in Figure 12b, the temperature and its gradient are reduced near the core center while a little increase is noted near the core boundary, which matches the local power distribution.The change of temperature gradient causes a lower local entropy production near the core center and the one higher near the core boundary, which is shown in Figure 12c.The cloud pictures Figures 9 and 11 share the same colorbar scale as Figures 5-7, respectively.Thus the differences between these two groups can be clearly observed.Comparing Figure 9 with Figure 5, it can be found that the power flattening effect is better after the control rod adjustment.The central control rod has an obvious influence on the power distribution, which causes a power valley at the core center, so the temperature at the center line of x = y = 0, is not so high as the temperature in the standard problem, which can be seen in Figures 6 and 10.The distribution of local entropy production is similar between Figures 7 and 11, but the entropy production values are reduced after control rod adjustment.
In order to show some details of the optimization results, the local power, temperature and local entropy production at the diagonal line on midplane are shown in Figure 12.It can be seen in Figure 12a that the local power decreases near the core center and a power valley forms between the distance of 40 cm and 50 cm, as the control rods in FA.1 and FA.18 are inserted deeply.At the same time, with the control rods in FA.5 and FA.31 rise, the local power of the optimization scheme is higher than that of the standard problem in the distance range of [65 cm, 110 cm].As a result, local power is flatter than in the standard problem and the power peak factor is reduced.For the temperature distribution, as shown in Figure 12b, the temperature and its gradient are reduced near the core center while a little increase is noted near the core boundary, which matches the local power distribution.The change of temperature gradient causes a lower local entropy production near the core center and the one higher near the core boundary, which is shown in Figure 12c.
is higher than that of the standard problem in the distance range of [65 cm, 110 cm].As a result, local power is flatter than in the standard problem and the power peak factor is reduced.For the temperature distribution, as shown in Figure 12b, the temperature and its gradient are reduced near the core center while a little increase is noted near the core boundary, which matches the local power distribution.The change of temperature gradient causes a lower local entropy production near the core center and the one higher near the core boundary, which is shown in Figure 12c.As a result, the evaluation indexes of the standard problem and optimization scheme are listed in Table 9.An obvious decrease can be seen, which is satisfactory for the optimization.

Conclusions
In this paper, a typical PWR reactor core based on the 3D-IAEA problem is taken to be analysed.The neutron flux field, temperature and local entropy production distribution in this reactor core are calculated by FVM.The calculation results are illustrated in cloud pictures in order to observe them clearly.For evaluating the safety and energy dissipation, three characteristics are selected as evaluation indexes, which are the power peak factor P max , the maximum temperature T max , and the total entropy production S tot .Then the RSM is introduced to optimize the control rod insertion positions in order to get a lower P max , T max and S tot , which means higher safety and lower energy dissipation.Three directions of optimization are taken, and a final optimization scheme is obtained.The comparison of cloud pictures between the standard problem and optimization scheme shows that the central control rod has a great influence on the power distribution.In total, P max , T max and S tot are reduced by 23%, 8.7% and 16% after the adjustment of the control rod insert positions, which is satisfactory.RSM may be a quick and useful method in reactor optimization.

Figure 1 .
Figure 1.Horizontal cross section of the 3D-IAEA problem.

Figure 2 .
Figure 2. Vertical cross section of the 3D-IAEA problem.

Figure 2 .
Figure 2. Vertical cross section of the 3D-IAEA problem.

Figure 3 .
Figure 3. Fast neutron flux at the diagonal line at the level of z = 195 cm.Figure 3. Fast neutron flux at the diagonal line at the level of z = 195 cm.

Figure 3 .
Figure 3. Fast neutron flux at the diagonal line at the level of z = 195 cm.Figure 3. Fast neutron flux at the diagonal line at the level of z = 195 cm.Entropy 2017, 19, 63 7 of 16

Figure 4 .
Figure 4. Thermal neutron flux at the diagonal line at the level of z = 195 cm.

Figure 4 .
Figure 4. Thermal neutron flux at the diagonal line at the level of z = 195 cm.

Figure 5 .
Figure 5. Local power distribution of standard problem (left: vertical cross section cloud picture at y = 0; right-top: horizontal cross section cloud picture at z = 315 cm; right-bottom: horizontal cross section cloud picture at z = 195 cm.Pavg represents the average local power of standard problem).

Figure 5 .Figure 6 . 3 Figure 6 .
Figure 5. Local power distribution of standard problem (left: vertical cross section cloud picture at y = 0; right-top: horizontal cross section cloud picture at z = 315 cm; right-bottom: horizontal cross section cloud picture at z = 195 cm.P avg represents the average local power of standard problem).Entropy 2017, 19, 63 8 of 16

Figure 6 .Figure 7 .
Figure 6.Temperature distribution of standard problem (left: vertical cross section cloud picture at y = 0; right-top: horizontal cross section cloud picture at z = 315 cm; right-bottom: horizontal cross section cloud picture at z = 195 cm).

Figure 7 .
Figure 7. Local entropy production distribution of standard problem (left: vertical cross section cloud picture at y = 0; right-top: horizontal cross section cloud picture at z = 315 cm; right-bottom: horizontal cross section cloud picture at z = 195 cm.S 0 represents the average local entropy production of standard problem).

Figure 8 .
Figure 8. Flow diagram of optimization procedure.

Figure 9 .
Figure 9. Local power distribution of optimization scheme (left: vertical cross section cloud picture at y = 0; right-top: horizontal cross section cloud picture at z = 315 cm; right-bottom: horizontal cross section cloud picture at z = 195 cm.Pavg represents the average local power of standard problem).

Figure 10 .
Figure 10.Temperature distribution of optimization scheme (left: vertical cross section cloud picture at y = 0; right-top: horizontal cross section cloud picture at z = 315 cm; right-bottom: horizontal cross section cloud picture at z = 195 cm).

Figure 9 .Figure 9 .
Figure 9. Local power distribution of optimization scheme (left: vertical cross section cloud picture at y = 0; right-top: horizontal cross section cloud picture at z = 315 cm; right-bottom: horizontal cross section cloud picture at z = 195 cm.P avg represents the average local power of standard problem).

Figure 10 .
Figure 10.Temperature distribution of optimization scheme (left: vertical cross section cloud picture at y = 0; right-top: horizontal cross section cloud picture at z = 315 cm; right-bottom: horizontal cross section cloud picture at z = 195 cm).

Figure 10 .
Figure 10.Temperature distribution of optimization scheme (left: vertical cross section cloud picture at y = 0; right-top: horizontal cross section cloud picture at z = 315 cm; right-bottom: horizontal cross section cloud picture at z = 195 cm).

Figure 11 .
Figure 11.Local entropy production distribution of optimization scheme (left: vertical cross section cloud picture at y = 0; right-top: horizontal cross section cloud picture at z = 315 cm; right-bottom: horizontal cross section cloud picture at z = 195 cm.S0 represents the average local entropy production of standard problem).

Figure 12 .
Figure 12.Local power, temperature and local entropy production of the standard problem and optimization scheme at the diagonal line on the midplane (the level of z = 195 cm).

Figure 11 .
Figure 11.Local entropy production distribution of optimization scheme (left: vertical cross section cloud picture at y = 0; right-top: horizontal cross section cloud picture at z = 315 cm; right-bottom: horizontal cross section cloud picture at z = 195 cm.S 0 represents the average local entropy production of standard problem).

Figure 12 .
Figure 12.Local power, temperature and local entropy production of the standard problem and optimization scheme at the diagonal line on the midplane (the level of z = 195 cm).Figure 12. Local power, temperature and local entropy production of the standard problem and optimization scheme at the diagonal line on the midplane (the level of z = 195 cm).

Figure 12 .
Figure 12.Local power, temperature and local entropy production of the standard problem and optimization scheme at the diagonal line on the midplane (the level of z = 195 cm).Figure 12. Local power, temperature and local entropy production of the standard problem and optimization scheme at the diagonal line on the midplane (the level of z = 195 cm).

.
Two-group nuclear data.
. Two-group nuclear data.

Figure 2. Vertical cross section of the 3D-IAEA problem.Table 1 .
Two-group nuclear data.

Table 3 .
Results on different mesh quantity.

Table 4 .
Calculation results of standard problem.

Table 5 .
Numerical experiments used to design response surface.

Table 6 .
Coefficients in the response surface.

Table 7 .
R-square of response surface fitting.

Table 8 .
The selected control rod positions.

Table 8 .
The selected control rod positions.

Table 9 .
Comparison between calculation results.