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:

Cross-section data (mono-energetic, cm-1)#

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
 5simulation = mcdc.Simulation("Moving pellet")
 6
 7# ======================================================================================
 8# Set model
 9# ======================================================================================
10
11# Set materials
12fuel = mcdc.Material.multigroup(
13    capture=np.array([0.5]),
14    fission=np.array([0.25]),
15    nu_p=np.array([1.5]),
16    speed=np.array([200000.0]),
17)
18air = mcdc.Material.multigroup(
19    capture=np.array([0.002]),
20    scatter=np.array([[0.008]]),
21    speed=np.array([200000.0]),
22)

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)

24# Set surfaces
25cylinder_z = mcdc.Surface.CylinderZ(center=[0.0, 0.0], radius=1.0)
26top_z = mcdc.Surface.PlaneZ(z=9.0)
27bot_z = mcdc.Surface.PlaneZ(z=6.0)
28
29# Move surfaces
30cylinder_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])

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)

32bot_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])
33
34# Set container cell surfaces
35min_x = mcdc.Surface.PlaneX(x=-5.0, boundary_condition="vacuum")
36max_x = mcdc.Surface.PlaneX(x=5.0, boundary_condition="vacuum")
37min_y = mcdc.Surface.PlaneY(y=-5.0, boundary_condition="vacuum")
38max_y = mcdc.Surface.PlaneY(y=5.0, boundary_condition="vacuum")
39min_z = mcdc.Surface.PlaneZ(z=-10.0, boundary_condition="vacuum")
40max_z = mcdc.Surface.PlaneZ(z=10.0, boundary_condition="vacuum")
41
42# Make cells
43fuel_pellet_region = +bot_z & -top_z & -cylinder_z
44fuel_cell = mcdc.Cell(region=fuel_pellet_region, fill=fuel)
45air_cell = mcdc.Cell(
46    region=~fuel_pellet_region & +min_x & -max_x & +min_y & -max_y & +min_z & -max_z,
47    fill=air,
48)
49simulation.set_model([fuel_cell, air_cell])
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    x=[2.0, 3.0],
57    y=[-0.5, 0.5],
58    z=[-0.5, 0.5],
59    isotropic=True,
60    energy=0,
61    time=[0.0, 9.0],
62)
63simulation.set_sources([source])

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# Tallies
70mesh = mcdc.MeshStructured(
71    x=np.linspace(-5, 5, 201),
72    z=np.linspace(-10, 10, 201),
73)
74tally = mcdc.Tally(mesh=mesh, scores=["fission"], time=np.linspace(0, 9, 46))
75simulation.set_tallies([tally])
76
77# Settings
78simulation.settings.N_particle = 100000
79simulation.settings.N_batch = 2
80simulation.settings.active_bank_buffer = 1000
81
82# Run (or visualize)
83visualize = False

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 = True to watch the geometry evolve with simulation.visualize_model(..., time=...).

  • Compare with moving_source to 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
  5simulation = mcdc.Simulation("Moving pellet")
  6
  7# ======================================================================================
  8# Set model
  9# ======================================================================================
 10
 11# Set materials
 12fuel = mcdc.Material.multigroup(
 13    capture=np.array([0.5]),
 14    fission=np.array([0.25]),
 15    nu_p=np.array([1.5]),
 16    speed=np.array([200000.0]),
 17)
 18air = mcdc.Material.multigroup(
 19    capture=np.array([0.002]),
 20    scatter=np.array([[0.008]]),
 21    speed=np.array([200000.0]),
 22)
 23
 24# Set surfaces
 25cylinder_z = mcdc.Surface.CylinderZ(center=[0.0, 0.0], radius=1.0)
 26top_z = mcdc.Surface.PlaneZ(z=9.0)
 27bot_z = mcdc.Surface.PlaneZ(z=6.0)
 28
 29# Move surfaces
 30cylinder_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])
 31top_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])
 32bot_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])
 33
 34# Set container cell surfaces
 35min_x = mcdc.Surface.PlaneX(x=-5.0, boundary_condition="vacuum")
 36max_x = mcdc.Surface.PlaneX(x=5.0, boundary_condition="vacuum")
 37min_y = mcdc.Surface.PlaneY(y=-5.0, boundary_condition="vacuum")
 38max_y = mcdc.Surface.PlaneY(y=5.0, boundary_condition="vacuum")
 39min_z = mcdc.Surface.PlaneZ(z=-10.0, boundary_condition="vacuum")
 40max_z = mcdc.Surface.PlaneZ(z=10.0, boundary_condition="vacuum")
 41
 42# Make cells
 43fuel_pellet_region = +bot_z & -top_z & -cylinder_z
 44fuel_cell = mcdc.Cell(region=fuel_pellet_region, fill=fuel)
 45air_cell = mcdc.Cell(
 46    region=~fuel_pellet_region & +min_x & -max_x & +min_y & -max_y & +min_z & -max_z,
 47    fill=air,
 48)
 49simulation.set_model([fuel_cell, air_cell])
 50
 51# ======================================================================================
 52# Set source
 53# ======================================================================================
 54
 55source = mcdc.Source(
 56    x=[2.0, 3.0],
 57    y=[-0.5, 0.5],
 58    z=[-0.5, 0.5],
 59    isotropic=True,
 60    energy=0,
 61    time=[0.0, 9.0],
 62)
 63simulation.set_sources([source])
 64
 65# ======================================================================================
 66# Set tallies, settings, and run MC/DC
 67# ======================================================================================
 68
 69# Tallies
 70mesh = mcdc.MeshStructured(
 71    x=np.linspace(-5, 5, 201),
 72    z=np.linspace(-10, 10, 201),
 73)
 74tally = mcdc.Tally(mesh=mesh, scores=["fission"], time=np.linspace(0, 9, 46))
 75simulation.set_tallies([tally])
 76
 77# Settings
 78simulation.settings.N_particle = 100000
 79simulation.settings.N_batch = 2
 80simulation.settings.active_bank_buffer = 1000
 81
 82# Run (or visualize)
 83visualize = False
 84if not visualize:
 85    simulation.run()
 86else:
 87    colors = {
 88        fuel: "red",
 89        air: "blue",
 90    }
 91    simulation.visualize_model(
 92        vis_plane="xz",
 93        y=0.0,
 94        x=[-5.0, 5.0],
 95        z=[-10, 10],
 96        pixels=(100, 100),
 97        colors=colors,
 98        time=np.linspace(0.0, 9.0, 19),
 99        save_as="figure",
100    )

How to Run#

From inside examples/moving_pellet run:

python 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.