Skip to content

Commit

Permalink
Merge pull request #18 from paulromano/tally-power-normalization
Browse files Browse the repository at this point in the history
Add tally power normalization notebook
  • Loading branch information
pshriwise authored May 23, 2024
2 parents fdc8386 + 472aaeb commit e44fdce
Showing 1 changed file with 366 additions and 0 deletions.
366 changes: 366 additions & 0 deletions tally-power-normalization.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,366 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "1c9efafb",
"metadata": {},
"source": [
"# Tally Power Normalization\n",
"\n",
"In this notebook, we demonstrate how to normalize tally values in OpenMC assuming you have a reactor operating at a known power. We'll begin by creating a very simple \"reactor\" model consisting of a sphere of U235 surrounded by water."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "aa112441",
"metadata": {},
"outputs": [],
"source": [
"import openmc"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "d3e05e9d",
"metadata": {},
"outputs": [],
"source": [
"u235 = openmc.Material()\n",
"u235.add_nuclide('U235', 1.0)\n",
"u235.set_density('g/cm3', 10.0)\n",
"\n",
"water = openmc.Material()\n",
"water.add_components({'H': 2.0, 'O': 1.0})\n",
"water.add_s_alpha_beta('c_H_in_H2O')\n",
"\n",
"fuel_sph = openmc.Sphere(r=6.5)\n",
"outer_sph = openmc.Sphere(r=15.0, boundary_type='vacuum')\n",
"\n",
"fuel = openmc.Cell(fill=u235, region=-fuel_sph)\n",
"moderator = openmc.Cell(fill=water, region=+fuel_sph & -outer_sph)\n",
"geometry = openmc.Geometry([fuel, moderator])\n",
"model = openmc.Model(geometry=geometry)\n",
"\n",
"model.settings.batches = 100\n",
"model.settings.inactive = 10\n",
"model.settings.particles = 1000\n",
"model.settings.source = openmc.IndependentSource(space=openmc.stats.Point())"
]
},
{
"cell_type": "markdown",
"id": "4346c07d",
"metadata": {},
"source": [
"Let's run OpenMC to see what we get for $k_\\text{eff}$:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "e3669be6",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"k-effective = 0.9943+/-0.0030\n"
]
}
],
"source": [
"sp_path = model.run(output=False)\n",
"with openmc.StatePoint(sp_path) as sp:\n",
" keff = sp.keff\n",
"print(f'k-effective = {keff}')"
]
},
{
"cell_type": "markdown",
"id": "a48ca59a",
"metadata": {},
"source": [
"Now let's say we want to know the fission reaction rate in \\[fissions/sec\\] in the fuel. We'll create a tally for this using the \"fission\" score:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "72905d10",
"metadata": {},
"outputs": [],
"source": [
"fission_tally = openmc.Tally()\n",
"fission_tally.filters = [openmc.CellFilter([fuel])]\n",
"fission_tally.scores = ['fission']\n",
"model.tallies = [fission_tally]"
]
},
{
"cell_type": "markdown",
"id": "6f5817ae",
"metadata": {},
"source": [
"We then run OpenMC and get the result from the statepoint:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "149b4d65",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Fission rate = 0.3955335940378975 [reactions/source]\n"
]
}
],
"source": [
"sp_path = model.run(output=False)\n",
"with openmc.StatePoint(sp_path) as sp:\n",
" tally = sp.tallies[fission_tally.id]\n",
" fission_rate = tally.mean.ravel()[0]\n",
" \n",
"print(f'Fission rate = {fission_rate} [reactions/source]')"
]
},
{
"cell_type": "markdown",
"id": "cc37fdde",
"metadata": {},
"source": [
"So, the number we get here is not quite what we wanted—it's the fission reaction rate *per source particle*. We want the reaction rate *per second*. In order to determine that, we need to figure out the source rate, i.e., the number of source particles emitted per second. The method for doing this is described in the OpenMC [user's guide](https://docs.openmc.org/en/stable/usersguide/tallies.html#normalization-of-tally-results) but here we'll give a practical demonstration.\n",
"\n",
"We need two things in order to determine our desired normalization factor in units of \\[source/sec\\]. First, we need to know the power in \\[W\\] or \\[J/sec\\]. This is not something OpenMC or any other neutronics code will provide for you—that is something you have to provide as a user. In addition to the power, we need the observed heating rate from the simulation in \\[J/source\\]. To get that, we can add another tally for the heating rate. We'll use the `fission-q-recoverable` score but there are some pitfalls as to which score you choose (see further discussion below)."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "48bd4b14",
"metadata": {},
"outputs": [],
"source": [
"heating_tally = openmc.Tally()\n",
"heating_tally.scores = ['fission-q-recoverable']\n",
"model.tallies = [fission_tally, heating_tally]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "4c804fc7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Heating rate = 76235017.82874997 [eV/source]\n"
]
}
],
"source": [
"sp_path = model.run(output=False)\n",
"with openmc.StatePoint(sp_path) as sp:\n",
" tally = sp.tallies[heating_tally.id]\n",
" heating_rate_ev = tally.mean.ravel()[0]\n",
" \n",
"print(f'Heating rate = {heating_rate_ev} [eV/source]')"
]
},
{
"cell_type": "markdown",
"id": "310229cd",
"metadata": {},
"source": [
"The OpenMC tally gives us units of \\[eV/source\\] and we need to convert this to \\[J/source\\]:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "e684b4c9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Heating rate = 1.2214222086486662e-11 [J/source]\n"
]
}
],
"source": [
"joule_per_ev = 1.60218e-19\n",
"heating_rate = heating_rate_ev * joule_per_ev\n",
"print(f'Heating rate = {heating_rate} [J/source]')"
]
},
{
"cell_type": "markdown",
"id": "1e9c1002",
"metadata": {},
"source": [
"At this point, let's assume that our spherical reactor has a power of 30 kW. That is:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "a7f4ae45",
"metadata": {},
"outputs": [],
"source": [
"power = 30e3 # [J/s]"
]
},
{
"cell_type": "markdown",
"id": "03d1211b",
"metadata": {},
"source": [
"The tally normalization factor is simply the power in \\[J/s\\] divided by the heating rate in \\[J/source\\]:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "ef03e05b",
"metadata": {},
"outputs": [],
"source": [
"source_per_sec = power / heating_rate"
]
},
{
"cell_type": "markdown",
"id": "f6989d40",
"metadata": {},
"source": [
"Whenever we want our tallied quantities in units of \"per sec\" rather than \"per source neutron\", all we have to do is multiply by the normalization factor. For example, our fission rate would be:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "218ad621",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Fission rate = 9.715e+14 [reactions/sec]\n"
]
}
],
"source": [
"fission_rate_sec = fission_rate * source_per_sec\n",
"print(f'Fission rate = {fission_rate_sec:.3e} [reactions/sec]')"
]
},
{
"cell_type": "markdown",
"id": "19dd5d6c",
"metadata": {},
"source": [
"As a quick sanity check, we can multiply by 200 MeV/reaction to see how much heating we'd expect from this many fission reactions:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "a6f5943a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Expected power = 31130.1 [W]\n"
]
}
],
"source": [
"expected_power = fission_rate_sec * 200e6 * joule_per_ev\n",
"print(f'Expected power = {expected_power:.1f} [W]')"
]
},
{
"cell_type": "markdown",
"id": "f3364440",
"metadata": {},
"source": [
"The value is fairly close to the expected 30 kW but not quite exact since we used a crude estimate of 200 MeV/reaction. The actual amount of heat being deposited per reaction can be determined by dividing the heating rate by the fission rate."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "19fbe062",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Heat per reaction = 192.7 [MeV]\n"
]
}
],
"source": [
"mev_per_reaction = heating_rate_ev / fission_rate * 1e-6\n",
"print(f'Heat per reaction = {mev_per_reaction:.1f} [MeV]')"
]
},
{
"cell_type": "markdown",
"id": "d7b63bb5",
"metadata": {},
"source": [
"How much heat you observe will depend on which [score](https://docs.openmc.org/en/stable/usersguide/tallies.html#scores) you choose to estimate the heating. OpenMC has a variety of scores that allow you to tally heating that differ in subtle ways, so you should exercise some caution when picking a score. A few notes about the available scores:\n",
"\n",
"- `\"heating\"` — This score gives you the heating from *all* reactions, not just fission, so it gives you a way of accounting for heat from (n,$\\gamma$) reactions. It will also include photon energy deposition but only if you are running a coupled neutron–photon calculation; for a neutron-only calculation, you should use the `\"heating-local\"` score to account for photon energy deposition.\n",
"- `\"heating-local\"` — This score is the same as `\"heating\"` except that it assumes photon energy is deposited locally and thus can be used in a neutron-only calculation. This score requires special cross section data (MT=901) that is not present in all data libraries, particularly those that were converted straight from ACE files. The \"official\" data libraries from https://openmc.org do have the necessary data included.\n",
"- `\"fission-q-recoverable\"` — This score represents the recoverable energy release from fission, which is basically the total fission energy release less neutrinos. Do note that it does not account for heating from non-fission reactions.\n",
"- `\"kappa-fission\"` — This older score is basically equivalent to `\"fission-q-recoverable\"` but does not include an energy-dependence on the incident neutron energy."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "35f457bf",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}

0 comments on commit e44fdce

Please sign in to comment.