Good evening everyone, i try to simulate a simplified model of HTR-10 reactor to obtain the neutron flux density in the core. And The neutron flux in the fuel region is very low and does not produce any fission. How can i fix this problem??
“The code and the graph obtained after the calculation are here:
import openmc
import openmc.model
import numpy as np
import os
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from scipy.interpolate import RegularGridInterpolator
###bibliotheque
os.environ[‘OPENMC_CROSS_SECTIONS’] = ‘/home/fandresena/bibliotheque/openmc/endfb80/endfb-viii.0-hdf5/cross_sections.xml’
###Creation de materiel
#Les elements pour les combustibles TRISO
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’, 10.5)
coeur_comb.color = ‘blue’
PyC = openmc.Material(name =‘Buffer layer’)
PyC.add_element(‘C’, 1.0)
PyC.set_density(‘g/cm3’, 1.02)
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)
OPyC = openmc.Material(name =‘Outter Pyrocarbon layer’)
OPyC.add_element(‘C’, 1.0)
OPyC.set_density(‘g/cm3’, 1.92)
#Reflecteur
graphite_refl= openmc.Material(name =‘Graphite’)
graphite_refl.add_element(‘C’, 1.0)
graphite_refl.set_density(‘g/cm3’, 1.8)
graphite_refl.add_s_alpha_beta(‘c_Graphite’)
graphite_refl.color = ()
#Moderateur
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’)
#Couche d’isolation
isolateur = openmc.Material(name = ‘isolateur’)
isolateur.add_element(‘C’, 0.95, ‘ao’)
isolateur.add_element(‘B’, 0.05, ‘ao’)
isolateur.set_density (‘g/cm3’, 1.8)
#Caloporteur
caloporteur = openmc.Material(name = ‘He_cal’)
caloporteur.add_element(‘He’, 1.0)
caloporteur.set_density(‘g/cm3’, 0.0001786)
#exportation vers xml
materials = openmc.Materials ([coeur_comb, PyC, IPyC, SiC, OPyC, graphite_refl, graphite_mod, isolateur, caloporteur])
materials.export_to_xml()
creation de la geometrie
#triso
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)
triso_p = openmc.Plot()
triso_p.universe = triso
triso_p.filename = “trisoParticle”
triso_p.pixels = (480, 480)
triso_p.width =(20, 20)
triso_p.color_by = “material”
triso_p.basis = ‘xz’
#creation d’un pebble
s6 = openmc.Sphere(r=3)
s5 = openmc.Sphere(r=2.5)
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.1, 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
lattice = openmc.model.create_triso_lattice(pebble_cells, lower_left, pitch, shape, graphite_mod)
zone_comb.fill = lattice
pebble = openmc.Universe(cells = [couche_graph, zone_comb])
pebble_p = openmc.Plot()
pebble_p.universe = pebble
pebble_p.basis = ‘xz’
pebble_p.filename = “pebble”
pebble_p.pixels = (1200, 1200)
pebble_p.width =(20, 20)
pebble_p.color_by = “material”
#Creation du coeur
gcyl = openmc.ZCylinder(r = 90)
dcyl = openmc.ZCylinder(r = 25)
plan_sup = openmc.ZPlane(z0=197)
plan_milieu1 = openmc.ZPlane(z0=50)
plan_milieu2 = openmc.ZPlane(z0=0.0)
plan_inf = openmc.ZPlane(z0=-197)
Gcyl = -gcyl & -plan_sup & +plan_milieu1
Dech = -plan_milieu2 & -dcyl & +plan_inf
rayon_pebble = 3
centers = openmc.model.pack_spheres(radius=rayon_pebble, region=Gcyl, pf=0.5, seed = 12)
fuel_cells =
for c in centers:
fuel_cell = openmc.model.TRISO(rayon_pebble, pebble, center=c) # Universe comme fill
fuel_cells.append(fuel_cell)
Coeur_sup = openmc.Cell(region=Gcyl)
lower_left = np.array((-93,-93,47))
upper_right = np.array((93,93,200))
shape = np.array((20, 20, 20))
pitch = (upper_right - lower_left)/shape
Gcyli = openmc.model.create_triso_lattice(
fuel_cells, lower_left, pitch, shape, caloporteur)
Coeur_sup.fill =Gcyli
centers = openmc.model.pack_spheres(radius=rayon_pebble, region=Dech, pf=0.5, seed = 13)
fuel_cells =
for c in centers:
fuel_cell = openmc.model.TRISO(rayon_pebble, pebble, center=c) # Universe comme fill
fuel_cells.append(fuel_cell)
Coeur_inf = openmc.Cell(region=Dech)
lower_left = np.array((-27,-27,-197))
upper_right = np.array((27,27,2))
shape = np.array((10, 10, 10))
pitch = (upper_right - lower_left)/shape
Decha= openmc.model.create_triso_lattice(
fuel_cells, lower_left, pitch, shape, caloporteur)
Coeur_inf.fill =Decha
c1 = openmc.ZCylinder(r = 77)
p1 = openmc.ZPlane(z0= 37.5)
c2 = openmc.ZCylinder(r = 64)
p2 = openmc.ZPlane(z0= 25)
c3 = openmc.ZCylinder(r = 51)
p3 = openmc.ZPlane(z0= 12.5)
c4 = openmc.ZCylinder(r = 38)
co1= -plan_milieu1 & -c1 & +p1
co2= -p1 & -c2 & +p2
co3= -p2 & -c3 & +p3
co4= -p3 & -c4 & +plan_milieu2
centers = openmc.model.pack_spheres(radius=rayon_pebble, region=co1, pf=0.4, seed = 15)
fuel_cells =
for c in centers:
fuel_cell = openmc.model.TRISO(rayon_pebble, pebble, center=c) # Universe comme fill
fuel_cells.append(fuel_cell)
ent1 = openmc.Cell(region=co1)
lower_left = np.array((-85,-85,30))
upper_right = np.array((85, 85, 60))
shape = np.array((10, 10, 10))
pitch = (upper_right - lower_left)/shape
entc= openmc.model.create_triso_lattice(
fuel_cells, lower_left, pitch, shape, caloporteur)
ent1.fill = entc
centers = openmc.model.pack_spheres(radius=rayon_pebble, region=co2, pf=0.3, seed = 16)
fuel_cells =
for c in centers:
fuel_cell = openmc.model.TRISO(rayon_pebble, pebble, center=c) # Universe comme fill
fuel_cells.append(fuel_cell)
ent2 = openmc.Cell(region=co2)
lower_left = np.array((-67,-67,23))
upper_right = np.array((67, 67, 40))
shape = np.array((10, 10, 10))
pitch = (upper_right - lower_left)/shape
entc= openmc.model.create_triso_lattice(
fuel_cells, lower_left, pitch, shape, caloporteur)
ent2.fill = entc
centers = openmc.model.pack_spheres(radius=rayon_pebble, region=co3, pf=0.3, seed =17)
fuel_cells =
for c in centers:
fuel_cell = openmc.model.TRISO(rayon_pebble, pebble, center=c) # Universe comme fill
fuel_cells.append(fuel_cell)
ent3 = openmc.Cell(region=co3)
lower_left = np.array((-60,-60,10))
upper_right = np.array((60, 60, 35))
shape = np.array((10, 10, 10))
pitch = (upper_right - lower_left)/shape
entc= openmc.model.create_triso_lattice(
fuel_cells, lower_left, pitch, shape, caloporteur)
ent3.fill = entc
centers = openmc.model.pack_spheres(radius=rayon_pebble, region=co4, pf=0.3, seed =18)
fuel_cells =
for c in centers:
fuel_cell = openmc.model.TRISO(rayon_pebble, pebble, center=c) # Universe comme fill
fuel_cells.append(fuel_cell)
ent4 = openmc.Cell(region=co4)
lower_left = np.array((-45, -45, -4))
upper_right = np.array((45, 45, 25))
shape = np.array((10, 10, 10))
pitch = (upper_right - lower_left)/shape
entc= openmc.model.create_triso_lattice(
fuel_cells, lower_left, pitch, shape, caloporteur)
ent4.fill = entc
#creation du cuve
plan_dessus1 = openmc.ZPlane(z0 =257,boundary_type =‘vacuum’)
plan_dessus2 = openmc.ZPlane(z0 =247)
plan_bas1 = openmc.ZPlane(z0 =-257, boundary_type =‘vacuum’)
plan_bas2 = openmc.ZPlane(z0 =-247)
cyl_ref = openmc.ZCylinder(r = 195)
cyl_lim = openmc.ZCylinder(r = 205, boundary_type =‘vacuum’)
refl_reg = -plan_dessus2 & -cyl_ref & +plan_bas2 & ~Gcyl
isol_reg = -plan_dessus1 & -cyl_lim & +plan_bas1 & ~Gcyl
reflect = openmc.Cell(fill =graphite_refl, region =refl_reg)
isolat = openmc.Cell(fill =isolateur, region =isol_reg)
cuve = openmc.Universe(cells =[reflect, isolat, ent1, ent2, ent3, ent4, Coeur_inf, Coeur_sup])
cuve_p = openmc.Plot()
cuve_p.universe = cuve
cuve_p.filename = ‘cuve_du reacteur’
cuve_p.width =(700,700)
cuve_p.pixels = (1200,1200)
cuve_p.basis = ‘xz’
cuve.colors_by =‘material’
geometry =openmc.Geometry(cuve)
geometry.export_to_xml()
plots = openmc.Plots([triso_p, pebble_p,cuve_p])
plots.export_to_xml()
openmc.plot_geometry()
tally
#maillage
mesh = openmc.CylindricalMesh(
r_grid=np.linspace(0, 205, 21),
z_grid=np.linspace(-257, 257, 157),
phi_grid=np.linspace(0, 2*np.pi, 50),
)
mesh_filter = openmc.MeshFilter(mesh)
tally = openmc.Tally(name=‘flux’)
tally.filters = [mesh_filter]
tally.scores = [‘flux’]
tallies = openmc.Tallies([tally])
tallies.export_to_xml()
settings = openmc.Settings()
settings.batches = 250
settings.inactive = 50
settings.particles = 5000
source = openmc.IndependentSource(space=openmc.stats.Point((0.0, 0.0, 90.0)))
settings.source = source
settings.export_to_xml()
openmc.run()
sp = openmc.StatePoint(‘statepoint.250.h5’)
tally = sp.get_tally(scores=[‘flux’])
print(tally)
flux = tally.get_slice(scores=[‘flux’])
print(flux)
reshape du flux en 3D (r, z, phi)
flux_3d = flux.mean.reshape(20, 156, 49)
grilles du mesh
r = np.linspace(0, 205, 20)
z = np.linspace(-257, 257, 156)
phi = np.linspace(0, 2*np.pi, 49)
créer un interpolateur en cylindrique
interp = RegularGridInterpolator(
(r, z, phi), flux_3d, bounds_error=False, fill_value=None
)
choisir une coupe en z fixe
z_fixed = 100 # ← tu mets la valeur que tu veux
z_line = None # juste pour info
paramètres de la ligne
R_max = r[-1]
num_points = 200
côté gauche : phi = pi
r_left = np.linspace(R_max, 0, num_points)
phi_left = np.full_like(r_left, np.pi)
côté droit : phi = 0
r_right = np.linspace(0, R_max, num_points)
phi_right = np.zeros_like(r_right)
concaténer gauche → centre → droite
r_line = np.concatenate([r_left, r_right])
phi_line = np.concatenate([phi_left, phi_right])
z_line = np.full_like(r_line, z_fixed)
interpolation
points = np.vstack([r_line, z_line, phi_line]).T
flux_line = interp(points)
créer un axe x équivalent (sans coordonnées cartésiennes)
x_line = np.concatenate([-r_left, r_right]) # gauche = négatif, droite = positif
tracer
plt.figure(figsize=(8,5))
plt.plot(x_line, flux_line, ‘-b’)
plt.xlabel(“Position x équivalente (cm)”)
plt.ylabel(“Flux moyen”)
plt.title(f"Variation du flux le long de r (z = {z_fixed}) en cylindrique")
plt.grid(True)
plt.savefig(“flux_line_z100.png”, dpi=300, bbox_inches=“tight”)
plt.close()
iz = np.argmin(np.abs(z - 50))
flux_r_core = flux_3d[:, iz, :].mean(axis=1)
plt.figure(figsize=(7,4))
plt.plot(r, flux_r_core, ‘o-’)
plt.axvline(90, color=‘k’, linestyle=‘–’, label=‘Limite cœur’)
plt.axvline(195, color=‘r’, linestyle=‘–’, label=‘Limite réflecteur’)
plt.xlabel(“r (cm)”)
plt.ylabel(“Flux moyen”)
plt.title(“Profil radial moyen du flux à z = 50 cm”)
plt.legend()
plt.grid()
plt.savefig(“flux_radial_z50.png”, dpi=300, bbox_inches=“tight”)
plt.close()
neutrons die out quickly → very low flux → no fission multiplication.