Task 1 - enter the names of materials

This activity allows you to put you neutronics knowledge to the test

We are going to simulate a tokamak where you get to decide the materials.

The goal is to:

  • maximize the Tritium Breeding Ratio

  • minimize the heating on the center colum

  • maximize the heat in the blanket

I recommend plotting the cross sections before you make your choices.

The example for plotting cross sections might be useful.

First import OpenMC and configure the nuclear data path

import openmc
from pathlib import Path
# Setting the cross section path to the correct location in the docker image.
# If you are running this outside the docker image you will have to change this path to your local cross section path.
openmc.config['cross_sections'] = Path.home() / 'nuclear_data' / 'cross_sections.xml'

Task 1 - enter the names of materials#

# options for coolant

water = openmc.Material(name='water')
water.add_element('H', 2)
water.add_element('O', 1)
water.set_density('g/cm3', 1)

helium = openmc.Material(name='helium')
helium.add_element('He', 1)
helium.set_density('g/cm3', 0.0014)

super_critical_co2 = openmc.Material(name='super critical co2')
super_critical_co2.add_element('C', 1)
super_critical_co2.add_element('O', 2)
super_critical_co2.set_density('g/cm3', 0.001)


# options for first wall

tungsten = openmc.Material(name="your code here")
tungsten.add_element('W', 1, percent_type='wo')
tungsten.set_density('g/cm3', 19.2)

steel = openmc.Material(name="your code here")
steel.add_element('Fe', 0.95, percent_type='wo')
steel.add_element('C', 0.95, percent_type='wo')
steel.set_density('g/cm3', 1)

silicon_carbide = openmc.Material(name="your code here")
silicon_carbide.add_element('Si', 1)
silicon_carbide.add_element('C', 1)
silicon_carbide.set_density('g/cm3', 3.21)

zirconium = openmc.Material(name="your code here")
zirconium.add_element('Zr', 1, percent_type='wo')
zirconium.set_density('g/cm3', 6.49)

# options for reflector

beryllium = openmc.Material(name='mat_reflector')
beryllium.add_element('Be', 1, percent_type='wo')
beryllium.set_density('g/cm3', 1.85)

tungsten = openmc.Material(name='mat_reflector')
tungsten.add_element('W', 1, percent_type='wo')
tungsten.set_density('g/cm3', 19.2)

graphite = openmc.Material(name='mat_reflector')
graphite.add_element('C', 1, percent_type='wo')
graphite.set_density('g/cm3', 2.2)

# options for breeder

# note you can change the enrichment
lithium_lead = openmc.Material(name="your code here")
lithium_lead.add_element('Pb', 84.2, percent_type='ao')
lithium_lead.add_element('Li', 15.8, percent_type='ao', enrichment=90.0, enrichment_target='Li6', enrichment_type='ao')
lithium_lead.set_density('g/cm3', 9.5)

# note you can change the enrichment
lithium_orthosilicate = openmc.Material(name="your code here")
lithium_orthosilicate.add_element('Si', 4, percent_type='ao')
lithium_orthosilicate.add_element('O', 1, percent_type='ao')
lithium_orthosilicate.add_element('Li', 4, percent_type='ao', enrichment=50.0, enrichment_target='Li6', enrichment_type='ao')
lithium_orthosilicate.set_density('g/cm3', 2.2)

# note you can change the enrichment
lithium_titanate = openmc.Material(name="your code here")
lithium_titanate.add_element('O', 12, percent_type='ao')
lithium_titanate.add_element('Ti', 5, percent_type='ao')
lithium_titanate.add_element('Li', 4, percent_type='ao', enrichment=40.0, enrichment_target='Li6', enrichment_type='ao')
lithium_titanate.set_density('g/cm3', 2.4)

# note you can change the enrichment
liquid_lithium = openmc.Material(name="your code here")
liquid_lithium.add_element('Li', 1, percent_type='ao', enrichment=7.0, enrichment_target='Li6', enrichment_type='ao')
liquid_lithium.set_density('g/cm3', 0.5)

# conductor material

mat_conductor = openmc.Material(name='mat_conductor')
mat_conductor.add_element('Nb', 1, percent_type='ao')
mat_conductor.add_element('Sn', 3, percent_type='ao')
mat_conductor.set_density('g/cm3', 8.96)

Task 2 - plot the tritium production cross sections of all the breeder options#

# your code goes here

Task 3 - plot the absorption of all the coolants and first wall options#

# your code goes here

Task 4 - select your materials#

# put in your chosen material option
my_materials = openmc.Materials([mat_conductor]) # extend this list with your code

Task 3 - Fill the cells with the materials#

# surfaces
central_column_surface_outer = openmc.ZCylinder(r=100)
central_column_surface_mid = openmc.ZCylinder(r=95)
central_column_surface_inner = openmc.ZCylinder(r=90)

inner_sphere_surface = openmc.Sphere(r=495)
middle_sphere_surface = openmc.Sphere(r=500) 
outer_sphere_surface = openmc.Sphere(r=505)
outer_outer_sphere_surface = openmc.Sphere(r=600)
edge_of_simulation_surface = openmc.Sphere(r=700, boundary_type='vacuum')

# regions
central_column_region = -central_column_surface_inner & -edge_of_simulation_surface
central_column_coolant_region = +central_column_surface_inner & -central_column_surface_mid & -edge_of_simulation_surface
central_column_fw_region = +central_column_surface_mid & -central_column_surface_outer & -edge_of_simulation_surface

inner_vessel_region = +central_column_surface_outer & -inner_sphere_surface

blanket_fw_region = -middle_sphere_surface & +inner_sphere_surface & +central_column_surface_outer
blanket_coolant_region = +middle_sphere_surface & -outer_sphere_surface & +central_column_surface_outer
blanket_breeder_region = +outer_sphere_surface & -outer_outer_sphere_surface & +central_column_surface_outer
blanket_reflector_region = +outer_outer_sphere_surface & -edge_of_simulation_surface & +central_column_surface_outer

# cells
central_column_cell = openmc.Cell(region=central_column_region, fill=mat_conductor)
central_column_coolant_cell = openmc.Cell(region=central_column_coolant_region) # add the fill argument with your selected material)
central_column_fw_cell = openmc.Cell(region=central_column_fw_region) # add the fill argument with your selected material)
inner_vessel_cell = openmc.Cell(region=inner_vessel_region)
blanket_fw_cell = openmc.Cell(region=blanket_fw_region) # add the fill argument with your selected material)
blanket_coolant_cell = openmc.Cell(region=blanket_coolant_region) # add the fill argument with your selected material)
blanket_breeder_cell  = openmc.Cell(region=blanket_breeder_region) # add the fill argument with your selected material)
blanket_reflector_cell = openmc.Cell(region=blanket_reflector_region) # add the fill argument with your selected material)

my_geometry = openmc.Geometry([
    central_column_cell,
    central_column_coolant_cell,
    central_column_fw_cell,
    inner_vessel_cell,
    blanket_fw_cell,
    blanket_coolant_cell,
    blanket_breeder_cell,
    blanket_reflector_cell,
])

task 4 - plot the geometry#

# your code here

Task 5 - plot the source location and energy#

# initialises a new source object
my_source = openmc.IndependentSource()
# the distribution of radius is just a single value
radius = openmc.stats.Discrete([300], [1])
# the distribution of source z values is just a single value
z_values = openmc.stats.Discrete([0], [1])
# the distribution of source azimuthal angles values is a uniform distribution between 0 and 2 Pi
angle = openmc.stats.Uniform(a=0., b=2* 3.14159265359)
# this makes the ring source using the three distributions and a radius
my_source.space = openmc.stats.CylindricalIndependent(r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0))
# sets the direction to isotropic
my_source.angle = openmc.stats.Isotropic()
# sets the energy distribution to a Muir distribution neutrons
my_source.energy = openmc.stats.muir(e0=14080000.0, m_rat=5.0, kt=20000.0)

my_settings = openmc.Settings()
my_settings.batches = 10
my_settings.inactive = 0
my_settings.particles = 500
my_settings.run_mode = 'fixed source'
my_settings.source = my_source

# plot the source term
# your code goes here

Task 7 - complete the tally scores and filters#

cell_filter_breeder = openmc.CellFilter(blanket_breeder_cell)
tbr_tally = openmc.Tally(name='TBR')
tbr_tally.filters = [cell_filter_breeder]
# tbr_tally.scores = [# your code here]

cell_filter_blanket_fw_coolant = openmc.CellFilter([
    blanket_breeder_cell,
    blanket_coolant_cell,
    blanket_fw_cell,
    central_column_coolant_cell,
    central_column_fw_cell,
])
blanket_heating_tally = openmc.Tally(name='heating')
blanket_heating_tally.filters = [cell_filter_blanket_fw_coolant]
blanket_heating_tally.scores = ['heating']

# material_filter_conductor = openmc.MaterialFilter() # add the material to the material filter
conductor_damage_tally = openmc.Tally(name='conductor_damage')
# conductor_damage_tally.filters = [material_filter_conductor]
# conductor_damage_tally.scores = [] # set the score of the tally

my_tallies = openmc.Tallies([tbr_tally, blanket_heating_tally, conductor_damage_tally])

model = openmc.model.Model(my_geometry, my_materials, my_settings) # add the tallies in here

# removes the old output files
!rm *.h5 *.xml

sp_filename = model.run()

                                %%%%%%%%%%%%%%%
                           %%%%%%%%%%%%%%%%%%%%%%%%
                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    %%%%%%%%%%%%%%%%%%%%%%%%
                                     %%%%%%%%%%%%%%%%%%%%%%%%
                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%
                ##################     %%%%%%%%%%%%%%%%%%%%%%%
                ###################     %%%%%%%%%%%%%%%%%%%%%%%
                ####################     %%%%%%%%%%%%%%%%%%%%%%
                #####################     %%%%%%%%%%%%%%%%%%%%%
                ######################     %%%%%%%%%%%%%%%%%%%%
                #######################     %%%%%%%%%%%%%%%%%%
                 #######################     %%%%%%%%%%%%%%%%%
                 ######################     %%%%%%%%%%%%%%%%%
                  ####################     %%%%%%%%%%%%%%%%%
                    #################     %%%%%%%%%%%%%%%%%
                     ###############     %%%%%%%%%%%%%%%%
                       ############     %%%%%%%%%%%%%%%
                          ########     %%%%%%%%%%%%%%
                                      %%%%%%%%%%%

                 | The OpenMC Monte Carlo Code
       Copyright | 2011-2025 MIT, UChicago Argonne LLC, and contributors
         License | https://docs.openmc.org/en/latest/license.html
         Version | 0.15.3-dev144
     Commit Hash | 14ab4d3b6f3d77b73e8bb480782bbe1bdd56483d
       Date/Time | 2025-05-14 14:07:12
  OpenMP Threads | 16

 Reading model XML file 'model.xml' ...
 Reading cross sections XML file...
 Reading Nb93 from /home/jon/nuclear_data/neutron/Nb93.h5
 Reading Sn112 from /home/jon/nuclear_data/neutron/Sn112.h5
 Reading Sn114 from /home/jon/nuclear_data/neutron/Sn114.h5

 Reading Sn115 from /home/jon/nuclear_data/neutron/Sn115.h5
 Reading Sn116 from /home/jon/nuclear_data/neutron/Sn116.h5
 Reading Sn117 from /home/jon/nuclear_data/neutron/Sn117.h5
 Reading Sn118 from /home/jon/nuclear_data/neutron/Sn118.h5
 Reading Sn119 from /home/jon/nuclear_data/neutron/Sn119.h5
 Reading Sn120 from /home/jon/nuclear_data/neutron/Sn120.h5

 Reading Sn122 from /home/jon/nuclear_data/neutron/Sn122.h5
 Reading Sn124 from /home/jon/nuclear_data/neutron/Sn124.h5
 Minimum neutron data temperature: 294.0 K
 Maximum neutron data temperature: 294.0 K
 Preparing distributed cell instances...
 Writing summary.h5 file...
 Maximum neutron transport energy: 20000000.0 eV for Sn112

 ===============>     FIXED SOURCE TRANSPORT SIMULATION     <===============

 Simulating batch 1
 Simulating batch 2
 Simulating batch 3
 Simulating batch 4
 Simulating batch 5
 Simulating batch 6
 Simulating batch 7

 Simulating batch 8
 Simulating batch 9
 Simulating batch 10
 Creating state point statepoint.10.h5...

 =======================>     TIMING STATISTICS     <=======================

 Total time for initialization     = 2.5735e-01 seconds
   Reading cross sections          = 2.4667e-01 seconds
 Total time in simulation          = 1.5212e-02 seconds
   Time in transport only          = 1.2982e-02 seconds
   Time in active batches          = 1.5212e-02 seconds
   Time accumulating tallies       = 2.4440e-06 seconds
   Time writing statepoints        = 2.0069e-03 seconds
 Total time for finalization       = 1.2130e-06 seconds
 Total time elapsed                = 2.7629e-01 seconds
 Calculation Rate (active)         = 328681.0 particles/second

 ============================>     RESULTS     <============================

 Leakage Fraction            = 1.02800 +/- 0.00251

Task 8 - complete the code to get the correct tally names#

sp = openmc.StatePoint(sp_filename)

# tbr_tally = sp.get_tally() # add the name of the tally to help find the tally

# print(f'TBR={tbr_tally.mean.flatten()[0]} with standard deviation of {tbr_tally.std_dev.flatten()[0]}')

# heating_tally = sp.get_tally() # add the name of the tally to help find the tally
# print(f'Heating={heating_tally.mean.flatten()[0]/1e6}MeV per source neutron with standard deviation of {heating_tally.std_dev.flatten()[0]}')

# damage_tally = sp.get_tally() # add the name of the tally to help find the tally
# print(f'damage={damage_tally.mean.flatten()[0]} with standard deviation of {damage_tally.std_dev.flatten()[0]}')
Contents