import openmc
import os
# MATERIALS
breeder_material = openmc.Material() # Pb84.2Li15.8
# breeder_material.add_element('Pb', 84.2, percent_type='ao')
breeder_material.add_element('Li', 1.) # natural enrichment = 7% Li6
breeder_material.set_density('g/cm3', 0.01) # around 11 g/cm3
copper_material = openmc.Material()
copper_material.set_density('g/cm3', 0.01)
copper_material.add_element('Li', 1.0)
eurofer_material = openmc.Material()
eurofer_material.set_density('g/cm3', 0.01)
eurofer_material.add_element('Li', 1)
mats = openmc.Materials([breeder_material, eurofer_material, copper_material])
# GEOMETRY
# surfaces
central_sol_surface = openmc.ZCylinder(r=100)
central_shield_outer_surface = openmc.ZCylinder(r=110)
vessel_inner_surface = openmc.Sphere(r=500)
first_wall_outer_surface = openmc.Sphere(r=510)
breeder_blanket_outer_surface = openmc.Sphere(r=610, boundary_type='vacuum')
# cells
central_sol_region = -central_sol_surface & -breeder_blanket_outer_surface
central_sol_cell = openmc.Cell(region=central_sol_region)
central_sol_cell.fill = copper_material
central_shield_region = +central_sol_surface & -central_shield_outer_surface & -breeder_blanket_outer_surface
central_shield_cell = openmc.Cell(region=central_shield_region)
central_shield_cell.fill = eurofer_material
inner_vessel_region = -vessel_inner_surface & +central_shield_outer_surface
inner_vessel_cell = openmc.Cell(region=inner_vessel_region)
# no material set as default is vacuum
first_wall_region = -first_wall_outer_surface & +vessel_inner_surface
first_wall_cell = openmc.Cell(region=first_wall_region)
first_wall_cell.fill = eurofer_material
breeder_blanket_region = +first_wall_outer_surface & -breeder_blanket_outer_surface & +central_shield_outer_surface
breeder_blanket_cell = openmc.Cell(region=breeder_blanket_region)
breeder_blanket_cell.fill = breeder_material
universe = openmc.Universe(cells=[central_sol_cell, central_shield_cell, inner_vessel_cell, first_wall_cell, breeder_blanket_cell])
my_geometry = openmc.Geometry(universe)
# SIMULATION SETTINGS
# Instantiate a Settings object
sett = openmc.Settings()
batches = 2
sett.batches = batches
sett.inactive = 0
sett.particles = 50000
sett.run_mode = 'fixed source'
# Create a DT point source
# initialises a new source object
my_source = openmc.Source()
# the distribution of radius is just a single value
radius = openmc.stats.Discrete([130], [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)
sett.source = my_source
# Create mesh which will be used for tally
mesh = openmc.RegularMesh().from_domain(
my_geometry, # the corners of the mesh are being set automatically to surround the geometry
# dimension=[10, 20, 30] # voxels in each axis direction (x, y, z)
dimension=[10, 30, 20] # voxels in each axis direction (x, y, z)
)
tallies = openmc.Tallies()
# Create mesh filter for tally
mesh_filter = openmc.MeshFilter(mesh)
# Create flux mesh tally to score flux
mesh_tally_1 = openmc.Tally(name='tbr_on_mesh')
mesh_tally_1.filters = [mesh_filter]
mesh_tally_1.scores = ['heating'] # where X is a wildcard
tallies.append(mesh_tally_1)
model = openmc.model.Model(my_geometry, mats, sett, tallies)
sp_filename = model.run()
# loads up the output file from the simulation
statepoint = openmc.StatePoint(sp_filename)
# extracts the mesh tally by name
my_tbr_tally = statepoint.get_tally(name='tbr_on_mesh')
# converts the tally result into a VTK file
mesh.write_data_to_vtk(
filename="tbr_tally_on_reg_mesh.vtk",
datasets={"mean": my_tbr_tally.mean} # the first "mean" is the name of the data set label inside the vtk file
)
import matplotlib.pyplot as plt
import numpy as np
f=my_tbr_tally.mean
ff=f.flatten()
ff2 = ff.reshape(mesh.dimension, order='F')
ffx = ff2.transpose(2, 1, 0)
plt.imshow(ffx[10])
plt.show()
ffy = ff2.transpose(1, 2, 0)
plt.imshow(ffy[15])
plt.show()
ffz = ff2.transpose(0,2,1)
plt.imshow(ffz[5])
plt.show()
something like this would help