*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 highend 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 assemblyby-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

#### *Nuclear Reactor Thermal Expansion Reactivity Effect Determination Using Finite Element… DOI: http://dx.doi.org/10.5772/intechopen.93762*

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 nonanalog 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 1012 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.

$$
\Delta \rho = \frac{\left(k\_2 - k\_1\right)}{k\_2 k\_1} \tag{3}
$$

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).

$$\mathbf{y} = m\boldsymbol{\omega} + \mathbf{b} \tag{4}$$

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.
