Tutorial: Shielding Analysis

Problem Setup

A 14 MeV D-T neutron point source sits at the origin behind a composite shield: 1.5 cm of steel followed by 30 cm of ordinary concrete. Detector points at 50 cm and 100 cm measure the attenuated flux beyond the shield.

The geometry is deliberately spare. A void sphere of radius 0.5 cm holds the source, so neutrons leave it without scattering and enter the steel at their birth energy. Steel fills the space out to the plane at x = 2.0 cm, which is 1.5 cm of thickness along the detector axis, and ordinary concrete carries from there to x = 32.0 cm — thirty centimeters, the sort of biological shield a small accelerator facility would actually have. A plane at x = −2.0 cm closes the steel behind the source so that it does not extend to the outer boundary, and a sphere of radius 200 cm terminates the problem.

Four questions make this worth running rather than estimating. How much does each layer attenuate, which is what the surface tallies answer. What flux survives out at the detector points, which needs a point detector because a cell tally 50 cm out in air would score almost nothing. How the spectrum changes shape, since a shield does not simply reduce a flux — it removes fast neutrons while creating thermal ones. And whether an importance ladder is enough variance reduction, or whether the deep bins need more.

Complete Input File

Hover over any highlighted section to see what it does and why.

Neutron Shielding Analysis Example
c Cell Cards
c Geometry: point source -> steel slab -> concrete slab -> air,
c all inside a large bounding sphere.
1  0         -1        imp:n=1   $ Source region (inside small sphere)
2  2  -7.85   1 4 -2 -5  imp:n=2 $ Steel slab, x = -2 to 2 cm
3  3  -2.3    2 -3 -5    imp:n=4 $ Concrete slab (x = 2 to 32 cm)
4  0          3 -5      imp:n=2   $ Air / detector region
6  0         -4 -5      imp:n=1   $ Air behind the source
5  0          5         imp:n=0   $ Outside world
 
c Surface Cards
1  so   0.5      $ Source sphere (at origin)
2  px   2.0      $ Steel / concrete interface
3  px   32.0     $ Concrete back face
4  px  -2.0      $ Steel front face (behind the source)
5  so   200.0    $ Outer boundary sphere
 
c Data Cards
c Materials
m2  26056.70c  -1.0       $ Steel (simplified, single isotope)
m3  1001.70c   -0.01      $ Portland concrete
    8016.70c   -0.532
    14028.70c  -0.337
    20040.70c  -0.044
    26056.70c  -0.014
    13027.70c  -0.034
    11023.70c  -0.029
 
c Source - 14 MeV neutron point source
sdef  par=n  erg=14  pos=0 0 0   $ par=particle, erg=energy(MeV)
 
c Tallies
f2:n  2 3                   $ Surface flux at surfaces 2 and 3
e2    1e-9 1e-6 1e-3 0.1 1 5 10 15   $ Energy bin upper bounds (MeV)
f5:n  50 0 0  0             $ Point detector at (50,0,0), R=0
      100 0 0  0            $ Point detector at (100,0,0), R=0
e5    1e-9 1e-6 1e-3 0.1 1 5 10 15
f4:n  2 3                   $ Volume-averaged flux in cells 2 and 3
 
c Physics and cutoffs
mode  n                     $ Neutron transport only
c No cut:n energy cutoff - keeps the lowest energy bins (see notes)
nps   1e6                   $ Number of source particles

Variance Reduction

The importance ladder (imp:n=1 → 2 → 4 → 2 → 0) built into the cell cards above is already a basic form of variance reduction. For very deep penetration problems (shield thickness of 10+ mean free paths) you may need additional techniques. A DXTRAN sphere around the detector region forces neutrons to scatter toward the detector.

c Add to the previous input:
c Variance reduction
c DXTRAN sphere around detector region (not at source)
c Syntax: DXT:n x y z ri ro
dxt:n   100 0 0  5.0  40.0  $ At x=100, inner=5 outer=40 cm
c Optional. 1e-8 MeV is 0.01 eV, so the thermal peak at 0.025 eV
c survives, but the lowest e2/e5 bin is zeroed. Drop this line if
c you care about the bottom of the spectrum, as the main deck does.
cut:n   J 1e-8              $ J=default time, 0.01 eV energy cutoff
nps     5e5

Running the Calculation

# Run the shielding calculation
mcnp6 i=shield_input o=shield_output

# Monitor progress
tail -f shield_output

# Check for completion
grep "mcnp     version" shield_output

Four things come out of the output file, and they are best read in this order. The ratio of the F2 tally on surface 2 to the one on surface 3 is the attenuation of the concrete, bin by bin. The F5 results give the flux at 50 and 100 cm, which is what a dose calculation would start from. The F4 spectrum shows what the shield did to the shape of the flux rather than merely to its size. And the statistical checks decide whether any of the first three mean anything: relative errors below 0.10, all ten checks passed, and the errors examined bin by bin rather than on the total, because the total is dominated by the fast bins that converge easily.

Typical Results

Work out what you should see before you run, because the surprise is instructive either way. The 1.5 cm of steel on the detector axis attenuates the 14 MeV flux modestly and degrades it through inelastic scattering and (n,2n), which is a small effect compared with what follows. Through the concrete, fast flux falls roughly exponentially, and the relaxation length for fast neutrons in ordinary concrete is on the order of 10 to 15 cm. Thirty centimeters is therefore two to three e-foldings, or a factor of about 8 to 20 on the fast group as a whole — considerably less than the orders of magnitude people expect, because neutrons scattering down from higher energies keep refilling the group behind the ones being removed.

The topmost energy bin behaves differently, and the contrast is the lesson. A neutron still at 14 MeV after 30 cm of concrete is one that has not interacted at all, so it attenuates on the total cross section rather than the removal cross section — about 0.12 cm⁻¹ for this concrete, which gives nearer a factor of 40 across the slab, and more again at the point detectors where 1/r² also applies. A large ratio in the 10 to 15 MeV bin is expected. A large ratio in the group total is a reason to go back and check the geometry.

Meanwhile the thermal flux rises inside the concrete even as the fast flux falls, because the hydrogen in it is moderating rather than merely removing. The low bins on e2should grow between surface 2 and surface 3 while the top bins shrink. Treat that shape as qualitative, though: m3 uses free-gas hydrogen with no S(α,β) card, so the deck cannot place the thermal peak accurately, and adding mt3 lwtr.20t is the first thing to do if you want to trust those bins. The two lowest bins will also be poorly converged without more variance reduction, since very few source neutrons thermalize and then find their way to a point detector.

Where to take it next

Three geometry changes each teach something this deck cannot. An air gap between the steel and the concrete shows how weakly a void attenuates and how much streaming it permits. Rebar inside the concrete introduces a heterogeneity that shifts the spectrum locally. And a duct or streaming path through the shield is the case where every intuition from slab attenuation fails — the flux downstream of a penetration is dominated by scattering off the duct walls, and it is where shielding calculations most often go wrong in practice.

On the physics side, the obvious next step is dose rather than flux. That means turning on photon transport with mode n p, since gamma rays produced by neutron capture in the shield frequently dominate the dose on the far side, and applying flux-to-dose conversion factors through a DE/DF pair. Do not hand-type those tables: the real neutron response spans a factor of about 55 across the energy range, so a table that looks nearly flat is one that was invented. Adding S(α,β) for the concrete hydrogen and replacing the monoenergetic source with a fission spectrum both move the model toward something a reactor shield calculation would recognize.