Abstract
In order to analyze the steady state and transient behavior of CNA-II, several tasks were required. Methods and models were developed in several areas. HELIOS lattice models were developed and benchmarked against WIMS/MCNP5 results generated by NA-SA. Cross-sections for the coupled RELAP5/PARCS calculation were extracted from HELIOS within the GenPMAXS framework. The validation of both HELIOS and PARCS was performed primarily by comparisons to WIMS/PUMA and MCNP for idealized models. Special methods were developed to model the control rods and boron injection systems of CNA-II. The insertion of the rods is oblique, and a special routine was added to PARCS to treat this effect. CFD results combined with specialized mapping routines were used to model the boron injection system. In all cases there was good agreement in the results which provided confidence in the neutronics methods and modeling. A coupled code benchmark between U of M and U of Pisa is ongoing and results are still preliminary. Under a LOCA transient, the best estimate behavior of the core appears to be acceptable.
1. Introduction
The objective of the work here was to develop and benchmark methods and models for the coupled neutronics and thermal-hydraulics analysis of the CNA-II reactor. This ongoing effort was performed in conjunction with NA-SA and the University of Pisa. The following section provides an overview of the specific objectives of the work, and the subsequent sections provide the details of the work performed. A new computer model was developed for CNA-II based on the U.S. NRC version of PARCS. All modeling was performed for the detailed design of the CNA-II reactor core and its particular control rod configuration and orientation. An illustration of CNA-II is shown in Figure 1.
(a)
(b)
(c)
The neutronics analysis of Atucha-I and II involved the development of methods and models using HELIOS [1] and PARCS [2]. The overall approach in this work was to model the CNA-II core using both Cartesian and Hexagonal geometries. The Cartesian model followed the same guidelines established by NA-SA [3, 4], while the Hexagonal model was independently developed to be consistent with the triangular lattice of the CNA-II core. The neutronics task was to develop and benchmark HELIOS to NA-SA-generated WIMS [5] and MCNP5 lattice models for both geometries and to perform consistent comparisons for the models with published results [3] as well as agreed upon MCNP5 [6] core benchmarks with NA-SA. The second neutronics task was to develop a PARCS model of CNA-II using the HELIOS-homogenized few group cross sections. These neutronics models were then coupled to the NRC version of RELAP5 [7] in order to analyze systems level plant transients. A benchmark was performed between University of Michigan and the University of Pisa considering both steady state and transient scenarios [8, 9]. Current efforts are underway to complete a Chapter 15 review of CNA-II while model development and improvements are ongoing.
2. Cross-Section Modeling
The cross section modeling task involved a series of subtasks in order to generate the homogenized few group cross sections for the CNA-II reactor. The work completed included three main steps: creation and validation of lattice model, control rod models, and reflector models. The ultimate goal was to interface these cross sections into PARCS using the GenPMAXS [10] framework. The work performed is described in detail below. Included results follow previous work [11].
2.1. HELIOS
All cross sections models were generated with HELIOS version 1.10. Data for the cross sections are formulated into a 47-group adjusted library from ENDF-VIB data. Resonance treatment is accomplished through a subgroup method, and the transport equation is solved with current coupling and collision probability in arbitrary geometry. When running HELIOS, the cross sections were collapsed to two groups with an energy boundary of 0.625 eV. For resonance treatment, category four was used. This option separates the isotopes into four groups: U-238, heavy metals and fission products, natural Zirconium, and all other isotopes. Data transferred to PARCS include the macroscopic cross sections, discontinuity factors, and control rod branch information. All cross section data used in PARCS result from HELIOS calculations.
2.2. Validation of HELIOS Cross Sections
For the purpose of validation, NA-SA created a benchmark which considers a modified cold condition, referred to as idealized. MCNP5 models were employed in order to study the transport effects more closely. When considering the geometry, some assumptions were made with regards to the detail of the model. The main changes were the homogenization of the channel wall and the assumed increase in the fuel pin radius. These changes were shown to have little impact on results. Instead of placing the temperatures at the reference cold conditions, the values used in the idealized cases were a fixed default library value for MCNP5. This decreases the total required time for the calculation. Lastly, the materials were modified as well. All D2O densities correspond to reference cold conditions of sixty degrees Celsius, though the temperature of the heavy water was set to 20.45°C. The compositions of some materials were adjusted as well; standard Zircalloy four was replaced with natural zirconium, and the natural uranium fuel contained no impurities.
2.3. Cross Section Modeling: Fuel Assembly
The model used most frequently in these analyses was the hexagon, which represented the most natural lattice model for a triangular pitch core. It was used to extract all cross sections, except for the control rod incremental data. Figure 2 illustrates the meshing within the pin which was chosen to be fine enough to resolve the self-shielding.
Several idealized lattice models were studied for use in PARCS and for comparison with WIMS. Though the CNA-I lattice is used in this analysis, the CNA-II lattice is very similar in most respects. HELIOS and WIMS results are compared to the circular and hexagonal geometry results from MCNP5. MCNP5 calculations were performed by NA-SA. Relatively good agreement was observed between MCNP, HELIOS, and WIMS for circular and hexagonal models as shown in Table 1. Larger discrepancies were observed with the rectangular models, Table 2. The HELIOS hexagonal model was used to generate homogenized cross sections for the PARCS Hexagonal nodal model and the HELIOS rectangular model was used to generate homogenized cross sections for the PARCS Cartesian model.
2.4. Cross Section Modeling: Control Rods
In order to model the control rods, three approaches were used. Model one was a six-bundle model with a control rod in the center. Bundles are laid out in two rows of three into a square lattice array. The edit region was taken as one whole assembly with the control rod and guide tube. This provides a one-to-one ratio of fuel to control rods and was used as the fully rodded condition. This is shown in Figure 3, (left). Aside from approximating the lattice shape as square, the bundle to control rod ratio was low. The actual ratio is approximately twenty five, and any increase would serve to improve second-order effects.
The second model used eight bundles in a hexagonal lattice array. This not only improved the bundle-to-rod ratio but also represented the actual geometry more accurately. An illustration is shown in Figure 3, (center). As with the first model, this model allowed the extraction of a control rod incremental cross section for the fully rodded unit cell only. This precludes the possibility of accurately modeling a partially rodded node. However, this problem was addressed in the third approach.
The third model was created from a variation on model two. The advantages of a hexagonal assembly layout and an improved rod ratio were combined with a dynamic edit region to produce the desired branches. As with the earlier models, the edit region conserves the actual lattice area. In the case of the hexagon, retaining the shape was not feasible. An example of the model is shown in Figure 3, (right). The edit region always includes the bundle near the bottom, but the amount of control rod included can be adjusted. Figure 3, (right) shows a one half rodded case. With this model, an interpolation is possible for a half rodded cell which can represent the control rod tip part way inserted into a node. Other cases were created to bring the total for each rod type to four. This approach was used in order to create incrementals for partially rodded control rod configurations.
Control rod supercell models were used to extract incremental cross sections for the control rods. Two approaches were investigated for modeling the rodded cell, but results for the second approach are given for brevity; similar agreement was observed for the first approach. Table 3 gives results case for a steel rod in the two geometries, and the agreement is very good. With the heavy absorber, hafnium, DRAGON [12], and HELIOS have more difficulty with the worth of the rod. However, agreement is still acceptable in both cases. Table 4 shows results for the Hafnium cases.
For the PARCS models, the approach shown in Figure 3, (right) was used to extract the incremental cross sections. No direct comparisons were made with NA-SA, but the model is very similar to the eight bundle hexagonal model, Figure 3, (center).
2.5. Cross Section Modeling: Reflector
Three regions are of interest in terms of reflector cross sections; these are the upper axial, lower axial, and radial reflectors. Previous work has shown that because of the long neutron mean free path, heavy water systems are very sensitive to the treatment of the core boundary interface. In most cases, the axial reflectors are of less importance than the radial reflector. Consequently, the radial reflector was studied in detail. The vessel and tank walls were analyzed and included in some models. Albedos were employed to represent these outer structures when modeling was not feasible. After some study, the final model as shown in Figure 4 was a typical lattice with attached reflector. The reflector length was calculated as the core average radial reflector thickness.
As noted earlier, this model was improved with albedos. The model used for this purpose had an improved representation of the core boundary in which both the tank wall and downcomer were added to improve treatment of the transport effects. The approach was used only for the radial reflector.
Assembly discontinuity factors were also used to better represent the boundary interface. Implementation was performed by trial and error because of the irregular core boundary, and because PARCS only uses one ADF for all sides of the hexagon. Application of the ADF in some cells helped to improve the channel power distribution, but most cells had no discontinuity factors. Calculation of this ADF is shown in detail in the GenPMAXS manual. For both the upper and lower axial reflectors, the model shown in Figure 4 was used. Because the reflector is relatively thick, ADFs and albedos caused only minor changes in results.
3. Reactor Core Modeling
The PARCS modeling included three main subtasks: development and validation or both Cartesian and hexagonal core models, treatment and validation of the control rods, and analysis of the reflector modeling.
3.1. Introduction to the Models
The natural core configuration of CNA-II is a triangular lattice, but it can be represented as hexagons or rectangles. The NA-SA approach uses a rectangular representation of the core, and in order to compare with NA-SA, a PARCS model was developed that is very similar to this representation. A second PARCS model was created in hexagonal geometry. Each model has distinct advantages with the rectangular model providing a better treatment of the boundary and it also represents the control rods with a finer mesh. In terms of geometry, the hexagon is more natural, which allows for a simpler model and lower computing time.
3.2. Assembly Model Representation
Between the two PARCS models, the hexagon is the simplest and most consistent representation since one lattice model can be used for each fuel bundle. In this case, the triangle pitch is measured between the faces. When using the rectangular representation, the assemblies are modeled with an irregular repeating Cartesian array in which the bundle is divided into four quadrants. Each contains twenty five percent of both the moderator and fuel assembly. A representation is given in Figure 5. Each face of each quadrant used a unique ADF and the cell area is conserved in both cases.
3.3. Reflector Representation
Different models were employed for representation of the reflector in CNA-I and CNA-II. In the case of CNA-I, ADFs were not used, but they were used for CNA-II. Between the Cartesian and hexagonal core models, different approaches were used. Because of the small size of each quadrant in the Cartesian model, a much better boundary representation was possible. The result was a much smoother edge, which resembled more closely the circular tank wall. Figure 6(a) shows the nodalization used in the NA-SA code PUMA [13] which was also used in the PARCS Cartesian model, and Figure 6(b) shows the nodalization used in the PARCS.
(a)
(b)
The hexagonal model used a coarse treatment with a layer of reflector hexagons surrounding the core. Another complication was the ADFs in the hexagonal case. PARCS currently only accepts one ADF, all faces in the hexagon. Some reflector cells had ADFs, and some did not. Aside from the differences in the assembly representation, the objective in the reflector modeling was to conserve the radial reflector area on the core. The upper and lower reflectors were assigned the same cross sections as the radial cases, but without ADFs. Both PARCS models were also adjusted to better match the MCNP results provided by NA-SA.
3.4. Control Rod Representation
In the case of CNA-II, the control rods are inserted at a variety of angles and depths into the core. The oblique nature of the rods increases the complexity of the control rod model in PARCS. A module was developed in PARCS to represent the rods and also maintain the dynamic functionalities, such as rod ejection and bank movement. The contribution to each node is provided by the xy assembly, z plane, number of rods in the cell, and the fraction of each rod. Figure 7 illustrates the approach. The PARCS model allows the rod to be divided into an arbitrary number of nodes. In this case, the rod is able to contribute to all nodes it traverses. The rod is divided up and moved slightly in the xy direction to maintain the true xyz path, and each portion of the rod contributes to only one cell. This approach accurately predicts the rod worth and channel power very accurately.
3.5. CNA-I Idealized Core Modeling
PARCS results are compared in Table 5 for both the Hexagonal and Cartesian models to MCNP5 and the PUMA Cartesian results. Some differences between the two models are expected due to transport effects, boundary treatment, and the nodal solver; however the overall agreement is acceptable. PUMA more accurately predicts the channel powers, while PARCS more accurately predicts the core eigenvalue.
Comparisons of PARCS, PUMA, and MCNP5 were also performed for a model of CNA-II. When compared with CNA-I, the treatment of the boundary has a larger impact, which is the result of the core layout. The largest error in power predictions is in these locations. Again, the Cartesian model predicts a slightly lower eigenvalue.
Comparisons of PARCS, PUMA, and MCNP5 results were made with an idealized rodded condition in CNA-I. Assembly corrections (ADF/CHC) were also not used in PARCS. For these cases, six rods were inserted half way into the core. Two different cases were run. The first was with steel rods and the second with hafnium rods. In the case of the steel rods, both PUMA and PARCS predict the rod worth very accurately. However, errors in the channel powers are larger in PARCS which is likely due to the large hexagonal node into which the control rod incremental is homogenized.
When the rods are changed to hafnium, similar results are observed, though the rod worths are slightly under predicted. Because hafnium is such a strong absorber, diffusion theory sometimes has difficulty properly modeling the transport effects. As before, the channel powers are predicted slightly less accurately by PARCS compared with PUMA.
The previously discussed hexagonal core model of CNA-II was expanded to model the control rods and guide tubes. MCNP5 cases were run by NA-SA for the nominal and fully inserted rod positions. PUMA and PARCS were compared directly for the guide tube only case. The control rod model is very similar to the actual control rods in CNA-II which is somewhat different from the CNA-I idealized rod case. Both grey and black rods are axially heterogeneous and enter the core at a variety of angles and positions.
The significant asymmetry of the CNA-II core decreased the accuracy of both PARCS and PUMA with respect to channel power prediction. This discrepancy in comparison with the MCNP5 result increased with rod insertion and is likely due to transport effects. The control rod worth was also predicted less accurately in the fully inserted case.
3.6. CNA-I Startup Criticality Comparisons
A limited amount of commissioning data was available from the CNA-I startup activities. Because of the similarities between the reactors, these data were considered relevant and useful for the purposes of validation. Criticality was compared at two different cold states. Boron or a combination of boron and control rods was used to achieve criticality. Results for the cold cases are detailed in Tables 11 and 12.
There is good agreement for the boron and control rods configuration. Because of the complexity of the control rods, and the difficulty in Modeling, the case was directly calculation and found by interpolation from two similar cases. The smaller difference in the PARCS case results from the partially rodded cross section branches.
4. Boron Lance Geometry
The boron injection system in CNA-II consists of four lances as shown in Figure 9. Each of these lances contains one jet coming out of the end of the lance (head jet) and fifty smaller jets coming from the side of the lance (side jet) [14]. The head jet, which carries roughly 55 percent of the total volume flow through the lance, is of particular interest in this work because of its significant size.
4.1. Boron Jet Equations
Modeling the boron jets includes both theory and numerical data taken from reports by Siemens [15, 16], respectively. The resulting equations, as given by Drzewiecki et al. [17], are shown as follows.
Boron Jet Equations are as follows:
These equations give the spatial and temporal poison distribution of a jet. Some of the variables are fixed by the geometry and other design constraints. This leaves , , and as unknowns. Conservation of mass provides an extra equation which constrains as shown in [17]. The remaining equations for and were determined using CFD.
Mass Conservation Constraint is as follows:
The CFD work by Drzewiecki [18] focuses on the size and shape of the head jet. The CFD solution is computed using a two-dimensional axisymmetric domain using the porous media model in FLUENT [19]. From the CFD result the variable is extracted from the centerline profile and is extracted from the radial profile at several points. It should be noted that is found to be almost constant at different axial positions. Since the CFD work focuses only on the head jet, the side jets values are assumed to be the same as found in literature [16].
4.2. Constants in Boron Injection Equation
The constants for both the head and side jet are given by Drzewiecki et al. [17]. The coefficients of the head jet are recalculated to better match the trends on the short time scales instead of at the asymptotic limit of the jet. The new coefficients are provided in Table 13.
5. Boron Jet Mapping Methodology
In order to generate a discrete spatial and temporal boron distribution for the continuous boron equations introduced above, a series of mapping routines were developed.
5.1. Two-Dimensional Map Generation
The boron distribution equations for the axisymmetric jet are two-dimensional in space. To represent the jet in a discrete grid, the domain is broken into an by grid. The boron concentration equations must be integrated over the range Δ and Δ to conserve the boron mass from the jet. The two-dimensional grid is generated by dividing the domain into by equal segments. The length of the domain is determined by the penetration distance at the maximum time of interest.
Equation Discretization is as follows:
The radius of the domain is determined by the radial position where the jet meets some tolerance, . This radius is found by taking the derivative of the boron concentration equation in the diffusive region with respect to the axial direction and setting the derivative equal to zero as shown in (4).
Radial Domain Derivation is as follows:
This relationship between the axial location and the radial location can be inserted into the original equation to provide the radial position where is obtained as shown in (5).
Geometric Relationship is as follows:
Once the bounds of the domain are defined, it can be subdivided into equal axial and radial segments. In addition to the axial subdivisions, a subdivision is defined at the jet penetration distance which defines every time step. In the radial direction two methods were considered, equal radius spacing and equal volume spacing. After testing both cases it was found that the equal radius spacing conserves the total mass of boron in the system for the same number of radial segments. This is because the jet structure changes much more frequently near the center axis of the jet. The equal volume spacing does not have sufficient structure near the axis to capture all of these effects. The total boron concentration is conserved by integrating the boron distribution equations over each grid point as shown in (6).
Boron Concentration Discretization is as follows:
5.2. Three-Dimensional Grid Generation
Once the boron concentration for the two-dimensional grid has been calculated, the domain must be transformed into three dimensions. Since the jet is naturally symmetric in the azimuthal direction, a simple projection is adequate. The domain is divided into segments over 2 radians in the azimuthal direction. The resulting coordinate system is Cylindrical with every azimuthally symmetric position having the same boron concentration. For further analysis, the grid is projected into the Cartesian coordinates as shown in (7).
Cylindrical-to-Cartesian Coordinate Transformation is as follows:
Now that the three-dimensional domain has been generated, a series of transformations are used to position the jet in the correct direction and location in the reactor core. The three-dimensional domain in its current form has the origin of the jet placed at (0,0,0) with the axis of the jet pointing in the z-direction. The data for each grid point are arranged into a matrix of Cartesian points. For every grid point in the three-dimensional domain, a vector, , is defined to be the Cartesian coordinates of the grid point as shown in (8).
Transformation Vector is as follows:
Two rotations are needed to point the jet in the correct direction. For this task Given’s rotations are used. The first rotation is a clockwise rotation degrees around the y-axis, shown in (9). The second rotation is a counter-clockwise rotation degrees around the z-axis as shown in (10).
Rotation Matrix 1 is as follows:
Rotation Matrix 2 is as follows:
The last transformation is a coordinate transformation from the origin to the location of the jet entrance (Xo,Yo,Zo) as shown in (11). The new set of vectors, , define the spatial location in the reactor core of a single jet.
Transformation Matrix is as follows:
The three-dimensional domain defined by the vectors is a much finer grid structure than what is used in the PARCS model. The PARCS model consists of 20 axial levels and a hexagon that defines the region around each of the fuel assemblies in the x-y plane. To determine the average boron concentration in each of the PARCS nodes a mapping routine is developed to ensure that the total mass of boron is conserved. The mapping is done by assuming that the entire volume of the specific grid point can be assigned by the centroid of the grid point. At the edge of a PARCS node, this assumption does not hold, but in the fine mesh limit of the grid, the error introduced goes to zero. The averaging of the boron in each PARCS node is done by volume weighting the concentration in the grid as shown in (12).
Volume Weighting is as follows:
The volume weighting serves two purposes. It provides a way to determine the effect of different grid points that map into the same PARCS node. It also provides a way to quantify the mixing effect if the entire PARCS node is not enclosed inside the grid. This mapping routine is repeated for all jets and all time steps to generate an input file for PARCS that provides the spatial and temporal nodal boron concentrations. The overlapping of jets is treated by linearly superimposing the jets over each other.
5.3. Graphical Representation of Jets in Core
To validate the methodology above, it is important to examine shape of the boron concentration once it is mapped into the reactor. The first case was to examine a single lance and then to examine the case with multiple lances. This case was examined at different time steps to verify the propagation of the jets into the core. Lastly, all four lances will be shown at the same times as the two lance case. The single lance is shown in Figure 11. Two iso-surfaces are shown: 10 ppm in yellow and 500 ppm in blue. The head jet points down and to the left. The 50 side jets create a sheet of boron.
The two-lance case is used in the analysis of the critical break LOCA. Two lances are used because of a single failure of a safety system, and since the lances are connected to one of two pressurized boron tanks located on the north and south sides of the core, the worst single fault case is the failure of the pressurized tank which will cause two lances to fail and not to inject. Figure 12 shows the boron cloud at 0.25, 0.5, 0.75, and 1 second after the boron injects using an isometric view of the reactor.
The four-lance case is used to analyze the large break LOCA. Since this case is considered beyond design basis, the analysis allows for all safety systems to be available during the accident. Figure 13 shows the isometric view of the jets at 0.25, 0.5, 0.75, and 1.0 seconds after injection.
5.4. Computational Validation
As noted above, the free jet equation was used as the basis for this approach. CFD was used to correct the ideal equation for the actual configuration of the jet lance. Therefore, it is not applicable to validate the final model with CFD. The original design process validated the jet modeling with a neutronic calculation. KWU used their core simulator, IQSBOX, to calculate the boron jet reactivity worth [15]. A similar approach was used here. Along with comparison to the original modeling approach and results, a comparison was made with an analogous MCNP core model. For brevity, the MCNP results were excluded, but agreement was generally less than 300 pcm in boron cloud worth.
To validate the boron reactivity predicted by PARCS a series of comparisons were performed with Siemens data from KWU B311. Boron distribution maps as shown in Figure 14 were provided by KWU at 12 axial planes in the core at four times during the boron injection.
A PARCS model was built that assigns the appropriate boron concentration to the PARCS node that falls inside each of the KWU regions. The following conditions outlined in the report were used in the HELIOS lattice calculations to determine the cross sections for PARCS.(i)Average Burnup: 4.2 MWd/kgHM,(ii)Average Moderator Temperature: 180°C,(iii)Average Coolant Temperature: 296°C,(iv)Average Fuel Temperature: 637°C,(v)Average Coolant Density: 0.8003 g/cc,(vi)Average Moderator Density: 0.9903 g/cc,(vii)G10/S10 Insertion Depth: 53 cm,(viii)G20/G30 Insertion Depth: 371 cm.
HELIOS calculations were performed with six boron branches at 25, 100, 250, 500, 1000, and 3000 ppm. Between these branch cases PARCS performs a second-order interpolation to determine the exact cross-section incremental for the given boron concentration.
A PARCS calculation was first performed by simply mapping the KWU data at each of the four time points into the corresponding PARCS nodes. Table 14 shows the comparison of both reactivity and axial offset calculated in PARCS for all four time steps. There is reasonably good agreement in both the reactivity and the axial offset prediction. The calculation of axial offset is calculated as shown in the KWU report and reproduced in (13).
Axial Offset Equation is as follows:
A PARCS calculation was then performed with the boron distributions from the equations coming from the CFD result. The reactivity and axial offset predicted by PARCS is shown in Table 15. Again, the agreement between the two cases is reasonably good. The largest error in reactivity is 148 pcm at 3.75 seconds. The 3.75 second and 4.85 second cases differ from the KWU distribution because the equations conserve the mass that would flow out of the top of the upper reflector and average this into the upper reflector. The KWU distribution does not take this into account and so larger differences are expected at the longer times. The axial offset also has good agreement with the largest error at 4.85 seconds. This is likely due to the propagation distance of the jet since a jet entering from the top and traveling downward would cause the axial offset to be strongly negative and slowly reach zero when the boron jet is uniformly distributed along the axis of the reactor.
6. Coupled Benchmark
In order to study the hot full power behavior of the reactor, a coupled thermal-hydraulic and neutronics model was used. This benchmark is collaboration between NA-SA/PISA and ARN/UMICH. Modeling data provided by NA-SA are for the purposes of comparison only and results do not necessarily represent actual real plant behavior. Work on this benchmark and the overall reactor analysis is ongoing and results given here are preliminary. Further publications may allow for complete code-to-code comparisons.
The RELAP5 model was developed separately from an original input for CNA-II. This original KWU RELAP5 model (which contained 3 channels with 6 axial segments to represent the core region) was updated to contain five channels with twenty axial nodes to model the fuel. This RELAP5 model was coupled to the 451 channel, twenty axial node PARCS model to perform the CNA-II LOCA transients. The channel wall design in Atucha is represented by the RELAP5 pipe component which contains the coolant and receives heat from the fuel elements by convection and radiation. Heat structures are used to represent the fuel elements which are the basic component for solving the conduction and radiation heat transfer equations for delivering heat to the coolant. As noted above, only five coolant pipes were used in the RELAP5 model, but there are physically 451 channels in the core, each of which was represented explicitly in the PARCS neutronics model. A lumping approach was used in which each of the 451 neutronics channels was mapped into the five RELAP5 TH pipes. These five channels represent the minimum needed to properly model the orificing of the coolant channels. The pipe and heat structure have the same height as the actual core. Twenty axial segments were used in both the PARCS model and RELAP5 model; therefore the coupling was exact.
The mapping of the RELAP5 and PARCS nodalization for CNA-II is shown in Figure 15 for the radial model. The central 253 PARCS channels transfer their power to the first RELAP5 channel, and the coolant density and fuel temperature from this channel are then transferred back to PARCS. Therefore the central 253 channels in PARCS receive identical information. The same approach is used for the remaining RELAP5 channels in mapping the remaining PARCS nodes. In the axial direction, the data are passed explicitly between nodes, and no averaging is needed. Axial information is transferred in the same way as described for the radial mapping.
6.1. Steady State: CR in Critical Configuration with TH Feedback and HFP Conditions
This case helped to identify differences in the control rod cross section model and core flux solution with the control rods. The calculations were performed simulating the reactor at HFP conditions with control rods inserted in the critical configuration supplied by NA-SA. Total reactor power was imposed at 2160 MW. The main physical parameters are summarized in Table 13, and the radial axial power distributions are shown in Figures 16 and 17.
6.2. Transient: LOCA with 0.6-Second Boron Injection Delay Time
The coupled LOCA calculations are shown here for the best estimate case only. The initial condition for this case was calculated as indicated earlier. The break in the cold leg begins at 0.0 second and is fully opened at 0.00001 second. This initiates a full blow down with nearly complete emptying of the coolant inventory in the core. Combined with the negative Doppler feedback, the boron injection system activates 0.6 second after the break in the cold leg to shut down the core. The driving mechanism of this event is the coolant voiding because of the large positive coolant void reactivity coefficient in CNA-II. The coolant density in the core is given in Figure 18. The coolant density decreases dramatically and is less than 50% of nominal after only 0.5 second.
The components of the reactivity feedback are given in Figure 19. The CVR is driven by the coolant voiding and is nearly linear until about 0.75 second. The asymptotic CVR is about $2, which corresponds to earlier PARCS estimates. At 0.75 seconds, the coolant voiding slows as the inventory is depleted. After 0.5 second, the power increases from the superprompt reactivity causing the fuel temperature to rise. This can be seen as the negative Doppler reactivity is introduced. At 0.6 second, the boron injection is activated and reactivity begins to decrease. The power continues to increase until the total reactivity is less than $1 which occurs at about 0.7 second.
The core and peak nodal powers are shown in Figure 20. Given the superprompt reactivity insertion, the power increase is initially exponential. At about 0.7 second, when the reactivity drops below $1, the power begins to decrease. Given the similar trends of the core and nodal powers, the power change is primarily in magnitude.
7. Summary and Conclusions
The objective of the work here was to develop and to benchmark methods and models for the coupled neutronics and thermal-hydraulics analysis of the CNA-II reactor. A new computer model was developed for CNA-II based on the U.S. NRC versions of PARCS and RELAP5. All modeling was performed for the detailed design of the CNA-II reactor core and its unique control rod configuration and orientation. Validation was performed in collaboration with NA-SA and the University of Pisa using both lattice and core level models with the codes WIMS, PUMA, and MCNP. Acceptable agreement was observed in all benchmark cases.
Transient calculations were then performed with RELAP5/PARCS for a Loss of Coolant Accident (LOCA). Results show the expected behavior for a superprompt critical event driven by the large positive coolant void reactivity in CNA-II in which the power increases until the Doppler and boron injection are able to shut the reactor down. The computer models developed for this analysis appear to predict the LOCA event correctly and could be applied to other plant events and transient. However, further code and model development will be necessary to analyze other transients, particularly those in which the reactor control system is active.
Nomenclature
| ADF: | Assembly Discontinuity Factor |
| ARI: | All Rods In |
| ARN: | Autoridad Regulatoria Nuclear |
| ARO: | All Rods Out |
| CHC: | Cell Heterogeneity Correction |
| CFD: | Computational Fluid Dynamics |
| CAN: | Central Nuclear Atucha |
| CVR: | Coolant Void Reactivity |
| IQSBOX: | KWU core simulator |
| DRAGON: | 3D Neutron Transport Lattice Code |
| HELIOS: | 2D Neutron Transport Lattice Code |
| HFP: | Hot Full Power |
| HZP: | Hot Zero Power |
| KWU: | Kraftwerk Union |
| MCNP: | Monte Carlo N-Particle Code |
| NASA: | Nucleoelectrica Argentina S.A. |
| NRC: | Nuclear Regulatory Commission |
| PARCS: | Purdue Advanced Reactor Core Simulator |
| PHWR: | Pressurized Heavy Water Reactor |
| PISA: | University of Pisa |
| PUMA: | 3D Finite Difference Core Simulator |
| RELAP5: | Reactor Excursion & Leak Analysis Program (NRC) |
| TH: | Thermal Hydraulics |
| UMICH: | University of Michigan |
| WIMS: | Winfrith Improved Multigroup Scheme. |
Acknowledgments
The authors would like to thank Timothy Drzewiecki, Brian Mount, Professor Marin Bertodano, Autoridad Regulatoria Nuclear, Nucleoelectrica Argentina, and the University of Pisa for their collaboration on this project.