Nuclear Reactor Thermal Expansion Reactivity Effect Determination Using Finite Element Analysis Coupled with Monte Carlo Neutron Transport Analysis

Abstract

The energy released from the nuclear fission process drives thermal expansion and mechanical interactions in nuclear reactors. These phenomena cause changes in the neutron chain reaction which results in further changes in thermal expansion and mechanical interactions. Coupling finite element analysis with Monte Carlo neutron transport analysis provides a pathway to simulate the thermal expansion and mechanical interaction to determine a fundamental parameter, namely, thermal expansion temperature coefficient of reactivity. Knowing the coefficient value allows predictions of how a reactor will behave under transient conditions. Using the coupling of finite element analysis and Monte Carlo neutron transport analysis, the thermal expansion temperature coefficient of reactivity was determined for the Godiva-IV reactor (−2E−05 Δk/k/°C) and the Experimental Breeder Reactor-II (EBR-II) (−1.4E−03 $/°C). The Godiva-IV result is within 3% of the measured result. The thermal expansion and mechanical interactions within EBR-II are sufficiently complex that experimentally measuring the isolated coefficient of reactivity was not possible. However, the calculated result fits well with the integral EBR-II reactivity coefficient measurements. Coupling finite element analysis with Monte Carlo neutron transport analysis provides a powerful technique that gives reactor operators and designers greater confidence in reactor operating characteristics and safety margins.

1. Introduction

Nuclear reactors exhibit remarkably complicated behavior ultimately originating from the energy released through the nuclear fission process. The complicated behavior involves many phenomena including nuclear, thermal, and mechanical. Individually, these phenomena involve processes that are challenging to quantify, measure, and model. When interactions between these phenomena are considered, the quantification, measurement, and modeling challenges become daunting. This chapter describes finite element analysis (FEA) coupled with Monte Carlo analysis as a methodology for quantification of a particularly important nuclear parameter which is primarily influenced by thermal and mechanical phenomena present in nuclear reactors.

1.1 Background

The multiplication factor, k, is used to quantify the fission chain reaction in nuclear reactors. Numerous definitions exist for k, with each definition applying to a particular situation. A simple definition of k is that it represents the ratio of the number of fissions in one generation to the number of fissions in the preceding generation. Through this definition, one can see that if k is less than unity, the number of fissions declines over time, and if k is greater than unity, the number of fissions increases over time. A unique situation exists when k is exactly equal to one. In that case, the number of fissions remains constant over time and is referred to as critical.

A companion parameter to k is reactivity, ρ. Reactivity represents the deviation from the critical state, as shown in Eq. (1).

$$\rho =\frac{\left(k-1\right)}{k}$$E1

The decimal form of reactivity can be converted to units of $ by dividing the decimal value by the fraction of delayed neutrons resulting from the fission process. Delayed neutrons are those neutrons emitted during the decay of select radioactive fission products rather than being emitted at the moment of fission. For uranium-235, the delayed neutron fraction is 0.0065.

When operating a nuclear reactor, frequently, one is interested in knowing the change in reactivity resulting from various activities such as control rod movements. Other changes resulting from thermal and mechanical phenomena can produce reactivity changes. Frequently, these reactivity changes are quantified in terms of the change in reactor temperature. The result is known as a temperature coefficient of reactivity defined by Eq. (2).

$${\alpha }_T=\frac{\mathrm{\Delta }\rho }{\mathrm{\Delta }T}$$E2

The temperature coefficient of reactivity can be further subdivided into explicit subjects such as coolant temperature, fuel temperature, and even thermally driven reactor geometry changes.

From a reactor safety perspective, a negative temperature coefficient is indicative of inherently stability. If a reactor transient was initiated that results in a temperature increase, the resulting change in reactivity will necessarily be negative, which means the multiplication factor will be reduced. Eventually, the temperature increase will produce a sufficient reduction in k such that the reactor will shut down. Contrarily, a positive reactivity temperature coefficient is indicative of inherent instability. With a positive coefficient, a transient resulting in a reactor temperature increase will result in a positive reactivity change and a resulting increase in the multiplication factor. The increased multiplication factor will be accompanied by an increase in the number of fissions and resulting heat release and corresponding temperature increase which will subsequently produce an additional positive change in reactivity. The reactor will continue on this path until it is acted upon by a more dominate negative action or the reactor will ultimately be damaged or even destroyed.

Reactivity coefficients can be determined for numerous phenomena. For example, reactivity coefficients can be established for changes in reactor power. Thus, as the reactor power is increased, the reactivity change needed to compensate for the power change can be identified. Another interesting phenomenon that has a significant reactivity effect centers on bubble or void formation as the result of coolant boiling. Reactor designers must pay particular attention to the reactivity effect associated with coolant bubble formation because it can have a significant safety impact.

In the case of reactors that use water as a coolant, the water has a significant effect on the overall neutron energy spectrum in the reactor. Neutrons tend to be born at high energies on the order of several million electron volts. As the neutrons collide with various nuclei in the reactor, they tend to lose energy in a process called moderation. As the neutrons lose energy, they become more likely to be absorbed in uranium-235, which can then fission and release additional neutrons. Similar to the three different regimes for k, there are three regimes for moderation: under-moderation, optimum-moderation, and over-moderation. If a reactor is designed with under-moderation, the loss of coolant through bubble formation will result in a reactivity decrease because fewer neutrons will be slowed to energies where they are more likely to be absorbed in uranium-235. In the case of a reactor designed with over-moderation, the formation of bubbles will tend to result in a reactivity increase because more neutrons will be slowed to the point where they will be absorbed by uranium-235 causing an increase in the number of fissions.

The most dramatic and tragic demonstration of positive reactivity due to bubble formation was seen in the 1986 Chernobyl accident. When operated at low power, the Chernobyl reactor had a positive void reactivity coefficient. Thus, if the reactor coolant began to boil, the bubbles created by the coolant boiling led to a positive reactivity change thereby driving an increase in the multiplication factor and a corresponding increase in the number of fissions occurring in the reactor. The heat released from the additional fissions led to additional coolant boiling which drove a very rapid power increase and subsequent steam explosion and reactor destruction.

Reactivity coefficients tied to geometry changes are of interest in certain situations because they are typically fast acting and can have important safety implications. In many cases, thermally driven geometry changes are coupled with resulting mechanical interactions that severely complicate quantification and modeling approaches. While the change in geometry causes the reactivity change, typically temperature is used to quantify the reactivity coefficient since it is a change in temperature that causes thermal expansion and mechanical interaction. Thus, a reactor may have a thermal expansion temperature coefficient of reactivity or even more specifically a thermal expansion/mechanical interaction temperature coefficient of reactivity. The thermal expansion temperature coefficient of reactivity in two reactors is described below followed by a demonstration of using finite element analysis to model the thermally driven geometric changes followed by use of a Monte Carlo simulation to determine the corresponding multiplication factor value.

1.2 Godiva-IV and Experimental Breeder Reactor-II

Two reactors serve as the test bed for evaluating the analysis approach described in this chapter. One reactor a uses comparatively simple design and the other is significantly more complicated. The simple reactor design is called Godiva-IV. The Godiva-IV reactor, see Figure 1, is unique in that it is designed to provide a burst of neutrons rather than being designed for extended steady state power production. The Godiva-IV reactor is very compact with a simple cylindrical shape with a 178 mm diameter and a 156 mm height. The reactor design uses a solid construction of approximately 66 kg of 93% enriched uranium alloyed with 1.5 wt.% molybdenum. No active cooling arrangement is used. The reactor construction is somewhat more complicated than a monolithic cylinder of enriched uranium. The Godiva-IV reactor uses six ostensibly equal rings. Three stacked cylinders of differing heights are located within the six rings to complete the overall cylindrical shape. Three large C-clamps are attached to the outer radius for the reactor to restrain the fuel movement during burst operations.

Figure 1.

Godiva-IV reactor [1].

When the Godiva-IV reactor is operated, a large power pulse occurs and heat from the fission process is deposited in the uranium alloy. The heat causes a temperature increase and subsequent thermal expansion. As the individual components of the reactor expand, they mechanically interact. As the reactor components expand, neutron leakage from the reactor increases which leads to a decrease in the multiplication factor and subsequent termination of the reactor power pulse. Thus, Godiva-IV has a negative reactivity temperature coefficient. That is, as the reactor temperature increases, the resulting reactivity change is negative which provides an inherent shutdown mechanism.

The other reactor used to evaluate the analysis approach described in this chapter is the Experimental Breeder Reactor-II (EBR-II) [2]. The EBR-II design is significantly more complicated than Godiva-IV, see Figure 2. EBR-II uses liquid sodium metal as the coolant. The fuel is 67% enriched uranium metal alloyed with a collection of various metals totaling 5 wt.%. The fuel is formed into individual 3.3-mm diameter pins along with stainless steel cladding. The fuel portion of the pins is 343 mm long while the cladding portion is 638 mm long. The additional length of the cladding allows for the containment of fission product gasses. A collection of 91 fuel pins are arranged into a hexagonal configuration which is commonly referred to an assembly. The 91 fuel pins in each assembly are contained within a stainless-steel hexagonal duct. The EBR-II reactor core consists of an arrangement of 637 assemblies. The core is fundamentally divided into two regions, a driver region containing the fissile material, and a blanket region containing depleted uranium. Within the driver region there are approximately 100 assemblies including control rods, experimental assemblies, stainless steel dummy assemblies, stainless steel reflector assemblies, and assemblies that use reduced fuel content. Surrounding the driver region is a collection of approximately 500 assemblies constructed of depleted uranium. The depleted uranium assemblies absorb neutrons that leak for the driver region to transmute depleted uranium to plutonium to breed new reactor fuel.

Figure 2.

EBR-II.

As EBR-II ascends to its operating power, heat from the fission process causes the fuel pins, hexagonal ducts, and all other components of the reactor to expand due to the temperature increase. The components undergo a complicated process involving thermal expansion and mechanical interaction. While the Godiva-IV thermal expansion process which leads to comparatively simple thermal expansion and mechanical interaction, the thermal expansion and mechanical interactions in EBR-II are significantly more complicated and must be subdivided into different areas. One area that is comparatively simple to understand and evaluate is the spacing of the assemblies in the hexagonal arrangement. As the reactor temperature increases during the reactor assent to power, the grid plate that holds the fuel assemblies thermally expands and the spacing between the fuel assemblies increases which results in increased neutron leakage and a decrease in the multiplication factor. A much more complicated process involves the thermal expansion and mechanical interaction of the stainless-steel hexagonal assembly ducts. Measuring and calculating the reactivity effect of the hexagonal duct thermal expansion and mechanical interaction is particularly challenging. The analysis method described in this chapter is used to evaluate the reactivity coefficient associated with the thermal expansion driven spacing of the assemblies along with the much more complicated reactivity coefficient associated with the thermal expansion and mechanical interaction of the fuel assembly hexagonal ducts.

2. Method

One of the difficulties with quantifying a geometric temperature coefficient is the complexity of the thermal expansion. Thermal expansion coefficients are nominally nonlinear, leading to different rates of expansion depending on how hot the geometry is in a given location. This leads to nonlinear thermal expansion. This is an important concept to understand because it drives the necessity for using more complex structural analysis techniques than first principles expansion. This is especially true when geometric expansion is mechanically restrained by other expanding materials.

The key to successfully quantifying a thermal expansion derived temperature coefficient is not the calculation of the coefficient itself, but more the mechanical model that is used to derive the geometry changes. To that end, finite element analysis is used to provide a high fidelity mechanical input into the Monte Carlo simulation [3]. Figure 3 shows the generalized process for quantifying the temperature coefficient.

Figure 3.

General process flow.

2.1 Finite element analysis

Regardless of the source of the geometry information, whether an existing CAD model is defeatured or built from scratch, a simplified CAD geometry should be generated. The simplified geometry should contain enough information such that any complex expansion is captured, but simple enough to reduce the overall element count. A common example is removal of bolts and generally any small features from large geometries. FEA models in general run the risk of being too-large-to-compute without using resources unavailable to the typical engineer. Keeping total element count to a minimum is a driving factor when constructing an FEA model.

2.1.1 Mesh size

Exceeding 10 million nodes in a given model almost certainly means the model cannot be executed on a workstation in any reasonable amount of time. The reason for this is the sheer size of the data generated. A 10 million-node model requires 10 million positions (x, y, z), temperatures, and displacements (x, y, z) for one solution step. A double precision number requires eight bytes leading to each node requiring seven-, eight-byte numbers (56 bytes). 56 bytes per node applied to a 10 million-node mesh leads to 560 MB to store just the results of the model for 1 timestep or substep, not including the other required parameters for the solution, heat flux, power, boundary conditions, structural support, etc. Assuming the model requires several hundred timesteps and thousands of substeps, the total data requirement becomes multiple terabytes that needs to be loaded into memory. At the time of writing, several terabytes of memory was only available on very high-end workstations and was problematic for large HPC machines due to the memory allocation per CPU.

Given the difficulties noted above, simplifying the CAD model to reduce node count is critically important. The limited memory should prioritize the expansion effects not necessarily geometric fidelity. Even with these reductions, the model might take weeks of runtime to complete.

2.1.2 Boundary conditions

Boundary conditions are required as inputs into the FEA models. The boundary conditions are what simulate the reactor state that causes the structural change. They are divided into two types, thermal and structural.

The thermal boundary conditions consist of heat transfer coefficients and a thermal load. The thermal load will nominally be the fission source distribution based upon total power output. Determining the heat transfer coefficients can be done either using calculated values, measured values or a combination of both. Applying this method to a real system generally necessitates both. For example, a coolant flow rate is known for the entire reactor but is unknown on an assembly-by-assembly basis. The main goal of determining boundary conditions is to take what is known, and calculate what is needed for the FEA simulation. A similar problem exists with the thermal boundary conditions as with the mesh size, creating individual boundary conditions for every assembly and coolant channel can lead thousands of boundary conditions. It will be up to the user to determine if a thermal hydraulic simulation is required to calculate the heat transfer coefficients, or if hand calculations can suffice.

The structural boundary conditions in many ways are easier and less numerous. This is because the boundary conditions are only needed to simulate real restrains and conditions to aid in FEA convergence. An example of a real restraint is supporting the body in space such that it does not fall forever due to gravity. This is known as rigid body movement.

Convergence aids are sometimes required such that the simulation will converge. Convergence aids can be limiting body movement to a reasonable amount or declaring that a body cannot move during a particular substep. The primary pitfall with structural boundary conditions will be to overconstrain the model such that whatever subtle structural effect that drives the thermal coefficient is not nulled by the boundary conditions.

The boundary conditions need to be expansive enough to both simulate the reactor state and give enough information to allow the model to come to a solution. Additionally, the boundary conditions need to not overly constrain the model such that multiple solutions exist for a given state.

2.1.3 Timesteps

Selection of timestepping in FEA is another balance of model fidelity and analysis time. Finer timestepping leads to more information captured for a given effect, but can lead to longer computation times. For example, if the model has 100 steps that need a week of computation time to simulate 100 seconds of reactor time, then how are those steps distributed such that the necessary effects are captured? Are 90 steps used over 2 seconds of reactor time enough to capture the effect, with 10 steps used for 98 seconds of reactor time? Answering that question is highly problem dependent and takes multiple iterations to refine. Nominally, high fidelity stepping is required when the model is undergoing rapid geometry changes. For example, a pulse reactor can go from room temperature to several hundred degrees in a matter of milliseconds. During that time, many solution steps are required since the model is changing significantly between each solution step, whereas during cool-down of that same system, the model is changing slowly as conduction takes place.

2.2 Monte Carlo neutron transport

One of the problems stated in a previous section was the lack of modeling fidelity to capture subtleties in geometric changes. While the solution for the mechanical input utilizes unstructured mesh to define the geometry, Monte Carlo tools use more simplified geometry. Some Monte Carlo codes can take an unstructured mesh as an input to create their geometries, but nominally, some amount of translation will be needed to take the unstructured mesh and import it into a Monte Carlo code. Details that were required for the FEA may not be needed for the neutron transport.

Determining the multiplication factor using the Monte Carlo method requires four fundamental items, first, explicit geometric and material descriptions, second, detailed material nuclear property data, third, mathematical processes for sampling nuclear data, and fourth, a method for generating a string of numbers that satisfy rigorous tests for randomness. Using the four fundamental items, a simulation can be conducted for an individual neutron. The simulation begins by selecting the initial birth location for a neutron. The location is initially specified by the analyst but is later selected based on the location of fission locations from a prior generation. Once the birth location is known, the neutron energy can be randomly selected using numbers from the string of numbers that satisfy the randomness criteria and the mathematical distribution of possible neutron energies. The neutron direction can then be determined by randomly selecting an azimuthal and polar direction in the case of an isotropic direction assumption. With the neutron energy and direction being randomly selected, the distance the neutron travels before colliding with a nucleus can be randomly determined based on nuclear data associated with the probability of interaction, commonly referred to as a cross-section. Once the distance traveled is known, the type of interaction (e.g., scattering, absorption, and fission) can be randomly determined based on the ratio of cross-sections for the various interactions. It is also possible that the selected distance may result in the neutron leaking from the system and a new neutron must be generated. In the case of a scattering event, a new direction and neutron energy, based on collision mechanics and nuclear data, are randomly selected, and a new path length is selected. In the case of absorption, tracking of that neutron is discontinued, and a new neutron must be generated. If a fission event is selected, the location is recorded, and the number of neutrons produced by the fission event is randomly selected [4].

To accelerate the process, the analog simulation described above is modified with mathematically justified non-analog variance reduction techniques. These non-analog variance reduction techniques are selected based on the trade-off between computational time and a reduction in the statistical uncertainty of the result. For example, a process referred to as survival biasing is commonly applied where neutrons that are selected for absorption are only “partially” absorbed thereby allowing the remaining portion of the neutron to continue being tracked. The general idea is that it is more efficient to track a portion of a neutron than to track a neutron for an extended history only to have it eliminated in a meaningless reaction. As long as the variance reduction technique maintains a fair game, it can be used. The process, using analog and non-analog techniques, is repeated a great number of times, and then parameters of interest such as the multiplication faction can be inferred. The multiplication factor is inferred by the ratio of neutrons generated in one generation to the number of neutrons generated in the prior generation. A simulation can require more than 10¹² random numbers, billions of neutrons, and thousands of generations to obtain sufficient statistical confidence in the result. For the work discussed in this chapter, the Monte Carlo code MCNP® was used for the multiplication factor calculations [5].

2.2.1 FEA interface to neutron transport code

Nominally, the results generated from structural FEA will be nodes that have been displaced in space. These nodes will need to be translated to the Monte Carlo neutron transport geometry definition. Even using unstructured mesh as an input will require some modification and additional geometry because not all of those parts were required for the FEA, but might need to be in place for the neutron transport. As stated previously, the mesh size will need to be kept to a minimum, hence some amount of defeaturing took place. Some of those features will need to be restored for the Monte Carlo neutron transport such that the particle population is simulated correctly. An example of geometry that would be removed for the FEA, would be explicit detail of the fuel pins in a nuclear fuel assembly. The FEA nominally would not require such detail, opting instead for bulk heating of the channel. Adding geometry after the FEA presents the problem of fitting un-deformed geometry into deformed spaces. Resolving this issue requires the use of custom computer codes to perform the geometry translation and geometry checking to make sure there are no overlaps. These codes nominally are custom to the particular reactor or nuclear system. In the following section “A Complex Example,” EBR-II required a code called MICKA to perform the Monte Carlo geometry construction and translation [6].

2.2.2 Quasi-static snapshots

Nominally, thermal expansion reactivity coefficients that are nonlinear need to be analyzed through the whole power range of the reactor to capture all of the possible geometry states. The FEA will calculate the geometry displacement through the power-band and then export the data. That data will be exported and translated as a series of geometry snapshots with each snapshot representing the deformation of the reactor at particular power level.

2.2.3 Temperature coefficient calculations

The quasi-static snapshots are individually analyzed for their respective multiplication factor. Each individual snapshot is not intrinsically valuable because temperature coefficients in general represent a trend over a particular range. The important value to calculate over these snapshots is the change in multiplication factor, reactivity. Each reactivity point associated with a particular bulk temperature of the reactor is plotted. The slope of the linear fit of those reactivity points is the temperature coefficient. For nonlinear temperature coefficients, a set of linear fits are derived where each coefficient has an associated temperature band where the coefficient holds true. Before demonstrating the fitting process, change in multiplication from a noncritical state needs to be discussed.

Change in multiplication from critical was shown in Eq. (2). While change from critical does have a use, most real multiplying systems are never perfectly critical (k = 1), they are nominally slightly super or subcritical. This holds true for analyzed systems as well. Monte Carlo methods by definition have uncertainty associated with whatever values are calculated and rarely yield k = 1. A more common occurrence is k = 0.998 ± 0.004. The previous value would be considered critical in any real sense; however, from a calculation perspective it is noncritical. A modification to Eq. (2) is needed. Eq. (3) can be used for change in the reactivity.

$$\mathrm{\Delta }\rho =\frac{\left(k_2-k_1\right)}{k_2{k}_1}$$E3

With an understanding of change in reactivity, a linear temperature coefficient can be determined from the Monte Carlo analysis. A linear regression is applied to the temperature dependent reactivity. This yields Eq. (4).

$$\mathrm{y}=mx+\mathrm{b}$$E4

The slope of the previous equation is the temperature coefficient. If the particular coefficient is nonlinear, multiple regressions will be required. The coefficient can be expressed as a nonlinear equation; however, temperature coefficients are traditionally expressed as linear quantities over temperature ranges.

3. A simple example

The following sections demonstrate how finite element analysis can be applied to nuclear systems to calculate extremely complex phenomena. The tools used in these works were ANSYS, for the finite element analysis, and MCNP® was used for the neutron transport [5, 7].

3.1 Godiva-IV

Confirmation of the methodology described above begins with a relatively simple geometry reactor. While certainly not a homogeneous single component bare cylindrical reactor, the Godiva-IV reactor provides an excellent case to apply the methodology described above and compare the modeling results to measured results. In particular, the Godiva-IV reactor temperature coefficient of reactivity is dominated by thermal expansion with limited mechanical interaction effects. Furthermore, the Godiva-IV reactor has been thoroughly characterized and detailed descriptions are available to allow construction of both the FEA model as well as the Monte Carlo model. Finally, detailed temperature measurements and the corresponding reactivity have been recorded which allow for comparison with the modeling results.

3.1.1 Thermal analysis

The thermal analysis required several boundary conditions as inputs into the model. The most important was the temperature data taken from an experiment on the Godiva-IV. The temperature data provided the pulse shape and total magnitude of the temperature rise. It also conveyed the time dependency of the FEA model. Experimental measurements are important to creating accurate models. They provide a more accurate input, depending on the quality of the measurement, than necessarily calculating and input from first principles. The experiment input data were for a 1.029$ reactivity insertion with a 68°C temperature rise over 300 seconds.

The second type of boundary conditions applied were the heat transfer coefficients, primarily the convection coefficients. These were hand calculated from heat and mass transfer equations. Hand-calculating these coefficients is normally required because these values are not normally measured for these facilities. For Godiva-IV specifically, capturing the thermal expansion temperature coefficient means modeling the thermal conduction paths as well as the convection into the room. The temperature differential of all of the components as they heat up and subsequently cool due to convection to the room is the primary driver to the complexity of the expansion. Figure 4 shows the thermal FEA results and the temperature differentials. For more information specifically on the boundary conditions, see the reference [8].

Figure 4.

GODIVA-IV thermal analysis results.

Given the rapid structural response of the Godiva-IV pulse, the analysis type chosen was transient. The solution steps were focused on the pulse. Of the 400 solution steps, 300 steps surrounded the pulse that consisted of 10s, with the other 100 steps covering 290 seconds. The reason for this was that the model was rapidly changing during the initial nuclear heating, while the rest of the steps only contained relatively slow thermal conduction. After the temperature analysis model completed, the temperature data were exported to the structural analysis model.

3.1.2 Structural analysis

As stated previously, the structural boundary conditions are less numerous but can be more difficult to determine. For Godiva-IV, the primary structural boundary conditions were, support of the safety block, control rods, and providing a fixed support for the back side of the clamps. These boundary conditions were more straightforward than for typical models, weak springs were not required, and fixed supports were the only type of boundary conditions that were necessary. Figure 5 shows exaggerated displacement at the end of the temperature input data. The exaggeration was required because the structural displacement on average was 0.2 mm and imperceptible to the human eye.

Figure 5.

Exaggerated structural displacement.

The structural analysis had similar timestepping as the thermal analysis. The data generated from the FEA were a set of averaged displacements on particular surfaces. These surfaces surrounded the curved faces of the fuel rings. The fuel rings were the focus because only the fuel movement and expansion matters to the neutronics of the reactor.

Figure 6.

Godiva-IV temperature coefficient result.

3.1.3 Neutron transport

The exported data were applied to the neutron transport model where a series of models were created, each with a different set of displacements per an average temperature. These results are shown in Figure 6. The slope of the linear regression is −2E−05 Δk/k/°C. The comparison to the measured value is shown in Table 1. This is the temperature coefficient of the Godiva-IV reactor. The results demonstrate that a coupling method of FEA and Monte Carlo Neutron Transport has to be potential to accurately predict the temperature coefficient.

Figure 7.

EBR-II simplified simulation model.

Source	Temperature coefficient (Δk/k/°C)
Measured	−1.95E−05
Calculated	−2E−05

Table 1.

Measured and calculated Godiva-IV temperature coefficient.

4. A complex example

4.1 EBR-II

With a comparatively simple application providing excellent comparison results, a more challenging application is warranted. As noted above, the EBR-II design includes numerous fuel assemblies, molten sodium coolant, and a complicated thermal expansion and mechanical interaction process. Detailed characterization of the reactor components and materials along with measurements of control rod critical positions and corresponding bulk coolant temperatures are available [9]. These measured data allow confirmation of the methodology for certain aspects of the reactivity coefficients present in EBR-II such as the thermal expansion of the reactor grid plate. Extrapolation of the methodology can then occur for the more complicated thermal expansion and mechanical interaction of the assembly hexagonal flow ducts.

While no reliable method of measuring the reactivity coefficient associated with the hexagonal duct expansion and mechanical interactions is known to exist, the methodology described here can be applied and a reliable estimate of the reactivity coefficients can be obtained.

4.1.1 Thermal analysis

EBR-II required more extensive thermal boundary conditions than Godiva-IV which was considered the simple system because the heating was simple conduction through the materials and the ultimate heat sink was convection into the room air. EBR-II was more complicated because there was forced convection using liquid sodium that flowed over the fuel elements. The ultimate heat sink was a series of heat exchangers that cooled the sodium. Heating was also not symmetric from assembly to assembly. Modeling this complex behavior would require a complex thermal-hydraulic model to simulate the various coolant channels. Creating this model would have substantially complicated the thermal FEA analysis and additional would require input information that was not measured at EBR-II. Instead, a simple cooling model was developed for each assembly. The cooling model stated that the sodium inside of the duct entered the bottom of the coolant channel as cold sodium, and over the part of the channel where the fuel was located, the sodium was heated such that the outlet temperature match measurements taken at the EBR-II. Each fuel assembly type had a different cooling profile. Figure 7 shows the different assembly types in the EBR-II FEA model.

Figure 8.

EBR-II exaggerated structural displacement.

Simplifying the coolant channels in this manner was sufficient because previous work done on the EBR-II suggested that duct-bowing was entirely driven by the temperature profile of the duct material and not by the internal structures. Thus, only the duct needed to be heated correctly.

The power input for EBR-II was derived from a linear interpolation of the ascent to power. All of the heat generation inputs were linearly scaled over timesteps. The timing did not match the real ascent to power, but that was not necessary since the model would be in thermal equilibrium for each calculated step. The more important aspect was that the thermal model would be a series of steps, each step corresponding to a different power level.

4.1.2 Structural analysis

The structural analysis required a simple boundary condition to hold the model in place, as well as a boundary condition to fix the center duct. Fixing the center duct meant that it was not allowed to thermally expand and was considered a rigid body. This was necessary to achieve convergence. Without fixing the center duct, the model could not resolve the contact overlap that existed between the ducts on the first solution step.

The structural FEA required significantly more time to solve than Godiva-IV due to the sheer size of the model (~5 million nodes) and the complexity of the thermal expansion. Figure 8 shows an exaggerated displacement of the ducts. The southeast quadrant shows how the differences in the assembly, types, and powers can impact thermal expansion. Additionally, it demonstrates why FEA was necessary to capture all of the geometric detail of the duct-bowing temperature coefficient.

4.1.3 Neutron transport

Similar to the Godiva-IV model, the displacement data were exported out of ANSYS and imported into a series of MCNP® models [7]. The major difference was in the translation method. The Godiva-IV translation was averaging nodal thermal expansion and manually applying the change in radii and heights to the MCNP® input files. That approach was prohibitive for EBR-II because the resulting data exceeded 1 TB. A custom code called MCNP® Input Card and KCODE Architect (MICKA) was written to perform the node translation and MCNP® input construction. The MCNP® model for the EBR-II was itself expansive and required special data handling. More inf0rmation can be found in the reference [6].

One additional difficulty with the translation was that MCNP® cannot model a bowed-duct, only a straight hexagonal duct. To overcome this geometry limitation, the straight duct was divided into axial sections. Each axial section was moved in space to approximate a bowed-duct. Figure 9 shows an example of the axial slices in an assembly.

Figure 9.

Axial sections to simulate a bowed-duct in MCNP®.

After the translation of the nodal data from the FEA to MCNP®, the analysis process was similar to that of Godiva-IV. A series of snapshots at various bulk temperatures were taken and a linear regression was performed to calculate the slope of the points. Figure 10 shows the results of the reactivity change per degree. The coefficient was calculated to be −1.4E−03 $/°C. While the data had a clear linear trend, some nonlinearity existed in sets of data points at lower bulk temperatures. This was consistent with historical measurements at EBR-II where lower powers exhibited a nonlinear trend in the reactivity change.

Figure 10.

Temperature coefficient results for duct-bowing coefficient.

5. Conclusions

The energy released from the nuclear fission process drives complicated thermal expansion and mechanical interactions in nuclear reactors. These expansions and interactions subsequently cause changes in the neutron chain reaction balance within a reactor which results in further changes in thermal expansion and mechanical interactions. Measurement of these coupled phenomena occurring within a reactor has proven to be elusive. However, coupling finite element analysis with Monte Carlo neutron transport analysis provides a pathway to simulate the thermal expansion and mechanical interaction driven by the energy released in the neutron-induced fission process and then to subsequently determine fundamental nuclear parameters, namely, thermal expansion temperature coefficient of reactivity.

There are important safety implications associated with the thermal expansion temperature coefficient of reactivity and its relation to other temperature coefficients of reactivity. Knowing both the sign and magnitude of individual coefficients allows reactor designers to predict how a reactor will behave under transient conditions.

Using the coupling of finite element analysis and Monte Carlo neutron transport analysis, the thermal expansion temperature coefficient of reactivity was determined for the Godiva-IV reactor and found to be within 3% of the experimentally measured value.

The coupling technique was also used to determine the thermal expansion temperature coefficient of reactivity for EBR-II. The thermal expansion and mechanical interactions within EBR-II are sufficiently complex that experimentally measuring the isolated thermal expansion temperature coefficient of reactivity was not possible. However, using the coupling technique, a calculated value of −1.4E−03 $/°C was determined for the thermal expansion temperature coefficient of reactivity. This result fits well with integral EBR-II reactivity coefficient measurements.

With the Godiva-IV comparison results and the EBR-II results, it can be concluded that coupling finite element analysis with Monte Carlo neutron transport analysis provides a powerful technique for determining important reactor safety parameters. The technique can be applied to existing reactors and reactors proceeding through the design process which gives reactor operators and designers greater confidence in reactor operating characteristics and safety margins.