Time-Dependent Kobayashi Dog-Leg

Description

Time-dependent variant of the Kobayashi dog-leg shielding benchmark. See the steady-state Kobayashi example for geometry and material specifications; this variant drives the problem with a pulsed source and records time-resolved tallies.

Step-by-Step Walkthrough

This example reuses the exact geometry from the steady-state Kobayashi dog-leg benchmark. The key differences are highlighted below.

1. Source with a Time Window (line 54–60)

50mcdc.Source(
51    x=[0.0, 10.0],
52    y=[0.0, 10.0],
53    z=[0.0, 10.0],
54    isotropic=True,
55    energy_group=0,
56    time=[0.0, 50.0],
57)
58
59# ======================================================================================
60# Set tallies, settings, techniques, and run MC/DC

The source now has time=[0.0, 50.0], meaning particles are emitted over a 50 s window rather than instantaneously.

2. Time-Resolved Tallies (lines 66–69)

66mcdc.Tally(mesh=mesh, scores=["flux"], time=time_grid)
67mcdc.Tally(scores=["density"], time=time_grid)
68
69# Settings

A time grid is added to both the mesh tally and a global density tally. This creates a time-resolved \(\phi(x, y, t)\) dataset. Global Tally with scores=["density"] tracks total neutron population over time.

What to try:

  • Shorten the source time to create a short pulse and watch the neutron cloud propagate.

  • Add a finer time grid to capture early transient behaviour.

  • Compare with the steady-state Kobayashi results.

Full Input

Click here to view the input file: examples/kobayashi-TD/input.py.

The complete input used for this example is embedded below:

 1import numpy as np
 2import mcdc
 3
 4# ======================================================================================
 5# Set model
 6# ======================================================================================
 7# Based on Kobayashi dog-leg benchmark problem
 8# (PNE 2001, https://doi.org/10.1016/S0149-1970(01)00007-5)
 9
10# Set materials
11m = mcdc.MaterialMG(capture=np.array([0.05]), scatter=np.array([[0.05]]))
12m_void = mcdc.MaterialMG(capture=np.array([5e-5]), scatter=np.array([[5e-5]]))
13
14# Set surfaces
15sx1 = mcdc.Surface.PlaneX(x=0.0, boundary_condition="reflective")
16sx2 = mcdc.Surface.PlaneX(x=10.0)
17sx3 = mcdc.Surface.PlaneX(x=30.0)
18sx4 = mcdc.Surface.PlaneX(x=40.0)
19sx5 = mcdc.Surface.PlaneX(x=60.0, boundary_condition="vacuum")
20sy1 = mcdc.Surface.PlaneY(y=0.0, boundary_condition="reflective")
21sy2 = mcdc.Surface.PlaneY(y=10.0)
22sy3 = mcdc.Surface.PlaneY(y=50.0)
23sy4 = mcdc.Surface.PlaneY(y=60.0)
24sy5 = mcdc.Surface.PlaneY(y=100.0, boundary_condition="vacuum")
25sz1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="reflective")
26sz2 = mcdc.Surface.PlaneZ(z=10.0)
27sz3 = mcdc.Surface.PlaneZ(z=30.0)
28sz4 = mcdc.Surface.PlaneZ(z=40.0)
29sz5 = mcdc.Surface.PlaneZ(z=60.0, boundary_condition="vacuum")
30
31# Set cells
32# Source
33mcdc.Cell(region=+sx1 & -sx2 & +sy1 & -sy2 & +sz1 & -sz2, fill=m)
34# Voids
35channel_1 = +sx1 & -sx2 & +sy2 & -sy3 & +sz1 & -sz2
36channel_2 = +sx1 & -sx3 & +sy3 & -sy4 & +sz1 & -sz2
37channel_3 = +sx3 & -sx4 & +sy3 & -sy4 & +sz1 & -sz3
38channel_4 = +sx3 & -sx4 & +sy3 & -sy5 & +sz3 & -sz4
39void_channel = channel_1 | channel_2 | channel_3 | channel_4
40mcdc.Cell(region=void_channel, fill=m_void)
41# Shield
42box = +sx1 & -sx5 & +sy1 & -sy5 & +sz1 & -sz5
43mcdc.Cell(region=box & ~void_channel, fill=m)
44
45# ======================================================================================
46# Set source
47# ======================================================================================
48# The source pulses in t=[0,5]
49
50mcdc.Source(
51    x=[0.0, 10.0],
52    y=[0.0, 10.0],
53    z=[0.0, 10.0],
54    isotropic=True,
55    energy_group=0,
56    time=[0.0, 50.0],
57)
58
59# ======================================================================================
60# Set tallies, settings, techniques, and run MC/DC
61# ======================================================================================
62
63# Tallies
64time_grid = np.linspace(0.0, 200.0, 21)
65mesh = mcdc.MeshUniform(x=(0.0, 1.0, 60), y=(0.0, 1.0, 100))
66mcdc.Tally(mesh=mesh, scores=["flux"], time=time_grid)
67mcdc.Tally(scores=["density"], time=time_grid)
68
69# Settings
70mcdc.settings.N_particle = 100
71mcdc.settings.N_batch = 2
72
73# Techniques
74mcdc.simulation.implicit_capture()
75
76# Run
77mcdc.run()

How to Run

From the repository root run:

python examples/kobayashi-TD/input.py

Expected Output

Time-resolved mesh tally HDF5 file and an animation produced by process-output.py that visualises neutron density versus time.

References

See: Kobayashi et al. (2001), Progress in Nuclear Energy.