Serpent run yielded k$_{eff}$ of .79216 with
Dimensions based on buckling size given on pg. 39, MSRE Design and Operations, part iii, nuclear analysis
In [ ]:
k_nom = 1.0545
k_f_1144 = 1.04149
fuel_reactivity = (k_f_1144 - k_nom) / k_nom / 400
print(fuel_reactivity)
In [ ]:
k_f_g_1144 = 1.02315
total_reactivity = (k_f_g_1144 - k_nom) / k_nom / 400
print(total_reactivity)
graph_reactivity = total_reactivity - fuel_reactivity
print(graph_reactivity)
In [ ]:
from numpy import abs
print(abs(total_reactivity + 6.96e-5) / 6.96e-5 * 100)
print(abs(fuel_reactivity + 3.28e-5) / 3.28e-5 * 100)
print(abs(graph_reactivity + 3.68e-5) / 3.68e-5 * 100)
In [ ]:
from numpy import exp
alpha_fuel_faeh = 1.18e-4
alpha_fuel_kel = 1.8 * alpha_fuel_faeh
alpha_graph_faeh = 1.0e-5
alpha_graph_kel = 1.8 * alpha_graph_faeh
rho0_fuel = 2.146
rho0_graph = 1.86
T0 = 922
rho1144_fuel = rho0_fuel * exp(-alpha_fuel_kel * (1144 - T0))
rho1144_graph = rho0_graph * exp(-alpha_graph_kel * (1144 - T0))
print(rho1144_fuel)
print(rho1144_graph)
In [1]:
print(100 * .00272 / 1.02316)
My two benchmarks in trying to understand serpent buckling are msre_homogeneous_critical_question_mark and msre_homog_rad_57_b1. These should both have the exact same material definitions and reactor height and differ only in the radius of the reactor. This is in fact the case, with the former having a radius of 73.86 cm and the latter having a radius of 57.1544 cm. What did this mean for respective $k_{inf}$? Intuitively we would hope for the change to be relatively small, since $\sigma(\vec{r},E)$ is unchanged between the two simulations. However, we could get changes because of the change in the reactor flux. So:
msre_homogeneous_critical_question_mark: abs_gc_kinf = 1.54774msre_homog_rad_57_b1: abs_gc_kinf = 1.51718I would contend that this is a small enough change to be satisfactory. Moreover, as calculated in the sage_scratch notebook, the critical material buckling changes by even less:
msre_homogeneous_critical_question_mark: 1.32e-3msre_homog_rad_57_b1: 1.31e-3Below generation and fission rates taken from msre_homogeneous_critical_question_mark_res.m:
In [2]:
tot_gen_rate = 1.02331
tot_fiss_rate = 4.2e-1
tot_gen_rate / tot_fiss_rate
Out[2]:
Above value perfectly matches one-group nubar: 2.43645. It's obvious that in Serpent, the total loss rate is normalized to one, and consequently the generation rate is exactly equal to $k_{eff}$.
Relative loss mechanisms differ between the two different sizes as one would hope. Leakage more important in smaller reactor (oh really?):
57 cm:
TOT_ABSRATE (idx, [1: 2]) = [ 5.40995E-01 0.00021 ];
TOT_LEAKRATE (idx, [1: 2]) = [ 4.59005E-01 0.00025 ];74 cm:
TOT_ABSRATE (idx, [1: 2]) = [ 6.61174E-01 0.00017 ];
TOT_LEAKRATE (idx, [1: 2]) = [ 3.38826E-01 0.00033 ];I don't even think that the idea of a critical buckling "search" makes any bloody sense. That implies to me that you are going to change the materials! To me the search should consist of setting $k_{eff}$ equal to one, and then calculating the appropriate geometric buckling. That's what makes sense to me. I just don't think that the buckling search works in Serpent. I decreased the geometric buckling, attempting to reach a value of 1.31e-3, starting from 3.52e-3 for the 57 cm radius reactor but didn't get there...I only got down to 1.90e-3 with the reactor radius of 74 cm. However, that was plenty far enough to go below the material buckling resulting in a super-critical reactor. That's why I don't think the buckling method in Serpent works.
Ok, we're at assembly level (or some other arbitrary level). Typically, even at assembly level we use reflective boundary conditions at first to isolate the assembly from the rest of the reactor core. Treating the current level this way yields a $k_{inf}$ of sorts and a spectrum that we can call an infinite medium spectrum because of the reflective conditions. However, we want to determine a criticality spectrum. A critical system can be realized in a few ways. One way would be to add some amount of absorber (if $k_{inf}>1$; negative absorber if the system is sub-critical). However, a more physically realistic way is to add leakage; e.g. if the reactor core is critical but one piece's $k_{inf}$ is subcritical, then there will be neat leakage into that component to balance.
Two group discretization results in two eigenvalues and their associated eigenvectors (e.g. two eigenvectors total) as long as the matrix has full rank.
It might be critical to understand that in general reactor theory language when someone talks about the neutron spectrum, s/he is talking about the neutron energy distribution and may very well not be referring at all to the spatial distribution of neutrons. Spectrum -> energy. And so an infinite spectrum may refer only to the neutron energy distribution produced by an infinite medium, and consequently the critical spectrum may only refer to the change in the neutron energy distribution from the infinite spectrum. E.g. we construct a relation between the leakage and infinite spectra (Stammler pg. 356):
\begin{equation} r_L = r_{\infty} \frac{1+x}{1 + xk_{\infty}} \end{equation}where $r = \frac{\phi_1}{\phi_2}$, $r_{\infty} = \frac{\Sigma_{12}}{\Sigma_{a2}}$ and $x = L^2/\tau$ where:
\begin{equation} L^2 = \frac{D_2}{\Sigma_{a2}}\\ \tau = \frac{D_1}{\Sigma_{12} + \Sigma_{a1}} \end{equation}Note that the above equation for $r_L$ is valid only when the overall reactor is critical and only thermal neutrons cause fission events. A more general equation for the spectrum index is:
\begin{equation} r_i = \frac{r_{\infty}}{1+L^2 B_i^2} \end{equation}where i equals either 1 or 2 for the two group case and the equations for $B_1^2$ and $B_2^2$ can be found on page 355 of Stammler. Spectrum eigenvalue problem corresponds to material buckling (largest of which is $B_1^2$ and is called the fundamental or asymptotic mode; for two group problems $B_2^2$ is referred to as the transient mode; presumably for a G group problem, anything other than $B_1^2$ would be referred to as a transient mode.) Spatial eigenvalue problem corresponds to geometric buckling. Here's $B_1^2$:
\begin{equation} B_1^2 = \frac{\lambda k_{\infty} - 1}{\tau + \frac{\tau}{\tau*}L^2} \end{equation}With only thermal fission, $\tau* = \tau$. Note that $B_1^2$ is a function only of the material data and the multiplication factor of the reactor. In general, the corresponding Helmholtz equation for the spatial shape can be solved with appropriate boundary conditions determined by continuity of flux and current (assuming we are currently only examining one part of the reactor). However, if we require that the flux in all groups must vanish at the boundaries of a bare homogeneous system, then the Helmholtz equation for the spatial shape becomes an eigenvalue problem.
It should be noted from a math standpoint that eigenvalue/eigenfunction problems can only arise in differential equations when boundary conditions are homogeneous.
What I want to figure out tomorrow is exactly how the calculation of few group cross sections is implemented in serpent; e.g. look at the actual code calculation routines.
Homogenized cross section of reaction $i$:
\begin{equation} i = j + 1 \end{equation}\begin{equation} \Sigma_{i,g} = \frac{\int_V \int_{E_g}^{E_{g-1}} \Sigma_i(\vec{r}, E)\ \phi(\vec{r},E)\ d^3r\ dE}{\int_V \int_{E_g}^{E_{g-1}} \phi(\vec{r},E)\ d^3r\ dE} \quad eq. 17.3.1.1 \end{equation}Deterministic transport version:
\begin{equation} \Sigma_{i,g} = \frac{ \sum_j \left( V_j \sum_h \left( \Sigma_{i,j,h} \Phi_{j,h} \right) \right) } { \sum_j \left( V_j \sum_h \Phi_{j,h} \right) } \quad eq. 17.3.1.2 \end{equation}In a cell with reflective boundary conditions and homogeneous material composition, the flux has no spatial dependence and is thus constant with respect to the spatial coordinates. However, it will not be constant with respect to a continuous energy variable. In general the determinstic transport calculation proceeds like this, similar to that outlined in Figure I.3 in Stammler (pg. 42):
Or put another way:
Another very useful diagram is Stammler X.4. It's perhaps the best one I've seen
Understanding Serpent output:
Calculation of RES_MAT_BUCKL (where is the index looping over the few groups):
if ((div = BufVal(RES_FG_L2, i)) > 0.0)
{
val = (BufVal(RES_ABS_GC_KINF, i) - 1.0)/div;
AddStat(RES_MAT_BUCKL, i, val);
Note that just because BufVal takes an index doesn't mean that the first argument is a true array, e.g. RES_ABS_GC_KINF only has one value (e.g. one mean and one measure of variance).
Ok, here is how RES_ABS_GC_KINF is calculated:
if ((loss = abs - n2n) > 0.0)
{
keff = src/loss;
AddStat(RES_ABS_GC_KINF, i, keff);
}
In [5]:
b1_remxs = 2.37419e-3
b1_nsf = 3.66701e-3
b1_diff = 9.85681e-1
b1_mat_buckl = (b1_nsf - b1_remxs) / b1_diff
b1_k_inf = b1_nsf / b1_remxs
print(b1_mat_buckl)
print(b1_k_inf)
Above matches perfectly with Serpent output B1_BUCKLING. So how does this compare with the original calculation where the material buckling was actually higher and yet the reactor was still sub-critical??
In [6]:
remxs = 2.21079e-3
nsf = 3.35432e-3
diff = 5.31584e-1
mat_buckl = (nsf - remxs) / diff
k_inf = nsf / remxs
print(mat_buckl)
print(k_inf)
So to move the reactor to a state of criticality (from a sub-critical state), the medium was made less absorptive, less fissile, and more diffusive. More particularly, $k_{\infty}$ was increased and the medium was made more diffusive. Ok, that seems like it could work.
See sage_scratch right before the beginning of 2/24/17 notes for some notes on buckling results. I believe that the quite large discrepancy between predicted critical buckling in Serpent and the actual buckling values is entirely due to large discrepancies in diffusion coefficient values. Let's consider the results from msre_homog_rad_57_b1_res.m. We have:
If we create a cylindrical reactor that from diffusion theory should have a geometric buckling equal to the predicted critical buckling predicted in the 57 cm radius Serpent simulation, then we produce this set of diffusion coefficients with Serpent:
Here's also the important b1 results for the 73 cm radius (supposed to be critical but isn't) simulation:
In [8]:
b1_73_remxs = 2.38026e-3
b1_73_nsf = 3.67837e-3
b1_73_diff = 9.85825e-1
mat_buckl_73 = (b1_73_nsf - b1_73_remxs) / b1_73_diff
k_inf_73 = b1_73_nsf / b1_73_remxs
print(mat_buckl_73)
print(k_inf_73)
All the above results are very close to the b1 calculation results for the 57 cm reactor. So my question: how is the b1 diffusion coefficient calculated and how does it differ from calculation of DIFFCOEFF?
Serpent 1 and 2 results for the same input file are pretty much identical, which I suppose should be expected. However, the Serpent 2 result file does contain a lot more information.
Analog estimator: score the physical interactions (fission, capture, etc.) or event sequences that are simulated during the calculation. Analog estimates are directly related to the simulation process.
Note: in general when the chain of equalities ends, that generally means I've reached a point where all arguments are generated via Monte Carlo analog or implicit estimates, e.g. none of the arguments are derived quantities, e.g. if I search in collectresults.c, their first hit will be on the RHS of an equality sign.
Default calculation of $D$:
\begin{equation} D_g = \frac{\bar{r_g^2}}{6}\ \Sigma_{r,g} \end{equation}where $\bar{r_g^2}$ is the mean value of all scored square distances that neutrons have made within energy group $g$, which is an analog estimator. The starting points are the locations where neutrons have entered the group by fission or scattering from another (higher) energy group. The end points are the locations at which the neutrons are lost from the group either by absorption or scattering to another (lower) energy group.
Leakage calculation of $D$:
\begin{equation} D_g = \frac{leak}{B_m^2}\\ B_m^2 = \frac{k_{\infty} - 1}{L_g^2}\\ L_g^2 = \frac{\bar{r_g^2}}{6}\\ k_{\infty} = \frac{src}{loss} = \frac{\Sigma_{g,f}\bar{\nu}}{\Sigma_{g,f} + \Sigma_{g,c} - \Sigma_{n\rightarrow n}} \end{equation}where $leak$ comes directly from Monte Carlo scoring of neutrons exiting the simulation domain I believe.
P1 calculation of $D$:
\begin{equation} D_g = \frac{1}{3\Sigma_{g,tr}}\\ \Sigma_{g,tr} = \Sigma_{g,t} - \Sigma_{g,s1} \end{equation}B1 calculation of $D$ (note that this occurs in b1solver.c which is different from all the other calculations of $D$ that occur in collectresults.c; moreover the below equation comes from the Fridman and Lepannen PHYSOR paper as opposed to directly from the source code):
where $J_g$ is the group current and $|B|$ is the magnitude of the calculated critical buckling.
Ok, B1 calculations just don't make sense in multiple regions in Serpent. Serpent performs criticality searches in both regions independently. So that's out!
In [2]:
40 * 20000
Out[2]:
In [3]:
(381.99 + 160.60) / 4
Out[3]:
In [4]:
381.99 + 160.60
Out[4]: