Tutorial: Pin Cell Model

Overview

The fuel pin cell is the most fundamental unit of a pressurized water reactor core. This tutorial walks through constructing a PWR pin cell model in Serpent, from material definitions through geometry specification to criticality calculation and post-processing.

The model represents a single fuel pin surrounded by moderator, with reflective boundary conditions imposed on all faces to simulate an infinite square lattice of identical pins. This idealization eliminates neutron leakage and yields the infinite multiplication factor k-infinity, a quantity that characterizes the intrinsic reactivity of the fuel-moderator system independent of core size and shape.

Model Specifications

The pin cell geometry follows a standard Westinghouse 17x17 fuel assembly design. The fuel pellet is a cylinder of uranium dioxide with a radius of 0.4096 cm, surrounded by a thin helium-filled gap extending to an inner cladding radius of 0.4178 cm. The Zircaloy-4 cladding has an outer radius of 0.4750 cm. The square unit cell has a pitch of 1.26 cm, corresponding to the rod-to-rod spacing in the assembly lattice. The active fuel height is 365.76 cm, and the fuel enrichment is set to 4.5 atom percent U-235. Operating conditions are representative of hot full power, with a coolant temperature of approximately 310 degrees Celsius at a system pressure of 15.5 MPa. These parameters collectively define a fuel-to-moderator ratio that is slightly under-moderated, as is characteristic of commercial PWR designs optimized for negative moderator temperature coefficients.

Step 1: Material Definitions

The model requires three principal materials: uranium dioxide fuel, Zircaloy-4 cladding, and light water moderator, along with a helium fill gas in the pellet-cladding gap. Each material is specified with a physical density (negative sign convention in Serpent) and an operating temperature for Doppler broadening. The fuel operates near 900 K due to the large temperature gradient across the pellet, the cladding and gap sit near 600 K, and the moderator is at 574 K corresponding to the bulk coolant temperature.

The thermal scattering treatment for hydrogen bound in water is essential for correct neutron thermalization. Without the S(alpha, beta) library, the free-gas scattering kernel would overestimate the energy loss per collision at thermal energies, leading to a softer spectrum and a systematic bias in the multiplication factor on the order of several hundred pcm.

% ========================================
% PWR Pin Cell Model - Materials
% ========================================

% UO2 fuel pellet (4.5% enriched at 900K)
mat fuel -10.4 tmp 900
92235.09c  0.045        % 4.5 at% U-235
92238.09c  0.955        % 95.5 at% U-238  
8016.09c   2.0          % Stoichiometric oxygen (2 atoms per UO2)

% Helium gap (low pressure)
mat gap -0.0018 tmp 600
2004.09c -1.0           % Helium-4

% Zircaloy-4 cladding at 600K
% Simplified: using single isotopes instead of natural element compositions
mat clad -6.56 tmp 600
40090.09c -0.9845       % Zirconium-90 (98.45%)
50120.09c -0.0155       % Tin-120 (1.55%)

% Light water moderator at 574K (301°C)
mat water -0.714 tmp 574 moder lwtr 1001
1001.09c  2.0           % Hydrogen-1
8016.09c  1.0           % Oxygen-16

% Thermal scattering library for H in light water
therm lwtr lwj3.22t

Step 2: Geometry Definition

The geometry is constructed using Serpent's combinatorial solid geometry (CSG) approach, in which regions of space are defined as intersections and unions of half-spaces bounded by elementary surfaces. For the pin cell, three concentric cylinders define the fuel pellet outer surface, the gap outer surface (equivalently the cladding inner surface), and the cladding outer surface. Six planar surfaces define the faces of the rectangular unit cell. Each cell card assigns a material to a region described by a Boolean combination of surface senses, where a negative sense indicates the interior of a surface and a positive sense indicates the exterior.

The moderator cell is defined as the region outside the cladding but inside all six bounding planes. The outside cell captures all space beyond the bounding surfaces and is required by Serpent to close the geometry. Reflective boundary conditions applied later ensure that neutrons reaching the outer planes are specularly reflected back into the model, simulating the infinite lattice environment.

% ========================================
% PWR Pin Cell Model - Geometry
% ========================================

% Surface definitions (all dimensions in cm)
surf fuel_radius cylz 0.0 0.0 0.4096    % Fuel pellet outer radius
surf gap_radius  cylz 0.0 0.0 0.4178    % Gap outer radius (clad inner)  
surf clad_radius cylz 0.0 0.0 0.4750    % Clad outer radius

% Square lattice boundaries
surf left   px -0.63    % Left boundary (half pitch)
surf right  px  0.63    % Right boundary  
surf front  py -0.63    % Front boundary
surf back   py  0.63    % Back boundary
surf bottom pz -182.88  % Bottom boundary (half height)
surf top    pz  182.88  % Top boundary

% Cell definitions (cell NAME UNIVERSE MATERIAL SURFACES)
cell 1 0 fuel  -fuel_radius
cell 2 0 gap    fuel_radius -gap_radius  
cell 3 0 clad   gap_radius -clad_radius
cell 4 0 water  clad_radius left -right front -back bottom -top

% Outside boundary (implicit complement)
cell 5 0 outside

Step 3: Boundary Conditions and Settings

Setting bc 2 imposes reflective (mirror) boundary conditions on all outer surfaces of the root universe, simulating an infinite lattice from a single unit cell. The criticality source is configured with set pop, which takes three arguments: the number of neutrons per cycle, the number of active cycles over which tallies are accumulated, and the number of inactive cycles used to converge the fission source distribution before tallying begins.

A population of 20,000 neutrons per cycle with 200 active and 50 inactive cycles is sufficient for a simple pin cell. The Shannon entropy diagnostic monitors the spatial convergence of the fission source and should stabilize to a constant value well before the end of the inactive cycles. If the entropy is still drifting when active cycles begin, the number of inactive cycles should be increased.

% ========================================
% PWR Pin Cell Model - Boundary Conditions
% ========================================

% Reflective boundary conditions (infinite lattice)
set bc 2    % 2 = reflective boundaries on all universe boundaries

% Criticality calculation parameters
set pop 20000 200 50   % 20k neutrons/cycle, 200 active, 50 inactive

% Enable source convergence monitoring
set entr 10 10 1       % Shannon entropy mesh (NX NY NZ)

TMS (target motion sampling) Doppler broadening ensures that cross sections reflect the thermal motion of target nuclei at each material's assigned temperature. The unresolved resonance probability table treatment improves self-shielding accuracy in the intermediate energy range where individual resonances are not resolved in the evaluated data.

% Physics options for accurate PWR calculations
set title "PWR Pin Cell k-infinity Calculation"

% Neutron physics
% TMS is configured per-material via the tmp keyword on the mat card
set ures 1               % Unresolved resonance treatment

% Memory and performance
set memfrac 0.8          % Use 80% of available memory
% Threads set via command line: sss2 -omp 4 input

Step 4: Complete Input File

The following listing consolidates all preceding cards into a single input file. Save it as pin_cell.inp and execute with sss2 pin_cell.inp. On a modern workstation with four threads, the calculation should complete in a few minutes.

% ========================================
% PWR Pin Cell k-infinity Calculation
% Created for Nuclear Monte Carlo Guide
% ========================================

set title "PWR Pin Cell k-infinity Calculation"

% ========================================
% MATERIALS
% ========================================

% UO2 fuel pellet (4.5% enriched at 900K)
mat fuel -10.4 tmp 900
92235.09c  0.045        % 4.5 at% U-235
92238.09c  0.955        % 95.5 at% U-238  
8016.09c   2.0          % Stoichiometric oxygen (2 atoms per UO2)

% Helium gap 
mat gap -0.0018 tmp 600
2004.09c -1.0           % Helium-4

% Zircaloy-4 cladding at 600K
% Simplified: using single isotopes instead of natural element compositions
mat clad -6.56 tmp 600
40090.09c -0.9845       % Zirconium-90
50120.09c -0.0155       % Tin-120

% Light water moderator at 574K
mat water -0.714 tmp 574 moder lwtr 1001
1001.09c  2.0           % Hydrogen-1
8016.09c  1.0           % Oxygen-16

therm lwtr lwj3.22t     % Thermal scattering for H in water

% ========================================
% GEOMETRY
% ========================================

% Cylindrical surfaces (cm)
surf fuel_radius cylz 0.0 0.0 0.4096
surf gap_radius  cylz 0.0 0.0 0.4178  
surf clad_radius cylz 0.0 0.0 0.4750

% Square lattice boundaries
surf left   px -0.63
surf right  px  0.63    
surf front  py -0.63
surf back   py  0.63    
surf bottom pz -182.88  
surf top    pz  182.88  

% Cell definitions
cell 1 0 fuel  -fuel_radius
cell 2 0 gap   +fuel_radius -gap_radius  
cell 3 0 clad  +gap_radius -clad_radius
cell 4 0 water +clad_radius +left -right +front -back +bottom -top
cell 5 0 outside

% ========================================
% BOUNDARY CONDITIONS
% ========================================

set bc 2                 % Reflective boundaries

% ========================================  
% PHYSICS SETTINGS
% ========================================

% TMS is configured per-material via the tmp keyword on the mat card
set ures 1               % Unresolved resonance treatment

% ========================================
% CRITICALITY CALCULATION
% ========================================

set pop 20000 200 50     % Population, active cycles, inactive cycles
set entr 10 10 1         % Shannon entropy mesh (NX NY NZ)

% ========================================
% OUTPUT OPTIONS  
% ========================================

set his 1                % History output
set pcc 0                % Predictor-corrector calculation off

Step 5: Running the Calculation

Serpent prints cycle-by-cycle estimates of k-effective during the run, allowing you to watch source convergence during inactive cycles. Upon completion, the code generates several output files: the main output file (pin_cell.out) containing the final k-effective estimate and summary statistics, the results file (pin_cell_res.m) in MATLAB-readable format for post-processing, the random number seed file, and the cycle history file if enabled.

# Run Serpent calculation
sss2 pin_cell.inp

# Monitor calculation progress
tail -f pin_cell.out

# Key output files generated:
# pin_cell.out     - Main output with k-eff results
# pin_cell_res.m   - Results in MATLAB format  
# pin_cell.seed    - Random number seed
# pin_cell.his     - History file

For a fresh 4.5%-enriched pin cell with reflective boundaries, the expected k-infinity falls in the range of approximately 1.31 to 1.35, well above unity because the infinite lattice assumption eliminates all neutron leakage. The statistical uncertainty with 20,000 neutrons per cycle and 200 active cycles should be below 10 pcm (0.00010 in delta-k). The four-factor formula components provide physical insight into the result: the thermal utilization factor is typically 0.92 to 0.95, reflecting the competition between fuel absorption and parasitic capture in the cladding and moderator; the resonance escape probability is 0.86 to 0.89, governed by U-238 resonance absorption during neutron slowing down; and the fast fission factor is 1.02 to 1.04, a modest contribution from fissions induced by fast neutrons in U-238. In a real reactor, criticality is achieved at k-effective equal to unity through the combined effects of neutron leakage from the finite core, control rod insertion, and soluble boron in the coolant.

Step 6: Analysis and Visualization

The results file uses MATLAB syntax and loads directly into MATLAB or Octave. A plot of k-effective versus cycle number should show random fluctuations about a stable mean after the inactive cycles, with no visible drift or autocorrelation.

% Load results in MATLAB
run('pin_cell_res.m');

% Extract key parameters (workspace variables from Serpent output)
keff = IMP_KEFF(1);
keff_err = IMP_KEFF(2);

fprintf('k-infinity = %.5f ± %.5f\n', keff, keff_err);

% Load history data (Serpent writes .m scripts, not flat files)
run('pin_cell_his0.m');

% Plot k-eff convergence
figure;
cycles = 1:size(HIS_IMP_KEFF, 1);
plot(cycles, HIS_IMP_KEFF(:,1));
xlabel('Cycle Number');
ylabel('k-effective');
title('Source Convergence');
grid on;

% Check for convergence
final_cycles = HIS_IMP_KEFF(end-100:end, 1);  % Last 100 cycles
std_dev = std(final_cycles);
fprintf('Standard deviation (last 100 cycles): %.1e\n', std_dev);

The plot card produces a color-coded cross-section image of the geometry at a specified plane for verifying surfaces and cells before running transport. Detector cards define spatial and energy-resolved tallies for flux, reaction rates, and surface currents.

% Add to your input file for visualization
plot 3 500 500          % 500x500 pixel plot

% Energy grid for flux spectrum (energies in MeV)
ene flux_grid 1 1E-11 0.625E-6 20.0

% Detector for flux mapping
det flux_map dm fuel de flux_grid dx -0.63 0.63 50 dy -0.63 0.63 50

% Surface current detectors  
det fuel_current ds fuel_radius du 0
det clad_current ds clad_radius du 0

Extensions and Parameter Studies

Varying the U-235 enrichment from 3 to 5 atom percent reveals the near-linear relationship between enrichment and excess reactivity in the fresh fuel limit. Sweeping the fuel temperature from 600 K to 1200 K while holding all other parameters fixed isolates the Doppler reactivity coefficient, which should be negative due to increased resonance absorption in U-238 at elevated temperatures. Changing the pin pitch or pellet radius alters the hydrogen-to-heavy-metal ratio and demonstrates the concept of optimum moderation, where the multiplication factor passes through a maximum as the moderator fraction increases. Adding soluble boron to the moderator material at concentrations of 500 to 2000 ppm illustrates the linear worth of chemical shim and its effect on the moderator temperature coefficient.

Extending the pin cell to depletion requires a power normalization card, a depletion schedule, and optionally branch calculations for generating parameterized few-group cross sections.

% Burnup calculation extension
set power 15.0E3         % Pin power (15 kW)

dep daystep              % Burnup time steps (days)
1 5 10 30 100 300 500 1000

% For parametric studies, define branches using the branch card
% to systematically vary fuel temperature, boron concentration,
% and other operating conditions at each burnup step