Any fission in a pebble fuel

[Fandresena]

Hello everyone, this is my project: i tried to simulate the fissions reactions inside pebbles fuels and after calculation, no fission was simulated. I’ll paste my code here, can anyone help me?

###TRISO Particle materials
coeur_comb=openmc.Material(name=‘Noyau combustible’)
coeur_comb.add_nuclide(‘U234’, 0.0015355, percent_type=‘ao’)
coeur_comb.add_nuclide(‘U235’, 0.171796, percent_type=‘ao’)
coeur_comb.add_nuclide(‘U236’, 0.0007869, percent_type=‘ao’)
coeur_comb.add_nuclide(‘U238’, 0.825882, percent_type=‘ao’)
coeur_comb.add_element(‘O’, 2.0)
coeur_comb.set_density(‘g/cm3’, 11.0)

PyC = openmc.Material(name =‘Buffer layer’)
PyC.add_element(‘C’, 1.0)
PyC.set_density(‘g/cm3’, 1.02)
PyC.add_s_alpha_beta(‘c_Graphite’)

IPyC = openmc.Material(name =‘Inner Pyrocarbon layer’)
IPyC.add_element(‘C’, 1.0)
IPyC.set_density(‘g/cm3’, 1.92)
IPyC.add_s_alpha_beta(‘c_Graphite’)
SiC = openmc.Material( name = ‘SiC layer’)

SiC.add_element(‘Si’, 1.0)
SiC.add_element(‘C’, 1.0)
SiC.set_density(‘g/cm3’, 3.2)
SiC.add_s_alpha_beta(‘c_Graphite’)

OPyC = openmc.Material(name =‘Outter Pyrocarbon layer’)
OPyC.add_element(‘C’, 1.0)
OPyC.set_density(‘g/cm3’, 1.92)
OPyC.add_s_alpha_beta(‘c_Graphite’)

###Reflector
graphite_refl= openmc.Material(name =‘Graphite’)
graphite_refl.add_element(‘C’, 1.0)
graphite_refl.set_density(‘g/cm3’, 1.8)

###Moderator
graphite_mod = openmc.Material(name =‘Graphite pour moderateur’)
graphite_mod.add_element(‘C’, 1.0)
graphite_mod.set_density(‘g/cm3’, 1.8)
graphite_mod.add_s_alpha_beta(‘c_Graphite’)

###Coolant
caloporteur = openmc.Material(name = ‘He_cal’)
caloporteur.add_element(‘He’, 1.0)
caloporteur.set_density(‘g/cm3’, 0.0001786)
materials = openmc.Materials ([coeur_comb, PyC, IPyC, SiC, OPyC, graphite_refl, graphite_mod, isolateur, caloporteur])
materials.export_to_xml()

**
###Creation of a TRISO particle**
s= openmc.Sphere(r=0.025)
s1= openmc.Sphere(r=0.035)
s2= openmc.Sphere(r=0.039)
s3=openmc.Sphere(r=0.0415)
s4=openmc.Sphere(r=0.0455)

kernel_zone = openmc.Cell(fill=coeur_comb, region=-s, name =‘kz’)
buff_zone = openmc.Cell(fill=PyC , region=+s & -s1, name =‘bz’)
IPyC_zone = openmc.Cell(fill=IPyC , region=+s1 & -s2, name =‘iz’)
SiC_zone = openmc.Cell(fill=SiC , region=+s2 & -s3, name =‘sz’)
OPyC_zone = openmc.Cell(fill=OPyC , region=+s3 & -s4, name =‘oz’)
cellules = [ kernel_zone, buff_zone, IPyC_zone, SiC_zone, OPyC_zone]
triso = openmc.Universe(cells = cellules)

#Creation of a pebble fuel surrounded by coolant
s6 = openmc.Sphere(r=3)
s5 = openmc.Sphere(r=2)

bille_1 = -s5
bille_2 = -s6 & +s5

couche_graph = openmc.Cell(fill = graphite_mod, region = bille_2, name =‘Couche_graphite2’)

rayon_triso= 0.0455
centers = openmc.model.pack_spheres(radius=rayon_triso, region=bille_1,pf = 0.13, seed = 11)

pebble_cells =
for c in centers:
pebble_cell = openmc.model.TRISO(rayon_triso, triso, center=c) # Universe comme fill
pebble_cells.append(pebble_cell)

zone_comb = openmc.Cell(region=bille_1)
lower_left, upper_right = zone_comb.region.bounding_box
shape = (5, 5, 5)
pitch = (upper_right - lower_left)/shape
lattice1 = openmc.model.create_triso_lattice(
pebble_cells, lower_left, pitch, shape, graphite_mod)

zone_comb.fill = lattice1

px1 = openmc.XPlane(x0=-3.5)
px2 = openmc.XPlane(x0=+3.5)
py1 = openmc.YPlane(y0=-3.5)
py2 = openmc.YPlane(y0=+3.5)
pz1 = openmc.ZPlane(z0=-3.5)
pz2 = openmc.ZPlane(z0=+3.5)

cell_box = +px1 & -px2 & +py1 & -py2 & +pz1 & -pz2
pebble_surf = openmc.Sphere(r=3.0)
coolant_cell = openmc.Cell(
fill=caloporteur,
region=cell_box & +pebble_surf
)
pebble = openmc.Universe(cells = [couche_graph, zone_comb, coolant_cell])
p1 = openmc.XPlane(x0 =10)
p2 = openmc.XPlane(x0 =-10)
p3 = openmc.YPlane(y0 =10)
p4 = openmc.YPlane(y0 =-10)
p5 = openmc.ZPlane(z0 =10)
p6 = openmc.ZPlane(z0 =-10)

cube_region = -p1 & +p2 & -p3 & +p4 & -p5 & +p6

#Creation of a lattice containing 3×3×3 pebbles

lattice_cell = openmc.Cell(region=cube_region)

lower_left, upper_right = lattice_cell.region.bounding_box
cube_size = upper_right - lower_left

nx, ny, nz = 3, 3, 3

pitch = cube_size / np.array([nx, ny, nz])

lattice2 = openmc.RectLattice()
lattice2.lower_left = lower_left
lattice2.pitch = pitch
lattice2.dimension = [nx, ny, nz]

lattice_universes =
for k in range(nz):
layer =
for j in range(ny):
row = [pebble] * nx
layer.append(row)
lattice_universes.append(layer)

lattice2.universes = lattice_universes

lattice_cell.fill = lattice2

#Creation of a graphite cube encompassing everything

p11 = openmc.XPlane(x0 =20, boundary_type =‘vacuum’)
p12 = openmc.XPlane(x0 =-20, boundary_type =‘vacuum’)
p13 = openmc.YPlane(y0 =20, boundary_type =‘vacuum’)
p14 = openmc.YPlane(y0 =-20, boundary_type =‘vacuum’)
p15 = openmc.ZPlane(z0 =20, boundary_type =‘vacuum’)
p17 = openmc.ZPlane(z0 =-20, boundary_type =‘vacuum’)

cube_region2 = ~cube_region & -p11 & +p12 & -p13 & +p14 & -p15 & +p17

graphite_cell = openmc.Cell(fill = graphite_mod, region = cube_region2)

global_universe = openmc.Universe(cells=[lattice_cell, graphite_cell])

#Mesh and Tallies

mesh = openmc.RegularMesh()
mesh.dimension = [50, 50, 50] # résolution spatiale
mesh.lower_left = [-20, -20, -20] # borne inférieure du domaine
mesh.upper_right = [20, 20, 20] # borne supérieure

mesh_filter = openmc.MeshFilter(mesh)

mesh_tally = openmc.Tally(name=‘mesh_flux_fission’)
mesh_tally.filters = [mesh_filter]
mesh_tally.scores = [‘flux’, ‘fission’]
tallies = openmc.Tallies()
tallies.append(mesh_tally)
tallies.export_to_xml()

#Settings
source = openmc.IndependentSource()
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Maxwell(1.0e6)

settings = openmc.Settings()
settings.batches = 200
settings.inactive = 50
settings.particles = 50000
settings.source = source
settings.export_to_xml()

[wahidluthfi]

Hi Fandresena, welcome to the openmc forum.
I have checked your input, and I think your fission reaction was not all zero eventought most of your tally result shows zero fission rate.
I think it was caused by your mesh definition, which covers all graphite regions with no fissile material, making the fission reaction rate shown as zero in those regions.

With your setup, the fission reaction rate was shown only for the mesh overlapping with the fueled pebble region

I am adding a small tally to show the fission reaction only on the pebble lattice 3x3x3 with finer mesh, so that you could see the meshtally result

notebook: pebblefission.ipynb (395.4 KB)

[Fandresena]

Thanks for you help, in another project based in HTR-PM reactor, i have the same problems: The neutron flux in the reactor core is equal to zero and reaches its maximum in the reflector, and the number of fissions is very low.
This project is in the following file.
Excuse me, can you help me to solve this problem too
Yours faithfully

24BC.ipynb (1.2 MB)

[wahidluthfi]

Hi Fandresena, sorry for the late reply.
Your model was so big that our typical computing power couldn’t handle it, hehehehe.

Regarding your model above, I have updated some parts of it, and tallying with the regular rectangular mesh, using 1,000,000 particles, 50 inactive, and 250 total batches to cover all regions within the core. The power distribution (kappa fission) and total neutron flux can be seen below

Here I normalize the power and flux with 1 kW th power, so the flux is low, but I hope you have some data for validating your core model, since I haven’t done this kind of random pebble bed core calculation before.
notebook: 24BC.ipynb (1.4 MB)

[Fandresena]

Hi, thanks for your help and.
I just finished my simulation based on the correction of the file you have sended. I would try the same way with the second one.
I am sincerely grateful for times and intentions that you are accordind to me.

[Fandresena]

Hi, sorry after this fews day, i try to insert the control rods in the models and after, i obtained the following neutron flux, it’s weird.

24BC.ipynb (699.4 KB)

[wahidluthfi]

Hi Fandresena,
I checked your geometry and did some updates to it to make sure that no cell overlappping each other. On the other hand, I also updated the lattice shape to shorten the running time, especially for the huge cell with a pebble within it.

Since I haven’t done a 3D mesh tally, I am only using a 2D mesh tally covering the whole core. I think you could try to use a proper index if you want to read the 3D mesh tally, Z,Y,X or Y,Z,X.

I also ran the geometry debug to make sure that each control rod didn’t overlap the reflector, but for the high-particle simulation, I turned it off.
24BC.ipynb (1.5 MB)

[Fandresena]

Hi, thanks for these correction. When i try to define a properly temperature for each cell, then when i start the simulation, max and mini temperature are in the following pictures but i’m giving many temperature like in the files:

24BC_with_temp.ipynb (1.1 MB)