Pin-wise fuel temperature distributions are not reflected

[fshriver]

Hi,

I've created a standard 17x17 PWR assembly with all steady-state conditions EXCEPT for a changing fuel temperature distribution across the assembly. The fuel temperature distribution is implemented at the material level, i.e. each pin has its own separate fuel material with its own separate temperature. The problem is the pin powers I'm getting (collected by following the notebook here) don't reflect the temperature distribution at all, and stay relatively flat. As a sanity check I'm passing in a 1,000 Kelvin difference for the fuel temperature distribution (500 to 1500) but the pin powers still don't reflect any change.

I can see the material temperature values given in the materials.xml file, and I can also confirm from the xml files and the code output that I'm using the WMP option. However, nothing I've been able to do so far seems to cause any changes to the pin power distribution due to fuel temperature. Any help on this issue would be appreciated.

[gridley]

Hey Forrest, hope you're doing well man.

So.. here's my thought. You linked to the multigroup page. Are you running in multigroup mode? The C5G7 XS used there have no temperature dependence, so this might make sense if that's the case.

[fshriver]

Hey Gavin! I'm doing well and so is Fan, hope your PhD is going well (mine's almost done, light at the end of the tunnel).

Unfortunately I don't think that's it, I'm running OpenMC in continuous energy mode. I just used the multigroup notebook to figure out what I needed to do to set up the mesh filters and tally scores. Otherwise, all cross sections should be pulling from the default OpenMC ENDF/B-7.1, and I'm creating all materials loosely following the VERA benchmark description for their assembly problem (here).

This is all being done using OpenMC 0.12.0, by the way. I'm putitng up a small example script demonstrating the issue below, if you want to take a look at it. From plotting the results I maybe see some SLIGHT differences from one quarter of the assembly to the other, however there doesn't appear to be a strong trend matching the significantly changing temperature distribution I'm inputting.

import openmc
import numpy as np
import h5py

#####
# Pre-set constants (including geometry)
#####
active_batches = 1000
inactive_batches = 200
nparticles = 10000

cr_indices = [
    (5,2), (8,2), (11,2),
    (3,3), (13,3),
    (2,5), (5,5), (8,5), (11,5), (14,5),
    (2,8), (5,8), (8,8), (11,8), (14,8),
    (2,11), (5,11), (8,11), (11,11), (14,11),
    (3,13), (13,13),
    (5,14), (8,14), (11,14),
]

GT_IR = openmc.ZCylinder(r = 0.561)
GT_OR = openmc.ZCylinder(r = 0.602)

FP_OR = openmc.ZCylinder(r = 0.4096)
CD_IR = openmc.ZCylinder(r = 0.418)
CD_OR = openmc.ZCylinder(r = 0.475)

FT_MIN = 500
FT_MAX = 1500

#####
# Create a fuel temperature distribution
#####
ft_distribution = np.zeros(shape = (17,17))
for i in range(17):
    for j in range(17):
        scale = (i*j)/289 # Entirely artificial, just done to create a smooth distribution
        ft_distribution[i,j] = FT_MIN * (1 - scale) + FT_MAX * scale

#####
# Material definitions - mostly based off of VERA Core Physics Benchmark 
# Progression Problems
#####
fuel_proto = openmc.Material()
fuel_proto.set_density('sum')
fuel_proto.add_nuclide('U234', 6.11864E-06)
fuel_proto.add_nuclide('U235', 7.18132E-04)
fuel_proto.add_nuclide('U236', 3.29861E-06)
fuel_proto.add_nuclide('U238', 2.21546E-02)
fuel_proto.add_nuclide('O16',  4.57642E-02)

# Using shortcut functions to save screen space
cladding = openmc.Material()
cladding.set_density('g/cm3', 6.56)
cladding.add_element('Zr', 100.0)
cladding.temperature = 650

moderator = openmc.model.borated_water(700, temperature = 550, pressure = 15.16847)

#####
# Create the assembly geometry
#####
universes = np.ndarray(shape = (17,17), dtype = object)
for i in range(17):
    for j in range(17):
        if (i,j) in cr_indices:
            inner_mod_cell = openmc.Cell(fill = moderator, region = -GT_IR)
            clad_cell = openmc.Cell(fill = cladding, region = +GT_IR & -GT_OR)
            outer_mod_cell = openmc.Cell(fill = moderator, region = +GT_OR)

            pin = openmc.Universe()
            pin.add_cells([inner_mod_cell, clad_cell, outer_mod_cell])
            universes[i,j] = pin
        else:
            # Cloning the material so we can set the temperature separately
            fuel = fuel_proto.clone()
            fuel.temperature = ft_distribution[i,j]
            fuel_cell = openmc.Cell(fill = fuel, region = -FP_OR)
            gap_cell = openmc.Cell(fill = None, region = +FP_OR & -CD_IR)
            clad_cell = openmc.Cell(fill = cladding, region = +CD_IR & -CD_OR)
            moderator_cell = openmc.Cell(fill = moderator, region = +CD_OR)

            pin = openmc.Universe()
            pin.add_cells([fuel_cell, gap_cell, clad_cell, moderator_cell])
            universes[i,j] = pin

lattice = openmc.RectLattice()
lattice.pitch = (1.26, 1.26)
lattice.lower_left = (-21.42/2, -21.42/2)
lattice.universes = universes

#####
# Set up boundary conditions and the root cell/universe
#####
min_x = openmc.XPlane(x0 = -(21.42)/2, boundary_type = 'reflective')
max_x = openmc.XPlane(x0 = +(21.42)/2, boundary_type = 'reflective')
min_y = openmc.YPlane(y0 = -(21.42)/2, boundary_type = 'reflective')
max_y = openmc.YPlane(y0 = +(21.42)/2, boundary_type = 'reflective')
root_cell = openmc.Cell(fill = lattice)
root_cell.region = +min_x & -max_x & +min_y & -max_y

root_universe = openmc.Universe()
root_universe.add_cell(root_cell)

model = openmc.model.Model()
model.geometry.root_universe = root_universe

#####
# Model settings - use at your own risk, this is meant to be run on a compute
# node.
#####
model.settings.batches = active_batches
model.settings.inactive = inactive_batches
model.settings.particles = nparticles
model.settings.source = openmc.Source(space = openmc.stats.Box(
    [-21.42/2, -21.42/2, -1],
    [21.42/2, 21.42/2, 1],
    only_fissionable = True
))
model.settings.run_mode = 'eigenvalue'
model.settings.temperature = {
    'method': 'interpolation',
    'multipole': True
}

#####
# Set up tallies for pin powers
#####
tallies = openmc.Tallies()

mesh = openmc.RegularMesh()
mesh.dimension = [17, 17]
mesh.lower_left = [-21.42/2, -21.42/2]
mesh.upper_right = [+21.42/2, +21.42/2]
mesh_filter = openmc.MeshFilter(mesh)

tally = openmc.Tally(name = 'pin_powers')
tally.filters = [mesh_filter]
tally.scores = ['fission']

tallies.append(tally)
model.tallies = tallies

#####
# Create a plot for convenience - used to sanity check geometry/materials
#####
plot = openmc.Plot()
plot.origin = (0.0, 0.0, 0.0)
plot.width = (21.42, 21.42)
plot.pixels = (1000,1000)
plot.color_by = 'material'
model.plots.append(plot)

#####
# Export the model and run inside of the current Python session
#####
model.export_to_xml()
model.run(threads = 32)

#####
# Get the fission rates/pin powers from the statepoint file and export them as 
# well as the pinwise fuel temperature distribution to an HDF5 file, the most 
# useful of all file types. Use whatever plotting utility you like to examine 
# the results; I use seaborn with the viridis map, personally.
#####
statepoint = openmc.StatePoint('statepoint.{}.h5'.format(active_batches))
mesh_tally = statepoint.get_tally(name='pin_powers')
fission_rates = mesh_tally.get_values(scores=['fission'])
fission_rates.shape = [17, 17]
fission_rates /= np.mean(fission_rates[fission_rates > 0.])

outfile = h5py.File('issue1663.h5', 'w')
outfile.create_dataset('pin_powers', data = fission_rates)
outfile.create_dataset('ft_distribution', data = ft_distribution)

[gridley]

Thanks man, glad to hear it. I'm still working on a master's degree though. xD

Hm, can you paste your material.xml in too? I think that might also reveal some hints. Also, have you tried not using multipole? Personally I've not tried running with multipole before and have only used stochastic interpolation when looking at temperature-dependent problems.

[fshriver]

Sure, I've attached the xml file here. materials.txt

From examining it apparently I screwed up slightly on the fuel temperature generation indices (0-16, not 1-17, doh!) but overall the distribution is almost a delta of 1,000 Kelvin, and this is reflected in the xml file (scroll down a bit). I tried turning off multipole, however no dice; there still doesn't appear to be a significant change in the pin powers.

[gridley]

Yeah, either way, that's unexpected. To be honest, I don't have a feel for exactly how much the power distribution should change with this interesting temperature profile. Have you compared k_eff between the uniform temperature and linear profile temperature case? If you see no feedback there, that's a dead giveaway that something is fishy.

[gridley]

Hey @fshriver, what was the verdict on this?

[gridley]

closing due to inactivity