Two-Group Neutronics (Neutronics)¶

This benchmark uses automatic adaptivity to solve a system of weakly coupled elliptic PDEs describing diffusion of neutrons through given medium. It employs the simple (yet often used in practice) two-group approximation by which all neutrons are divided into two distinct groups according to their energy (speed). This leads to the system of two equations shown below.

Model problem¶

Equations solved:

(1) $- \nabla \cdot D_1 \nabla \phi_1 + \Sigma_{r1}\phi_1 - \nu\Sigma_{f1} \phi_1 - \nu\Sigma_{f2} \phi_2 - Q_1 = 0,\\ - \nabla \cdot D_2 \nabla \phi_2 + \Sigma_{r2}\phi_2 - \Sigma_{s,2\leftarrow 1} \phi_1 - Q_2 = 0.$

For the physical meaning of the various material parameters, see the example 4-Group Neutronics. Their numerical values for this benchmark are given below.

Domain of interest:

Material properties¶

Piecewise constant material properties for the four regions of the domain (reactor core) and each energy group are specified by the following code:

const double D[4][2]  = { {1.12, 0.6},
                          {1.2, 0.5},
                          {1.35, 0.8},
                          {1.3, 0.9}  };
const double Sr[4][2] = { {0.011, 0.13},
                          {0.09, 0.15},
                          {0.035, 0.25},
                          {0.04, 0.35}        };
const double nSf[4][2]= { {0.0025, 0.15},
                          {0.0, 0.0},
                          {0.0011, 0.1},
                          {0.004, 0.25}       };
const double chi[4][2]= { {1, 0},
                          {1, 0},
                          {1, 0},
                          {1, 0} };
const double Ss[4][2][2] = {
                             { { 0.0, 0.0 },
                               { 0.05, 0.0 }  },
                             { { 0.0, 0.0 },
                               { 0.08, 0.0 }  },
                             { { 0.0, 0.0 },
                               { 0.025, 0.0 } },
                             { { 0.0, 0.0 },
                               { 0.014, 0.0 } }
                           };

Boundary conditions¶

Typical conditions for nuclear reactor core calculations are used:

zero Neumann on left and top edge (axes of symmetry),
zero Dirichlet on bottom edge (neutron-inert medium around the reactor core),
Newton condition on right edge (neutron-reflecting medium around the core):

$-\frac{\partial D_1\phi_1}{\partial n} = \gamma_1 \phi_1, \quad\quad -\frac{\partial D_2\phi_2}{\partial n} = \gamma_2 \phi_2,$

where the reflector albedo $\gamma$ is given by the exact solution and is equal for both groups to 8.

Exact solution¶

Quite complicated, see the source code.

Right-hand side¶

Obtained by inserting the exact solution into the equation. The function get_material is used to obtain the material marker given the physical coordinates (see main.cpp).

The following picture shows the two right-hand side functions (distribution of neutron sources/sinks) - $Q_1$ is plotted on the left, $Q_2$ on the right.

Weak formulation¶

Weak formulation of the present two-group neutron diffusion problem with fixed source terms may be derived from the general multigroup formulation shown in the 4-Group Neutronics example. Concerning its implementation (see the file forms.cpp), it is worth noticing that we manually define a higher integration order for the volumetric linear forms to correctly integrate the non-polynomial source terms, although we may set it lower for the group-1 equations than for the group-2 equations as $Q_1$ is much smoother than $Q_2$ :

Sample results¶

The following figures show the computed distributions of neutron flux for both neutron groups.

Notice the largely different behavior of the two solution components, where the first one is quite smooth while the other one more oscillating. It reflects the typical behavior observed in real cases, which arises from the different rate of interactions of fast (1^st group) and slow (2^nd group) neutrons with surrounding nuclei. This makes multimesh a preferred choice for automatic adaptivity, as can be clearly seen from the first of the series of convergence comparisons presented below.

In each convergence comparison, the reported error is the true approximation error calculated wrt. the exact solution given above and measured in a H¹ norm. The calculation was ended when the energy error estimate (often used to guide adaptivity in real multiphysics problems where exact solution is not known) became lower than 0.1%.