Moving Fuel Pellet¶
Problem Description¶
A three-dimensional time-dependent problem in which a cylindrical fuel pellet traverses a box of air-like material along a piecewise-linear trajectory. Both the cylindrical surface and the bounding planes of the pellet move, demonstrating MC/DC’s moving-surface capability for transient geometry.
Geometry and Materials¶
The computational domain is a rectangular box: \(x \in [-5,5]\), \(y \in [-5,5]\), \(z \in [-10,10]\) cm, with vacuum boundary conditions on all faces.
A z-aligned cylindrical fuel pellet (radius 1 cm) is initially located at \(z \in [6,9]\) cm. The pellet moves in three consecutive phases:
Phase |
Cylinder velocity |
End-cap velocity |
Duration (s) |
1 |
\((-0.5,\; 0,\; 0)\) |
\((0,\; 0,\; -2)\) |
2 / 5 |
2 |
\((1,\; 0,\; 0)\) |
\((0,\; 0,\; 4)\) |
5 / 2 |
3 |
\((-2,\; 0,\; 0)\) |
\((0,\; 0,\; -10)\) |
1 / 1 |
Two materials are used:
Region |
\(\Sigma_c\) |
\(\Sigma_s\) |
\(\Sigma_f\) |
\(\nu\) |
Fuel pellet |
0.50 |
— |
0.25 |
1.5 |
Air |
0.002 |
0.008 |
— |
— |
Neutron speed: \(v = 2 \times 10^{5}\) cm/s (both regions).
Physical Assumptions¶
Mono-energetic (one-speed) neutron transport.
Isotropic scattering.
Time-dependent transport with moving geometry surfaces.
Fission in the pellet region only.
Numerical Setup¶
Spatial mesh (tally) |
\(201 \times 201\) in the \((x,z)\)-plane |
Time mesh (tally) |
46 equally spaced bins over \(t \in [0,9]\) s |
Tally score |
Fission rate |
Source particles |
\(10^{5}\) |
Batches |
2 |
Quantities of Interest¶
Time-resolved 2-D fission rate distribution in the \((x,z)\)-plane.
Animation of the fission rate tracking the moving pellet geometry.
Reference Solution¶
No analytical reference. The solution is validated by verifying that the fission rate follows the pellet trajectory and that particle conservation is maintained.
Step-by-Step Walkthrough¶
1. Materials (lines 1–22)
1import numpy as np
2
3import mcdc
4
5# ======================================================================================
6# Set model
7# ======================================================================================
8
9# Set materials
10fuel = mcdc.MaterialMG(
11 capture=np.array([0.5]),
12 fission=np.array([0.25]),
13 nu_p=np.array([1.5]),
14 speed=np.array([200000.0]),
15)
16air = mcdc.MaterialMG(
17 capture=np.array([0.002]),
18 scatter=np.array([[0.008]]),
19 speed=np.array([200000.0]),
20)
21
22# Set surfaces
A fissile fuel pellet (\(\Sigma_f = 0.25\), \(\nu = 1.5\)) and
an air-like background. Both include speed for time-dependent
transport.
2. Surfaces and Moving Geometry (lines 24–30)
24top_z = mcdc.Surface.PlaneZ(z=9.0)
25bot_z = mcdc.Surface.PlaneZ(z=6.0)
26
27# Move surfaces
28cylinder_z.move([[-0.5, 0.0, 0.0], [1.0, 0.0, 0.0], [-2.0, 0.0, 0.0]], [2.0, 5.0, 1.0])
29top_z.move([[0.0, 0.0, -2.0], [0.0, 0.0, 4.0], [0.0, 0.0, -10.0]], [5.0, 2.0, 1.0])
30bot_z.move([[0.0, 0.0, -2.0], [0.0, 0.0, 4.0], [0.0, 0.0, -10.0]], [5.0, 2.0, 1.0])
A z-cylinder and two z-planes define the pellet. The key feature:
surface.move(velocities, durations) makes these surfaces time-dependent.
The cylinder moves laterally while the endcaps move axially, simulating
a pellet traversing the domain.
3. Container and Cells (lines 32–50)
32# Set container cell surfaces
33min_x = mcdc.Surface.PlaneX(x=-5.0, boundary_condition="vacuum")
34max_x = mcdc.Surface.PlaneX(x=5.0, boundary_condition="vacuum")
35min_y = mcdc.Surface.PlaneY(y=-5.0, boundary_condition="vacuum")
36max_y = mcdc.Surface.PlaneY(y=5.0, boundary_condition="vacuum")
37min_z = mcdc.Surface.PlaneZ(z=-10.0, boundary_condition="vacuum")
38max_z = mcdc.Surface.PlaneZ(z=10.0, boundary_condition="vacuum")
39
40# Make cells
41fuel_pellet_region = +bot_z & -top_z & -cylinder_z
42mcdc.Cell(region=fuel_pellet_region, fill=fuel)
43mcdc.Cell(
44 region=~fuel_pellet_region & +min_x & -max_x & +min_y & -max_y & +min_z & -max_z,
45 fill=air,
46)
47
48# ======================================================================================
49# Set source
50# ======================================================================================
The fuel pellet region is defined by the intersection of the cylinder and the two planes. The air fills the complement inside the bounding box.
4. Source (lines 56–63)
56 isotropic=True,
57 energy_group=0,
58 time=[0.0, 9.0],
59)
60
61# ======================================================================================
62# Set tallies, settings, and run MC/DC
63# ======================================================================================
A small box source near the pellet’s initial position, active over the full simulation time \(t \in [0, 9]\) s.
5. Tallies, Settings, and Run (lines 69–83)
69)
70mcdc.Tally(mesh=mesh, scores=["fission"], time=np.linspace(0, 9, 46))
71
72# Settings
73mcdc.settings.N_particle = 100000
74mcdc.settings.N_batch = 2
75mcdc.settings.active_bank_buffer = 1000
76
77# Run (or visualize)
78visualize = False
79if not visualize:
80 mcdc.run()
81else:
82 colors = {
83 fuel: "red",
A structured mesh tally in the \((x,z)\)-plane with 46 time bins captures the fission rate as the pellet moves.
What to try:
Change the pellet velocities to create different trajectories.
Set
visualize = Trueto watch the geometry evolve withmcdc.visualize(..., time=...).Compare with
moving_sourceto see source motion vs. geometry motion.
Full Input¶
Click here to view the input file: examples/moving_pellet/input.py.
The complete input used for this example is embedded below:
1import numpy as np
2
3import mcdc
4
5# ======================================================================================
6# Set model
7# ======================================================================================
8
9# Set materials
10fuel = mcdc.MaterialMG(
11 capture=np.array([0.5]),
12 fission=np.array([0.25]),
13 nu_p=np.array([1.5]),
14 speed=np.array([200000.0]),
15)
16air = mcdc.MaterialMG(
17 capture=np.array([0.002]),
18 scatter=np.array([[0.008]]),
19 speed=np.array([200000.0]),
20)
21
22# Set surfaces
23cylinder_z = mcdc.Surface.CylinderZ(center=[0.0, 0.0], radius=1.0)
24top_z = mcdc.Surface.PlaneZ(z=9.0)
25bot_z = mcdc.Surface.PlaneZ(z=6.0)
26
27# Move surfaces
28cylinder_z.move([[-0.5, 0.0, 0.0], [1.0, 0.0, 0.0], [-2.0, 0.0, 0.0]], [2.0, 5.0, 1.0])
29top_z.move([[0.0, 0.0, -2.0], [0.0, 0.0, 4.0], [0.0, 0.0, -10.0]], [5.0, 2.0, 1.0])
30bot_z.move([[0.0, 0.0, -2.0], [0.0, 0.0, 4.0], [0.0, 0.0, -10.0]], [5.0, 2.0, 1.0])
31
32# Set container cell surfaces
33min_x = mcdc.Surface.PlaneX(x=-5.0, boundary_condition="vacuum")
34max_x = mcdc.Surface.PlaneX(x=5.0, boundary_condition="vacuum")
35min_y = mcdc.Surface.PlaneY(y=-5.0, boundary_condition="vacuum")
36max_y = mcdc.Surface.PlaneY(y=5.0, boundary_condition="vacuum")
37min_z = mcdc.Surface.PlaneZ(z=-10.0, boundary_condition="vacuum")
38max_z = mcdc.Surface.PlaneZ(z=10.0, boundary_condition="vacuum")
39
40# Make cells
41fuel_pellet_region = +bot_z & -top_z & -cylinder_z
42mcdc.Cell(region=fuel_pellet_region, fill=fuel)
43mcdc.Cell(
44 region=~fuel_pellet_region & +min_x & -max_x & +min_y & -max_y & +min_z & -max_z,
45 fill=air,
46)
47
48# ======================================================================================
49# Set source
50# ======================================================================================
51
52mcdc.Source(
53 x=[2.0, 3.0],
54 y=[-0.5, 0.5],
55 z=[-0.5, 0.5],
56 isotropic=True,
57 energy_group=0,
58 time=[0.0, 9.0],
59)
60
61# ======================================================================================
62# Set tallies, settings, and run MC/DC
63# ======================================================================================
64
65# Tallies
66mesh = mcdc.MeshStructured(
67 x=np.linspace(-5, 5, 201),
68 z=np.linspace(-10, 10, 201),
69)
70mcdc.Tally(mesh=mesh, scores=["fission"], time=np.linspace(0, 9, 46))
71
72# Settings
73mcdc.settings.N_particle = 100000
74mcdc.settings.N_batch = 2
75mcdc.settings.active_bank_buffer = 1000
76
77# Run (or visualize)
78visualize = False
79if not visualize:
80 mcdc.run()
81else:
82 colors = {
83 fuel: "red",
84 air: "blue",
85 }
86 mcdc.visualize(
87 "xz",
88 y=0.0,
89 x=[-5.0, 5.0],
90 z=[-10, 10],
91 pixels=(100, 100),
92 colors=colors,
93 time=np.linspace(0.0, 9.0, 19),
94 save_as="figure",
95 )
How to Run¶
From the repository root run:
python examples/moving_pellet/input.py
Expected Output¶
An HDF5 tally file with time-resolved fission rates and an optional animation created by the example’s post-processing script.