Abstract

CFD simulation for a PWR is an important part for the development of Numerical Virtual Reactor (NVR) in Harbin Engineering University of China. CFD simulation can provide the detailed information of the flow and heat transfer process in a PWR. However, a large number of narrow flow channels with numerous complex structures (mixing vanes, dimples, springs, etc.) are located in a typical PWR. To obtain a better CFD simulation, the challenges created by these structural features were analyzed and some quantitative regularity and estimation were given in this paper. It was found that both computing resources and time are in great need for the CFD simulation of a whole reactor. These challenges have to be resolved, so two schemes were designed to assist/realize the reduction of the simulation burden on resources and time. One scheme is used to predict the combined efficiency of the simulation conditions (configuration of computing resources and application of simulation schemes), so it can assist the better choice/decision of the combination of the simulation conditions. The other scheme is based on the suitable simplification and modification, and it can directly reduce great computing burden.

1. Introduction

Advanced simulation tools based on the latest physical models, computing condition, and numerical technology have been the research focuses in the international scope, such as VERA which is being developed by CASL in USA [1] and NERESIM/NURISP/NURESAFE [2, 3] which were developed by the cooperative countries in Europe. The development of these advanced products relies on the fine-scale multiphysics computational models [4]. With the same objective, Numerical Virtual Reactor (NVR) is another advanced simulation tool which is being developed by Harbin Engineering University in China, and NVR is expected to be used for the design, safety analysis, operation research, and education and training.

CFD simulation of a reactor is an important part for the development of NVR because CFD simulation can reflect the detailed flow and heat transfer which can be used in the further application or research. For instance, with the application of CFD simulation, In et al. [5] analyzed the principle of heat transfer enhanced by spacer grid with mixing vane, Kim and Seo [6, 7] gave the optimal design of the shape of mixing vane, Lee and Choi [8] study the influence of the region size of the research object, and so on. Though CFD simulation is so important, some challenges still exist in the application.

Firstly, a successful CFD simulation depends on the mesh quality and the size of a research object, which always bring great computing burden. In a typical pressurized water reactor, there are thousands of narrow flow channels and complex structures such as the spacer grid and mixing vane, so that great cell number is needed. Thus the computing resource is also great, or the research object is cut into a local region with fewer rods. Then the computing resource can be satisfied easily; however the region size can affect the accuracy of the simulation [8]. To obtain the high fidelity results, simulation for the whole or the large region of a reactor has to be done. So some optimal methods or schemes should be designed to reduce the computing burden.

Secondly, simulation time is also too long when the cell number is large because of the iterative numerical technology in a CFD computation. Though some optimal methods or schemes can reduce a certain cell number, the total cell number is still very large for a whole reactor or a large region of a reactor. The researches of physical process and the engineering application using CFD simulation are also impossible under low simulation efficiency, so the optimal methods or schemes are also needed to improve the efficiency.

Thirdly, the accuracy can be tested by the experiments; however it is confusing to have a high efficient CFD simulation without a quantitative guidance on the design and choice of simulation conditions which contains the configuration of computing resources and the design of simulation schemes (meshing methods, solution methods, physical models, etc.). Usually, a researcher has to try and compare many times to find a suitable combination. Thus lots of time is lost and much computing resources are occupied during this process. Furthermore, the combination set for the simulation is always just suitable, but not optimal. So the judgment standards are also needed for the efficiency improvement.

For the reasons above, this paper tries to obtain () the clear understanding of the challenges for the CFD simulation of a PWR reactor; () the quantitative estimation to choose the high efficient simulation conditions; () the advanced methods and schemes to reduce the burden in both computing resources and time.

2. Research Objects and Approaches

2.1. Research Objects

In the application of CFD simulation, the first and important step is the mesh analysis which is closely related to the accuracy. However, in most previous researches, the mesh analysis is only done for the whole research object, and the last chosen mesh may be still not good because poor quality mesh can exist in the local feature region. Many feature regions are in a typical PWR, so this research began from each feature region.

A typical PWR has hundreds of fuel assemblies, and each assembly contains hundreds of fuel rods; then numerous narrow flow channels are created. Fortunately, similar feature structures exist, and they are the simple subchannel (SC), spacer grid (SG), and mixing vanes (MV), as shown in Figures 1(a)~1(c). The geometry of SG is more complex. To suit the small local region, SG is divided into two regions (regions around dimple (DI) and spring (SP) as shown in Figures 1(d)~1(e)).

In addition, the number of the feature structures (SC, MV, DI, and SP) is dozens of thousands. Therefore, the cell number saved in one local region will be enlarged thousands of times for the whole reactor. So it is significant to find the better meshing method for each feature structure.

The geometrical parameters of the SC and MV come from the work of Karoutas et al. [9], as shown in Figure 2. The parameters of this object can be found in the reference from Navarro and Santos [10]. The object of SG is a part of AFA2G grid as shown in Figure 3, and the parameters can be obtained from Ma [11].

2.2. Numerical Approaches

In this paper, the research was performed by ANSYS Fluent 15.0. Three RANS models (SKE, RNG, and RSM) were used. The standard wall function method was used and the in the first layer of the mesh near the wall is kept between 30 and 60 for these RANS models. The spatial discretization methods in the simulations are second order for pressure term and second order upwind for other terms. The convergence criteria of each term are set as . The boundary condition is shown in Table 1. The physical properties of the coolant are calculated by IAPWS-IF97.

2.3. Analysis of Meshing Methods

To avoid the low accurate simulation, the mesh independent analysis was done for each meshing method (mesh type) for each feature region. Then the minimum cell number can be obtained for the meshing method for each feature region. Then more suitable meshing method can be found by the comparison among the minimum cell numbers for each feature region.

In addition, the pressure distribution, lateral flow distribution, the total simulation time (TST), the mean time per step (MTPS), the iteration step number (ISN), the cell number (CN), and so on were recorded, compared, and analyzed for each meshing method to obtain some regularities.

2.3.1. Mesh Analysis for Simple Subchannels

The simple subchannel region has the same crosswise structure along the vertical direction, so prismatic mesh types were compared since their high mesh quality and low need of cells. The triangular prism (TP), polyhedral prism (PP), and hexahedral prism (HP) were chosen for the subchannel region, as shown in Figures 4(a)4(c). In the mesh establishments, the same axial mesh size was used firstly to find the suitable crosswise cell size; then various axial sizes were used to find the suitable axial size.

The simulation results created by different mesh types and crosswise cell sizes are shown in Figure 5. Table 2 gives the simulation time for each mesh type with the least cell number. In the comparison, it can be seen that each mesh type can get the same result when the cell number is over a certain value, and both the cell number and the simulation time needed by hexahedral prism mesh are the least. Hexahedral prism also can be established by many meshing methods as shown in Figures 4(c)4(e). Hexa-prism2 (HP2) and hexa-prism3 (HP3) are two structured mesh types which were widely used in the previous simulations, such as the work of Horváth and Dressel [12] and Nematollahi and Nazifi [13]. Hexaprism (HP) is an unstructured mesh type. It is always known that structured types are better for the data storage and transfer, so the simulation speed of structured mesh is always faster and the iteration step number is smaller. However, HP2 and HP3 need more cell number from the comparison in Figure 6 and Table 2. It can be seen that the unstructured meshing method is more efficient with less cell number and simulation time. The reason may be that the unstructured mesh has strong geometrical adaption for the structural feature of the simple subchannel region; however the quality of the structured mesh will decline obviously when the cell number is low for this structural feature. Meanwhile, the unstructured mesh is also simple and convenient for the establishment, but the structured mesh types need more experience and time to design and modify. Furthermore, the establishment of the unstructured mesh is easier for a large region like the whole reactor. So it is better to choose HP for the simple subchannel region.

After the comparison of different mesh types using various crosswise cell sizes, the comparison of various axial sizes was done for HP mesh with the optimal crosswise cell size. As shown in Figure 7, HP4 (the red point in Figure 7) is a HP mesh with the optimal axial cell size.

From the comparison in Table 2, the suitable cell number is from 133k to 9k, and the simulation time is from 198 s to 2 s. This is also an evidence to show the significance of the optimal meshing method to reduce the simulation cost.

2.3.2. Mesh Analysis for Complex Structure Regions

The regions surrounding mixing vane, dimple, and spring have complex geometrical features. The polyhedral (PH) mesh was used because of its high adaptive ability, and PH has been proved better for these regions in PWR by Li and Gao [14]. Each mesh is shown in Figure 8.

The mesh independent analysis is shown in Figure 9, and Table 3 gives the least cell number and the least simulation time for these complex structure regions. In the analysis, only one single feature region was considered, and the cell number is still not small. Then the simulation burden will be enlarged obviously when the simulation object is a large region with a large number of complex regions.

2.3.3. Mesh Analysis for 5 × 5 AFA 2G Grid

Based on the research of meshing method for each feature region, the meshes were built and checked for a 5 × 5 AFA 2G grid, and the mesh type is shown in Figure 10.

As the previous estimation, the cell number for the inner region of 5 × 5 AFA 2G grid is nearly ((103 × 2 + 229)/2 + 85) × 5 × 5k (≈7.56 million). The mesh analysis was done for both pressure drop and lateral flow status as shown in Figures 11 and 12. It can be found that the difference on resistance distribution is not obvious; however 7.98 M CN can obtain better lateral flow rate. So the suitable CN is very near the previous estimation.

2.4. Validation on the Simulation

Besides the mesh independent analysis, the validation is also an important step for the application. In a validation, the mesh type, cell density, turbulence model, design of boundary layer, operation of geometry, numerical discrete scheme, solution scheme, and so on can be tested to see whether they are suitable or optimal.

The tested object and boundary condition come from the work of Karoutas et al. [9]. The mesh type of each feature region for the validation is shown in Figure 13. The hexahedral prism is used for the simple subchannel. The polyhedron is used for the spacer grid. The mesh sizes are based on the previous mesh independent analysis.

The spatial discretization methods are introduced in Section 2.2. The main boundary conditions are listed in Table 1. As shown in Figure 14, PB1 and PB2, PB3 and PB4, and PB5 and PB6 are designed with periodic boundaries. The transverse flow status is important for the mass, momentum, and energy exchange between adjacent coolant channels. Therefore, this status is chosen for the comparison and the parameter is defined as the rate of the lateral velocity and the bulk velocity. The cross sections with vertical position of 12.7 mm and 101.6 mm were chosen for the comparison. The red path line of each section for the comparison was designed as Figure 14. The lateral velocity is normal to the red pathline.

As shown in Figure 15, for the lateral flow status, the simulation results and the experimental data [9] coincide well when we use the mesh and other simulation conditions introduced in the previous sections.

Besides the lateral flow status, the flow resistance is another important status which can affect the economy and safety of a reactor. As shown in Figure 16, a semiempirical formulation [15] is used to analyze the simulations on the flow resistance, and this formulation has a 10% uncertainty. Based on the analysis on the values and tendencies of the pressure distributions, it can be found that the simulations are in well agreement with the formulation, especially for the region with spacer grid (Region A) and the region downstream of the vertical location 0.2 m (Region B). However, the tendency has a transition region (0–0.2 m) between Region A and Region B in the simulations, and the formulation is with a sudden change on the tendency at the outlet of spacer grid. So the simulations may be more reasonable.

With these validations and analysis, our research approach can be used for the following researches.

3. Challenge Analysis

3.1. Great Need of Computing Resources

Based on the mesh independent analysis above for the feature regions (MV, DI, SP, and SC), the cell number can be estimated for a whole PWR. Table 4 gives the parameters of the nuclear power plant Qinshan I of China. Based on these parameters, the estimation of the CN for this reactor is listed in Table 5.

In this estimation, the CN for the inner region of spacer grid per channel is based on the sensitivity analysis for the dimple region (DI) and spring region (SP) in Table 3. In the sensitivity analysis, the vertical length of DI is half the vertical length of inner grid region, and the vertical length of SP is equal to the latter region. Meanwhile, one spring generally matches with two dimples. So the mean CN per layer can be equal to (103 × 2 + 229)/2 (≈218k/layer). In addition, the CN of one simple channel per meter is based on the sensitivity analysis for the SC region in Table 2. In the sensitivity analysis, the vertical length of SC is 90 mm, so the CN per meter is equal to 9.2 × 1000/90 (≈102k/m).

Then the total CN is around 70 billion, which is so large. Furthermore, the estimation is just based on the regions of SC, MV, and SG. The cell number for the regions between adjacent assemblies and the upper/downer plenum is not included. So the CN for a whole reactor must be larger than 70 billion.

In other estimations, such as the work of Conner et al. [16], the CN for a 5 × 5 rod bundle with one layer of grid spacer is around 20 million. Then the CN for a whole reactor must be more than 20 × 15 × 15/(5 × 5) × 8 × 121 million (174 billion). So the computing cost for a whole PWR is really extraordinary. In addition, the computing time is also very long even if the advanced HPC equipment and technology are used. This is really a great challenge for the CFD simulation of a whole PWR.

3.2. High Demand on Fine Mesh

Some researches have been focused on the coupling technology between CFD and MOC for a whole reactor, and they have achieved very good results [17, 18]. However, in the part of hydraulic analysis, the mesh density is not fine enough. This kind of mesh may affect the accuracy of the flow status, although it can save great computing resources.

As shown in Figure 17, (a) is the mesh type used in [18], and (b) is the mesh type tested by previous mesh independent analysis. Simulations using these two types on a 2 × 2 subchannels (a part of the experimental object [9] and the experimental conditions were also used in the simulations) were done and compared. As shown in Figure 18, two cross sections with vertical locations, 12.7 mm and 101.6 mm, were used for the comparison. The red line in Figure 17 is chosen for analysis. The velocity at the direction is used as the parameter for the comparison. In Figure 18, it can be found that the coarse mesh is not suitable. Though mesh type 1 is also very fine, it still create large error. So the mesh must be kept fine enough for the CFD simulation of a PWR.

3.3. Large Cost of Computing Time

As shown in Figure 19, similar regularity is suitable for each feature object (SC, MV, DI, and SP). This regularity is that the mean time per iteration step (MTPS) and the cell number (CN) have a good linear relationship which can be written as (1). In (1), and are two constant parameters for the linear relationship.

Regularity also exists between the CN and the iteration step number (ISN). As shown in Figure 20, it is also a good linear relationship for each research region. Equation (2) can be used to reflect this regularity, and and are two constant parameters for the linear relationship.

The total simulation time is decided by the ISN and MTPS. The product of two linear relationships makes a parabolic relationship. So the relationship between the total simulation time (TST) and CN is a parabolic relationship which can be seen from Figure 21. This relationship can be expressed as (3), and , , and are the constant parameters.

From the previous research, it was known that the fine mesh is in great need for the CFD simulation of a PWR. Then the high mesh density and the parabolic relationship decide that the total simulation time will be very long for the CFD simulation of a whole PWR.

4. Design of the Simulation Scheme

From the analysis above, it can be known that the computing resources and time are in great need for a CFD simulation of a whole PWR. Then the improvement of the simulation efficiency is very important because it is meaningless if the simulation is without an end or the simulation needs too much computing resources.

To improve the efficiency, one usual approach is to design an optimal mesh scheme with the minimum cell number (related to computing resources) and the least simulation time, as the mesh analysis in Section 2.3. However, both the resources and time are still great for the whole PWR after the application of optimal mesh scheme. So some researches were also done to develop other schemes to assist/obtain the further reduction of computing resources and time.

4.1. Quantitative Estimation Scheme for Efficiency

To improve the efficiency, the estimation of the efficiency has to be done firstly. However, few researches have been focused on this field. Though there is a lack of knowledge about the quantitative estimation of the efficiency, one point can be sure that the efficiency can be affected by the simulation conditions (the configuration of computing resources and the simulation schemes), as shown in Figure 22. In one aspect, the computing resources contains many factors, such as the number of computing nodes, the type and the clock speed of CPU, and the memory size. In the other aspect, there are also many simulation schemes which can affect the efficiency. For instance, the schemes can focus on the establishment of cell meshing, design of boundary layer, choice of turbulence model, discretization of numerical model, design of solution methods, and so on [19, 20].

All these factors above can affect the efficiency. So the efficiency should be considered by the overall performance of these factors. Furthermore, the research and development of the configuration of computing resources and design of simulation schemes are carried out in the international scope. The optimal one is also difficult to choose, because there is no quantative standard to estimate each international simulation condition (computing resources and simulation schemes).

To obtain an optimal simulation with high efficiency, a quantitative standard will be useful if it can give a combination assessment for the efficiency status of the computing resources and simulation schemes. To obtain a successful CFD simulation of a whole PWR, some work has been done.

4.1.1. Preestimation Method

As the regularity discussed in Section 3.3, the total simulation time can be estimated by some equations. Meanwhile, the main coefficients are all constant parameters because of the liner or parabolic relationship, and these constant parameters should be decided by the combined effect of the simulation conditions. The constant parameters , , , and can be obtained by two test simulations because of the linear relationship. Then from (3), , , and can be obtained by two test simulations, because is equal to , is equal to , and is equal to . In other words, the total simulation time can be estimated no matter how many the cell number is used, when two simulation tests have been carried out.

This regularity will be useful to give the quantitative reference. We can predict the nearly time before the simulation with great cell number. So this method is called preestimation. The details of the application of this method are listed as follows.

Firstly, based on the linear fitting, we can obtain the slops of linear relationship for ISN and CN and MTPS and CN. The slop between ISN and CN is set as , and the slop between MTPS and CN is set as . Then the product of and should be important for the efficiency when CN is great because is the coefficient of quadratic term. The need of cell number is great for a whole PWR; therefore is meaningful to decide the nearly total simulation time for the CFD simulation of a PWR.

Furthermore, this quantitative method just needs nearly two simulation tests. Meanwhile, just small cell number is enough for the simulation tests, because and are the slops for linear relationship, and a line can be obtained by two points. So the cost of test will be very small if the cell number is controlled in the tests. Certainly, more tests will be better, because the solution of the slop is by the linear fitting. However, too many tests will increase the simulation burden, and too many tests will not create obvious superiority because of linear relationship.

The steps of preestimation method can be generalized as follows: () two or several test simulations need to be done with the target configuration of computing resources and simulation schemes, and the cell number should be small in the test simulations; () linear fitting needs to be done on the simulation results to obtain ; () with the comparison of for the simulations using different simulation conditions, the judgment can be done.

4.1.2. Analysis Using Preestimation

To give the further explanation of preestimation method, the examples are given as follows.

Each computer is a combination of a set of computing resources, such as CPU type, clock speed, memory, and nodes quality for parallel computation. So two computers were used to respect different configurations. The details of the configurations are shown in Table 6. In addition, many simulation schemes also should be decided for a high efficient CFD simulation. So some schemes were considered as Table 7. Then several combinations of these simulation conditions can be set as shown in Table 8. These combinations respect different configurations, parallel nodes, turbulence models, and discrete methods.

With the application of preestimation approach, the ISN and MTPS were recorded, and the results are shown in Figure 23. Based on the good liner relationship, the slops ( and ) were computed, respectively. Then the product can be got as shown in Table 8. The smaller means the better efficiency when the cell number is great. So the quantitative standard is obtained, and the judgment can be easier for various combinations of simulation conditions. It can be used to clarify a more suitable one when the difference is small between two conditions or when there is no experience to choose the optimal one.

4.2. SSMS Scheme to Improve Efficiency

The mixing vane is one typical structure of PWR. It has great effect on the lateral flow. The spacer grid is another important structure. Its structure is complex for the existence of springs and dimples. Complex structure makes it difficult to establish the mesh. The cell number will be reduced obviously if the simplification can be done on the structure of the spacer grid. However, the influence of the simplification has to be kept little enough. So some research was also done to discuss the feasibility of the structure simplification.

4.2.1. Flow Characteristics due to Feature Structures

As shown in Figure 24, the research object is a part of AFA 2G spacer grid. (a) is with real structure, and (b) is simplified without springs and dimples. The second region is simple enough to use prismatic mesh. Prismatic mesh is better at the computing resources, speed, and simulation stability. On the contrary, tetrahedron and polyhedron can be easily used for the real structure of a spacer grid. However, the mesh cells will be great more, and the simulation time will also be longer. Moreover, the simulations may fail due to the poor stability. Figure 25 shows mesh for structure research.

The lateral flow status and the pressure profile were chosen for the comparison between the simulations using these two structure types. In the simulations, there are two flow channels with 300 mm length at the upstream and downstream of the spacer grid as shown in Figure 26. Meanwhile, the inlet temperature is 20°C, and the inlet pressure is 0.1 MPa.

In the comparison, six-cross section downstream of spacer grid has been chosen, and the vertical locations are 50 mm, 100 mm, 150 mm, 200 mm, 250 mm, and 300 mm, respectively. The zero vertical location is at the bottom of the spacer grid. The pathline for the data collection of each cross section is designed as the red line in Figure 27, and the length of this red line is 1 pitch.

As shown in Figure 28, the difference of the lateral flow status is not obvious, especially for the region downstream of 100 mm. However, as shown in Figure 29, different pressure profiles exist at the inner region of spacer grid. Fortunately, the pressure profiles at other locations upstream and downstream of spacer grid are not affected by the inner structure of spacer grid, which can be seen from the gradient of the pressure distributions. This regularity makes it possible to modify the simulation results at the outlet of the grid region, when the simplification scheme is used to the structure of spacer grid.

4.2.2. Simulation Using SSSM Scheme

From the analysis above, we can know the simple grid without spring and dimple also can reflect the nearly real flow status; however the modification has to be done for the pressure distribution in the inner of a spacer grid. This scheme is named as structure simplification and simulation modification (SSSM). As shown in Figure 30, the pressure profile using SSSM is nearly as same as the profile using the real grid. In this comparison, the pressure modification was only done to the outlet of the grid region using

Before the modification, further research is still needed to obtain the resistance characteristic of the grid region () using CFD simulation or experiment. The effect from the hydraulic process has been discussed above, and the liner distribution regularity can be used to modify the pressure drop.

However, the effect of thermal process can also affect the regularity of flow resistance. So the heat transfer was also considered for a 5 × 5 fuel assembly with AFA2G spacer grids. In the simulation, the inlet pressure is 15.5 MPa, and the heat flux was based on the vertical location. The value is equal to  W/m2, and it will be equal to 0 if the computation is less than 0. As shown in Figure 31, the liner regularity is still suitable for the inner regions of 5 × 5 spacer grids under various thermal hydraulic (T-H) conditions along the vertical length. To give a clear analysis, the various thermal statuses of these spacer grids were listed in Table 9, and it can be found that the difference of the thermal statuses is really obvious.

As shown in Figure 32, the enlarged figure for the pressure distributions along the vertical length of these SGs was also illustrated. It can be seen that the liner regularity is very good, and the main difference is the slop under each T-H condition.

The reason may be that the vertical length of the inner region of spacer grid is very short (around 30~40 mm), so the changes of temperature and other fluid properties could not be very large; then the thermal effect would not be obvious for the pressure profile of each layer of spacer grid.

As a result, due to the good liner regularity at the inner region of spacer grid under various T-H conditions, it is reasonable and feasible to modify the pressure profile for the inner region of the spacer grid using the simple structure.

As shown in Table 10, the comparisons of the cell number and simulation time were done for the single grid region with the real structure and simple structure using SSSM scheme. The SSSM scheme can reduce nearly half cell number, and the time of this scheme is near the one-fifth times of the original simulation time. Therefore, the computing resources and time can be reduced greatly when SSSM scheme is used.

5. Conclusions

CFD simulation is needed urgently in the development of NVR in Harbin Engineering University of China. For the successful application of the CFD simulation for a PWR in NVR, some work has been done and the conclusions were obtained as follows.

Firstly, based on the analysis on meshing methods and regularities of simulation process for the local feature regions and 5 × 5 rod bundle, suitable research methods (mesh type, cell size, numerical sets, etc.) were validated and quantitative challenges on computing resources and time were estimated. Large number of narrow coolant channels with complex structures makes the cell size must be fine enough. Then near 70 billion of mesh cells should be used in the CFD simulation for a whole research reactor. Both the relationship for the simulation time and the cell number and the demand on the fine mesh indicate the great need of computing time for the CFD simulation of a whole reactor. To realize a successful CFD simulation of a whole PWR, computing burden on resources and time must be reduced.

Secondarily, the quantitative regularities were used in the design of preestimation method which can be applied in the design and choice of high efficient simulation conditions (high efficient configuration of computing resources and simulation schemes). Furthermore, a scheme with structure simplification and simulation modification (SSSM) was designed, and the feasibility of this scheme was validated. Both these two approaches can assist/realize the reduction of computing burden. They will be useful and important to overcome the great challenges in the CFD simulation of a whole reactor.

Competing Interests

The authors declare that they have no competing interests.