Kobayashi Dog-Leg Void Benchmark

Problem Description

A mono-energetic, three-dimensional shielding benchmark featuring a dog-leg vacuum channel embedded in a purely scattering/absorbing shield. The problem evaluates the ability of a Monte Carlo code to transport neutrons through deep-penetration streaming paths.

It is based on the NEA steady-state fixed-source benchmark problem suite by Kobayashi et al. [Kobayashi2001].

Geometry and Materials

The computational domain spans \(x \in [0,60]\), \(y \in [0,100]\), \(z \in [0,60]\) cm. Boundary conditions are reflective on the three symmetry planes (\(x=0\), \(y=0\), \(z=0\)) and vacuum elsewhere.

Three regions are defined:

  1. Source region — \(x \in [0,10]\), \(y \in [0,10]\), \(z \in [0,10]\) cm.

  2. Dog-leg void channel — an L-shaped duct connecting the source corner to the far side of the domain.

  3. Shield — the remaining volume.

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

Region

\(\Sigma_c\)

\(\Sigma_s\)

Shield

0.05

0.05

Void channel

\(5\times10^{-5}\)

\(5\times10^{-5}\)

Physical Assumptions

  • Mono-energetic (one-speed) neutron transport.

  • Isotropic scattering.

  • Steady-state (time-independent) fixed-source problem.

  • Implicit capture variance-reduction technique.

Numerical Setup

Spatial mesh (tally)

\(60 \times 100 \times 60\) uniform cells (1 cm spacing)

Tally score

Scalar flux

Source particles

\(10^{3}\) (demonstration; increase for production runs)

Batches

2

Quantities of Interest

  • Three-dimensional scalar flux distribution \(\phi(x,y,z)\).

  • Flux attenuation along the streaming channel and through the shield.

Reference Solution

Reference solutions are tabulated in [Kobayashi2001] for several axial slices. Post-processing in MC/DC generates \((x,y)\) flux maps at selected \(z\)-planes for direct comparison.

References

[Kobayashi2001] (1,2)

K. Kobayashi, N. Sugimura, and Y. Nagaya, “3D Radiation Transport Benchmark Problems and Results for Simple Geometries with Void Region,” Progress in Nuclear Energy, 39:2, 119–144 (2001). [link]

Step-by-Step Walkthrough

This section walks through the input file block by block.

1. Import and Materials (lines 1–13)

 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

Two mono-energetic multi-group materials are created: m for the shield (\(\Sigma_c = \Sigma_s = 0.05\)) and m_void for the dog-leg channel (\(10^{-4}\) total).

2. Surfaces (lines 15–30)

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

Fifteen planar surfaces define the 3-D bounding box and the internal partitions. Reflective conditions on sx1, sy1, sz1 exploit the quarter-symmetry; vacuum on the outer faces allows leakage.

3. Cells — CSG Region Definitions (lines 32–44)

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

Three cells cover the domain:

  • The source cell (a small corner cube) filled with shield material.

  • The void channel — four rectangular segments combined with the | (union) operator to form the L-shaped duct.

  • The shield — the full box minus the void channel, using the ~ (complement) operator.

4. Source (lines 50–57)

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)
57

An isotropic, uniformly distributed source fills the \(10 \times 10 \times 10\) cm corner cube.

5. Tallies, Settings, Techniques, and Run (lines 63–74)

63mesh = mcdc.MeshUniform(x=(0.0, 1.0, 60), y=(0.0, 1.0, 100), z=(0.0, 1.0, 60))
64mcdc.Tally(mesh=mesh, scores=["flux"])
65
66# Settings
67mcdc.settings.N_particle = 1000
68mcdc.settings.N_batch = 2
69
70# Techniques
71mcdc.simulation.implicit_capture()
72
73# Run
74mcdc.run()

  • A uniform \(60 \times 100 \times 60\) mesh tally records scalar flux.

  • 1 000 source particles in 2 batches (increase for production).

  • Implicit capture prevents particles from being absorbed prematurely.

  • mcdc.run() launches the simulation.

What to try:

  • Increase N_particle to \(10^5\) or more for smoother flux maps.

  • Change void-channel cross sections to see how attenuation changes.

  • Add a time grid to the tally for a transient variant (see the TD example).

Full Input

Click here to view the input file: examples/kobayashi/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)
57
58# ======================================================================================
59# Set tallies, settings, techniques, and run MC/DC
60# ======================================================================================
61
62# Tallies
63mesh = mcdc.MeshUniform(x=(0.0, 1.0, 60), y=(0.0, 1.0, 100), z=(0.0, 1.0, 60))
64mcdc.Tally(mesh=mesh, scores=["flux"])
65
66# Settings
67mcdc.settings.N_particle = 1000
68mcdc.settings.N_batch = 2
69
70# Techniques
71mcdc.simulation.implicit_capture()
72
73# Run
74mcdc.run()

How to Run

From the repository root run:

python examples/kobayashi/input.py

Expected Output

A mesh tally HDF5 file with 3-D flux data and example plotting using the companion process-output.py in the same directory.