% This LaTeX was auto-generated from MATLAB code.
% To make changes, update the MATLAB code and export to LaTeX again.

\documentclass{article}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{graphicx}
\usepackage{color}
\usepackage{hyperref}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{epstopdf}
\usepackage[table]{xcolor}
\usepackage{matlab}
\usepackage[paperheight=795pt,paperwidth=614pt,top=72pt,bottom=72pt,right=72pt,left=72pt,heightrounded]{geometry}

\sloppy
\epstopdfsetup{outdir=./}
\graphicspath{ {./exam1_media/} }

\begin{document}

\matlabtitle{EXAM 1 - NUCE 2101}

\begin{matlabcode}
clear all
\end{matlabcode}

\begin{par}
\begin{flushleft}
Dane Sabo
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
October 5, 2025
\end{flushleft}
\end{par}

\matlabheading{Problem 1}

\matlabheadingtwo{Part A}

\begin{par}
\begin{flushleft}
Usually, the in-hour equation looks like:
\end{flushleft}
\end{par}

\begin{par}
$$\rho =\omega \Lambda +\sum_{i=1}^6 \frac{\omega \beta_i }{\omega +\lambda_i }$$
\end{par}

\begin{par}
\begin{flushleft}
but for us with 2 groups and $\Lambda =50$microseconds, we have:
\end{flushleft}
\end{par}

\begin{matlabcode}
syms omega rho
sympref('FloatingPointOutput',1)
\end{matlabcode}
\begin{matlaboutput}
ans = logical
   1

\end{matlaboutput}
\begin{matlabcode}
Lambda = 5e-5; %s
beta = [0.00168272, 0.00476728];
lambda = [0.0254904, 0.268012]; %Hz


in_hour = rho == omega*Lambda + omega*beta(1)/(omega + lambda(1))...
    + omega*beta(2)/(omega + lambda(2))
\end{matlabcode}
\begin{matlabsymbolicoutput}
in\_hour = 

\hskip1em $\displaystyle \rho =\text{5.0000e-05}\,\omega +\frac{0.0048\,\omega }{\omega +0.2680}+\frac{0.0017\,\omega }{\omega +0.0255}$
\end{matlabsymbolicoutput}

\begin{par}
\begin{flushleft}
Now we can solve for the largest roots for different values of $\rho$:
\end{flushleft}
\end{par}

\begin{matlabcode}
rho_beta = [-10, -0.25, -0.1, 0, 0.1, 0.25]*(sum(beta))
\end{matlabcode}
\begin{matlaboutput}
rho_beta = 1x6
   -0.0645   -0.0016   -0.0006         0    0.0006    0.0016

\end{matlaboutput}
\begin{matlabcode}
for i = rho_beta
    in_hour_eq = subs(in_hour, rho, i); % Substitute rho with current value of i
    soln = solve(in_hour_eq, omega);
    largest_root = max(double(soln)); % Convert symbolic solution to double and find the largest root
    fprintf('Largest root for rho = %.4f is: %.5f\n', i/sum(beta), largest_root);
end
\end{matlabcode}
\begin{matlaboutput}
Largest root for rho = -10.0000 is: -0.02484
Largest root for rho = -0.2500 is: -0.01156
Largest root for rho = -0.1000 is: -0.00613
Largest root for rho = 0.0000 is: 0.00000
Largest root for rho = 0.1000 is: 0.00998
Largest root for rho = 0.2500 is: 0.03847
\end{matlaboutput}

\begin{par}
\begin{flushleft}
Compared to the class solution for the six group in-hour equation on page 15 of Fundamental Kinetics Ideas V17 (FKI), we achieve an interesting result.
\end{flushleft}
\end{par}

\begin{enumerate}
\setlength{\itemsep}{-1ex}
   \item{\begin{flushleft} For $\rho =-10$, two groups significantly underpredicts the largest root (roughly twice the magnitude!). \end{flushleft}}
   \item{\begin{flushleft} For $\rho =-0.25$, two groups underpredicts the largest root. \end{flushleft}}
   \item{\begin{flushleft} For $\rho =-0.1$, two groups underpredicts the largest root. \end{flushleft}}
   \item{\begin{flushleft} For $\rho =0$, two groups matches the six-group solution exactly (both predict a marginally stable system with $\omega = 0$). \end{flushleft}}
   \item{\begin{flushleft} For $\rho =0.1$, two groups is close, but still underpredicts the largest root. \end{flushleft}}
   \item{\begin{flushleft} For $\rho =0.25$, two groups underpredicts the largest root. \end{flushleft}}
\end{enumerate}

\begin{par}
\begin{flushleft}
It appears that our solution always underpredicts the six-group equation, and suffers more error further away from $\rho =0$.
\end{flushleft}
\end{par}

\matlabheadingtwo{Part B}

\begin{par}
\begin{flushleft}
\boxed{A=\left\lbrack \begin{array}{ccc}
\frac{\rho -(\beta_1 +\beta_2 )}{\Lambda } & \lambda_1  & \lambda_2 \\
\beta_1 /\Lambda  & -\lambda_1  & 0\\
\beta_2 /\Lambda  & 0 & -\lambda_2 
\end{array}\right\rbrack}
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
\boxed{X_0 =\frac{n_0 }{\Lambda }\left\lbrack \begin{array}{c}
\Lambda \\
\beta_1 /\lambda_1 \\
\beta_2 /\lambda_2 
\end{array}\right\rbrack}
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
\boxed{S=\left\lbrack \begin{array}{c}
S_0 \\
0\\
0
\end{array}\right\rbrack}
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
with $X=\left\lbrack \begin{array}{c}
n(t)\\
C_1 (t)\\
C_2 (t)
\end{array}\right\rbrack$
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
These results are easily obtained by modifying the matrices on page 17 of the Fundamental Kinetics Ideas V17 handout. For our case, we can trim the matrices down assuming that $\beta_3 ....\beta_6$ and $\lambda_3 ...\lambda_6$ are all zero, and then remove the padded rows and columns of zeros to obtain our result.
\end{flushleft}
\end{par}

\matlabheadingtwo{Part C}

\begin{par}
\begin{flushleft}
Our system in part B can be used to create the system of differential equations:
\end{flushleft}
\end{par}

\begin{par}
$$\frac{dX}{dt}=AX+S$$
\end{par}

\begin{par}
\begin{flushleft}
This is a simple first order system, and as such has standard solution:
\end{flushleft}
\end{par}

\begin{par}
$$\boxed{X(t)=e^{At} X_0 +e^{At} \int_0^t e^{-A\tau } S(t)d\tau}$$
\end{par}

\begin{par}
\begin{flushleft}
as found on page 17 of FKI.
\end{flushleft}
\end{par}

\matlabheadingtwo{Part D}

\begin{matlabcode}
A = [(0.25*sum(beta) - sum(beta))/Lambda, lambda(1), lambda(2);
    beta(1)/Lambda,                       -lambda(1), 0;
    beta(2)/Lambda,                       0,          -lambda(2)]
\end{matlabcode}
\begin{matlaboutput}
A = 3x3
  -96.7500    0.0255    0.2680
   33.6544   -0.0255         0
   95.3456         0   -0.2680

\end{matlaboutput}
\begin{matlabcode}
%The modal matrix is P, while the eigenvalues are the diagonal entries
%of D.
[P, D] = eig(A)
\end{matlabcode}
\begin{matlaboutput}
P = 3x3
   -0.6915    0.0016    0.0009
    0.2399    0.8608   -0.9104
    0.6814    0.5090    0.4138

D = 3x3
  -97.0230         0         0
         0    0.0385         0
         0         0   -0.0590

\end{matlaboutput}
\begin{matlabcode}
% Calculate the matrix exponential of A multiplied by tau
%Here is e^A*tau
tau = 1e-4; %s
part_d_matrix = P*expm(D*tau)*inv(P)
\end{matlabcode}
\begin{matlaboutput}
part_d_matrix = 3x3
    0.9904    0.0000    0.0000
    0.0033    1.0000    0.0000
    0.0095    0.0000    1.0000

\end{matlaboutput}
\begin{matlabcode}
%expm is used because techinically D is a matrix here.
%exp(trace(D) * tau) is equivalent
proof = P*exp(trace(D) * tau)*inv(P)
\end{matlabcode}
\begin{matlaboutput}
proof = 3x3
    0.9903   -0.0000   -0.0000
         0    0.9903    0.0000
    0.0000         0    0.9903

\end{matlaboutput}

\matlabheadingtwo{Part E}

\begin{matlabcode}
%Initial conditions
n_0 = 1; %this is an assumption
X_0 = n_0/Lambda*[Lambda; beta(1)/lambda(1); beta(2)/lambda(2)];
S = [0;0;0];

% Define the time span for the ODE solution
tspan = [0, 30]; % seconds

% Define the ODE function
part_e = @(t, X) A * X + S;

% Solve the ODE using ode45
[t_e, X] = ode45(part_e, tspan, X_0);

% Plot the first state of X (neutron population) over time
figure;
semilogy(t_e, X(:, 1)/n_0);
xlabel('Time (s)');
ylabel('Neutron Population (log scale)');
title('Neutron Population over Time');
subtitle('(Two Precursors)');
grid on;
\end{matlabcode}
\begin{center}
\includegraphics[width=\maxwidth{56.196688409433015em}]{figure_0.png}
\end{center}

\matlabheadingtwo{Part F}

\begin{matlabcode}
lambda_eff = 0.1 %Hz
\end{matlabcode}
\begin{matlaboutput}
lambda_eff = 0.1000
\end{matlaboutput}
\begin{matlabcode}
rho = 0.25*sum(beta)
\end{matlabcode}
\begin{matlaboutput}
rho = 0.0016
\end{matlaboutput}

\begin{par}
\begin{flushleft}
The prompt jump makes an assumption about the value of the precursor, and assumes that prompt neutrons immediately equalize. Using the formula given in the exam:
\end{flushleft}
\end{par}

\begin{matlabcode}
prompt_jump = n_0 * sum(beta)/(sum(beta)-rho)
\end{matlabcode}
\begin{matlaboutput}
prompt_jump = 1.3333
\end{matlaboutput}

\begin{par}
\begin{flushleft}
And then we use this value as our our initial condition for the neutron population.
\end{flushleft}
\end{par}

\begin{matlabcode}
% Define the ODE function
part_f = @(t, X) lambda_eff*rho/(sum(beta) - rho);

% Solve the ODE using ode45
[t_f, n] = ode45(part_f, tspan, prompt_jump);

% Plot the first state of X (neutron population) over time
figure;
semilogy(t_f, n/n_0);
xlabel('Time (s)');
ylabel('Neutron Population (log scale)');
title('Neutron Population over Time ');
subtitle('(Single Precursor, Prompt Jump)');
grid on;
\end{matlabcode}
\begin{center}
\includegraphics[width=\maxwidth{56.196688409433015em}]{figure_1.png}
\end{center}

\matlabheadingtwo{Part G}

\begin{par}
\begin{flushleft}
I've rewritten SimplePower.m using ode45 for clarity. This is the six group equivalent solution, compared to the other two solutions I made previously.
\end{flushleft}
\end{par}

\begin{matlabcode}
% Physical parameters for 6-group model
beta_6 = [0.00021, 0.00141, 0.00127, 0.00255, 0.00074, 0.00027];  
% delayed neutron fractions
lambda_6 = [0.01246403, 0.03052863, 0.11141479, 0.30130435, 1.13606557, 3.01304348];  
% decay constants [1/sec]

% Build system matrix A_6
A_6 = zeros(7, 7);
A_6(1, 1) = (rho - sum(beta_6)) / Lambda;
A_6(1, 2:7) = lambda_6;
for i = 1:6
    A_6(i+1, 1) = beta_6(i) / Lambda;
    A_6(i+1, i+1) = -lambda_6(i);
end

% Initial conditions (steady state)
X_0_6 = n_0 / Lambda * [Lambda; beta_6(1)/lambda_6(1); 
beta_6(2)/lambda_6(2);
beta_6(3)/lambda_6(3);
beta_6(4)/lambda_6(4);
beta_6(5)/lambda_6(5);
beta_6(6)/lambda_6(6)];

% Source term (zero for this problem)
S_6 = zeros(7, 1);

% Define the ODE function
part_g = @(t, X) A_6 * X + S_6;

% Solve the ODE using ode45
[t_6, X_6] = ode45(part_g, tspan, X_0_6);
% Combined plot of all three models
figure;
hold on;

% Two-group model (from Part E)
semilogy(t_e, X(:, 1)/n_0, 'b-', 'LineWidth', 1.5, 'DisplayName', 'Two Groups');

% Single group prompt jump (from Part F)
semilogy(t_f, n/n_0, 'r--', 'LineWidth', 1.5, 'DisplayName',...
'Single Group (Prompt Jump)');

% Six-group model (from Part G)
semilogy(t_6, X_6(:, 1)/n_0, 'g-.', 'LineWidth', 1.5, 'DisplayName', 'Six Groups');

hold off;
xlabel('Time (s)');
ylabel('Neutron Population (normalized, log scale)');
title('Comparison of Point Kinetics Models');
subtitle(sprintf('\\rho = %.3f\\beta', rho/sum(beta)));
legend('Location', 'best');
grid on;
\end{matlabcode}
\begin{center}
\includegraphics[width=\maxwidth{56.196688409433015em}]{figure_2.png}
\end{center}

\matlabheadingtwo{Part H}

\begin{par}
\begin{flushleft}
We can create a two-group prompt jump assumption model by adapting the six group prompt jump model on page 34 of FKI:
\end{flushleft}
\end{par}

\begin{par}
$$\boxed{N(t)=\frac{\Lambda \left(\sum_{i=1}^2 \lambda_i C_i +S\right)}{\beta -\rho }}$$
\end{par}

\begin{par}
\begin{flushleft}
where each precursor is a first order differential equation:
\end{flushleft}
\end{par}

\begin{par}
$$\boxed{\frac{dC_i (t)}{dt}=\frac{N(t)\beta_i }{\Lambda }-\lambda_i C_i (t)}$$
\end{par}

\begin{par}
\begin{flushleft}
We have two precursors, so we then can build a matrix form of the two equations together:
\end{flushleft}
\end{par}

\begin{par}
$$\boxed{\left\lbrack \begin{array}{c}
{\dot{C} }_1 \\
{\dot{C} }_2 
\end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
\beta_1 \lambda_1 -\lambda_1 (\beta -\rho ) & \beta_1 \lambda_2 \\
\beta_2 \lambda_1  & \beta_2 \lambda_2 -\lambda_2 (\beta -\rho )
\end{array}\right\rbrack \left\lbrack \begin{array}{c}
C_1 \\
C_2 
\end{array}\right\rbrack +\left(\frac{S}{\beta -\rho }\right)\left\lbrack \begin{array}{c}
\beta_1 \\
\beta_2 
\end{array}\right\rbrack}$$
\end{par}

\begin{par}
\begin{flushleft}
where the initial condition for each precursor is:
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
with initial condition $$\boxed{{\vec{C} }_0 =\frac{N_{S/D} }{\Lambda }\left\lbrack \begin{array}{c}
\beta_1 /\lambda_1 \\
\beta_2 /\lambda_2 
\end{array}\right\rbrack}$$ where $N_{S/D}$ is derived from the shutdown equilibrium.
\end{flushleft}
\end{par}


\matlabheading{Problem 2}

\begin{matlabcode}
clear all
\end{matlabcode}

\matlabheadingtwo{Part A}

\begin{par}
\begin{flushleft}
Let's assume for this problem, the 'sudden insertion' of reactivity is an instant addition. In that case, the prompt neutrons will quickly increase the total neutron population in the first few milliseconds. After they very quickly stabilize, the prompt neutrons will play a back-seat role to the delayed neutrons. A fraction of these neutrons and fission reactions are becoming precursors, which will release delayed neutrons on a much slower timescale. These delayed neutrons will continue to increase power, as the reactor is delayed critical with this reactivity addition. 
\end{flushleft}
\end{par}

\matlabheadingtwo{Part B}

\begin{par}
\begin{flushleft}
The equation becomes 
\end{flushleft}
\end{par}

\begin{par}
$$\boxed{\frac{N(0_+ )}{N(0_- )}=\frac{\beta -\rho (0_- )}{\beta -\rho (0_+ )}}$$
\end{par}

\begin{par}
\begin{flushleft}
Usually we just assume that $\rho (0_- )=0$, but with sources involved to maintain prompt neutron steady state, $\rho (0_- )$ is actually negative. This is supported by page 32 of FKI.
\end{flushleft}
\end{par}

\matlabheadingtwo{Part C}

\begin{par}
\begin{flushleft}
For the initial jump power, we can use the prompt jump assumption from part B:
\end{flushleft}
\end{par}

\begin{matlabcode}
P_0_left = 1.0e-8; %assume neutron count and power count are proportional
P_0_right = P_0_left *(1-(-10))/(1-(-9)) %beta factors out and cancels
\end{matlabcode}
\begin{matlaboutput}
P_0_right = 1.1000e-08
\end{matlaboutput}

\begin{par}
\begin{flushleft}
This is the initial jump power. 
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
To calculate the steady state value, we recognize that $\rho$ is still strongly subcritical. For that reason, the steady state value is eventually going to return to a shutdown equilibrium power. We can use the shutdown equilibrium equation on page 35 of FKI to find the steady state power:
\end{flushleft}
\end{par}

\begin{par}
$$P_{S/D} \propto-\frac{S\Lambda }{\rho }$$
\end{par}

\begin{par}
\begin{flushleft}
We assume that the source and prompt neutron period are constant, and thus can be solved for.
\end{flushleft}
\end{par}

\begin{par}
$$S\Lambda \propto-1.0\times 10^{-8} \times -10\beta$$
\end{par}

\begin{par}
\begin{flushleft}
and then plug this value back into the equation but with using $\rho =-9\beta$:
\end{flushleft}
\end{par}

\begin{par}
$$P_{S/D} (0_+ )=\frac{-(1.0\times 10^{-8} \times 10\beta )}{-9\beta }$$
\end{par}

\begin{par}
$$\boxed{P_{S/D} (0_+ )=1.111\times 10^{-8}}$$
\end{par}

\begin{par}
\begin{flushleft}
This is the final steady state power. Very slightly more than the prompt jump (\textasciitilde{}1\%).
\end{flushleft}
\end{par}

\matlabheadingtwo{Part D}

\begin{par}
\begin{flushleft}
The length of time for the transient to settle is dependent on the precursor pole locations. In this problem, we are not given values for $\lambda$, but I can make an educated guess. Using the table on page 15 of FKI, we can extract the slowest period for the in-hour equation where $\rho /\beta =-10$. This is $\omega_0 =-0.012423$. This corresponds to a half-life $t_{1/2} =\frac{0.693}{|-0.012423|}=55.7836[\textrm{s}]$.
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
For five half lives, this value is scaled by 5:
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
$\boxed{t_{ss} =5t_{1/2} =278.9181[\textrm{s}]}$.
\end{flushleft}
\end{par}


\matlabheading{Problem 3}

\begin{matlabcode}
clear all
\end{matlabcode}

\matlabheadingtwo{Part A}

\begin{par}
\begin{flushleft}
Subcritical multiplication is the process where neutron sources play a significant role in the fission chain. This is usually for low power or shutdown reactors, as the reactivty is so low the steady state power of the reactor is a function of the source activities.
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
Subcritical multiplication is also present when a reactor core that has been operating is shut down. Certain fission decay chains have products that create additional neutrons in different timescales. First, hyrdogren interacting with gamma rays produces neutrhons for about the first day, then heavier elements dominate source production thereafter.
\end{flushleft}
\end{par}

\matlabheadingtwo{Part B}

\begin{par}
\begin{flushleft}
Assuming the precursors are in steady-state, 10\% of the neutrons in the fission chain must come from sources to maintain steady state.
\end{flushleft}
\end{par}


\matlabheading{Problem 4}

\begin{matlabcode}
clear all
\end{matlabcode}

\matlabheadingtwo{Part A}

\begin{par}
$$SUR=26.06[\textrm{dpm}\textrm{-}\sec ]\frac{\dot{\rho} +\lambda_{eff} \rho }{(\beta -\rho )}$$
\end{par}

\begin{par}
\begin{flushleft}
First, I'm going to find the values for $SUR$ at certain points relating to the rod movement:
\end{flushleft}
\end{par}

\begin{matlabcode}
SUR = @(rho_dot, lambda_eff, rho, beta) 26.06*(rho_dot + lambda_eff*rho)/(beta-rho);

%before rod pull
rho_dot = 0;
lambda_eff = 0.1; %Hz
rho = 0;
beta = 0.00645; %assumed

SUR(rho_dot, lambda_eff, rho, beta)
\end{matlabcode}
\begin{matlaboutput}
ans = 0
\end{matlaboutput}
\begin{matlabcode}
%now rod pull starts
rho_dot = 0.00160/15;
SUR(rho_dot, lambda_eff, rho, beta)
\end{matlabcode}
\begin{matlaboutput}
ans = 0.4310
\end{matlaboutput}
\begin{matlabcode}
%end of rod pull
rho = 0.00160;
SUR(rho_dot, lambda_eff, rho, beta)
\end{matlabcode}
\begin{matlaboutput}
ans = 1.4329
\end{matlaboutput}
\begin{matlabcode}
%rods not moving, just chillin'
rho_dot = 0;
SUR(rho_dot, lambda_eff, rho, beta)
\end{matlabcode}
\begin{matlaboutput}
ans = 0.8597
\end{matlaboutput}
\begin{matlabcode}
%start putting rods back in
rho_dot = -0.00160/15;
SUR(rho_dot, lambda_eff, rho, beta)
\end{matlabcode}
\begin{matlaboutput}
ans = 0.2866
\end{matlaboutput}
\begin{matlabcode}
%end of rod push
rho = 0;
SUR(rho_dot, lambda_eff, rho, beta)
\end{matlabcode}
\begin{matlaboutput}
ans = -0.4310
\end{matlaboutput}
\begin{matlabcode}
%rods back in to zero
rho_dot = 0;
SUR(rho_dot, lambda_eff, rho, beta)
\end{matlabcode}
\begin{matlaboutput}
ans = 0
\end{matlaboutput}

\matlabheadingtwo{\includegraphics[width=\maxwidth{60.81284495735073em}]{image_0}}

\matlabheadingtwo{Part B}

\begin{par}
\begin{flushleft}
I would expect power to reach its maximum value at about 32+-1 seconds from the beginning of the operation. That is when SUR becomes negative.
\end{flushleft}
\end{par}

\matlabheadingtwo{Part C}

\begin{par}
\begin{flushleft}
We need to find rho\_dot that makes SUR 0 at the beginning of rod push.
\end{flushleft}
\end{par}

\begin{par}
$$SUR=26.06[\textrm{dpm}\textrm{-}\sec ]\frac{\dot{\rho} +\lambda_{eff} \rho }{(\beta -\rho )}$$
\end{par}


\vspace{1em}
\begin{par}
$$0=26.06[\textrm{dpm}\textrm{-}\sec ]\frac{\dot{\rho} +\lambda_{eff} \rho }{(\beta -\rho )}$$
\end{par}

\begin{par}
$$0=\dot{\rho} +\lambda_{eff} \rho$$
\end{par}

\begin{par}
$$\dot{\rho} =-\lambda_{eff} \rho$$
\end{par}

\begin{matlabcode}
rho = 0.0016;
rho_dot = -lambda_eff * rho
\end{matlabcode}
\begin{matlaboutput}
rho_dot = -1.6000e-04
\end{matlaboutput}

\begin{par}
\begin{flushleft}
We would need a rate of 1.6 pcm/s reactivity removal.
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
Check:
\end{flushleft}
\end{par}

\begin{matlabcode}
SUR(rho_dot, lambda_eff, rho, beta)
\end{matlabcode}
\begin{matlaboutput}
ans = 0
\end{matlaboutput}

\matlabheadingtwo{Part D}

\begin{par}
\begin{flushleft}
This use of the start up rate equation has significant numerical errors that limit its usefulness in reactor operations, limiting the basic SUR equation to only be used for getting a feel for transient behavior, not for exact operations. The big error is the assumption that $\lambda_{eff}$ isn't changing with time. This can be corrected by using the additional term $\frac{{\dot{\lambda} }_{eff} }{\lambda }$ as shown on page 51 of FKI, but doing this requires solving for all six precursors continuously.Using modern computational methods, this is not difficult, and achieavable with real-time simulation.
\end{flushleft}
\end{par}


\matlabheading{Problem 5}

\begin{par}
\begin{flushleft}
$\lambda_{eff}$ changes significantly with the rate of power change because the dominant precursor effect is dependent on that rate of power change and the duration of the change in power. When power is rising, the simplification of $\lambda_{eff}$ must either neglect the behavior of long time constant precursors or short time constant precursors. If the rate of power change is fast, $\lambda_{eff}$ also needs to be fast to model the dominating effect of the fast responding precursors. If the rate of power change is slow, $\lambda_{eff}$ must be slow to accurately capture the dominating behavior of the slow precursors. Generally, $\lambda_{eff}$ must keep up with the fast precursors on power increase. On power decrease, the slower precursors can drag power into the mix for longer, requiring a slower $\lambda_{eff}$.
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
This is further complicated by the duration of the power change. If the duration of a power change is long enough, a fast $\lambda_{eff}$ is not going to correctly model the contribution of large time constant precursors (and vice versa for the short duration power change case). This is true for both increases and decreases in power.
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
Basically, there are too many exponential decays and growths to be modelled accurately in all cases by one $\lambda_{eff}$. 
\end{flushleft}
\end{par}


\matlabheading{Problem 6}

\begin{par}
\begin{flushleft}
The adjoint flux $\phi^*$ represents the importance of a neutron at a given location and energy for contributing to the reactor's power generation or other quantities of interest (such as reactivity or detector response). It is highest in locations like the core center, where a neutron at the right spot is most likely to cause fission events that contribute to the global reactor behavior, and lowest near the edges of the reactor, where neutron escape becomes more likely. It's a way to account for the fact that a constant amount of flux everywhere does not create a constant power everywhere. 
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
This is super useful if there's mixed fuel rods in the reactor. The adjoint flux is lower around more spent fuel rods compared to fresh rods.
\end{flushleft}
\end{par}


\matlabheading{Problem 7}

\begin{par}
\begin{flushleft}
The equation developed in class is:
\end{flushleft}
\end{par}

\begin{enumerate}
\setlength{\itemsep}{-1ex}
   \item{\begin{flushleft} $\displaystyle \frac{d}{dt}\left\lbrack \frac{\phi (\vec{r} ,E,t)}{v}\right\rbrack =$ \end{flushleft}}
   \item{\begin{flushleft} $\displaystyle \nabla \cdot D(\vec{r} ,E)\nabla \phi (\vec{r} ,E,t)$ \end{flushleft}}
   \item{\begin{flushleft} $\displaystyle -\Sigma_t (\vec{r} ,E,t)\phi (\vec{r} ,E,t)$ \end{flushleft}}
   \item{\begin{flushleft} $\displaystyle +\int_0^{\infty } \Sigma_s (\vec{r} ,E^{\prime } \to E,t)\phi (\vec{r} ,E^{\prime } ,t)dE^{\prime }$ \end{flushleft}}
   \item{\begin{flushleft} $\displaystyle +q(\vec{r} ,E^{\prime } ,t)$ \end{flushleft}}
   \item{\begin{flushleft} $\displaystyle +\nu (1-\beta )\chi_p (E)\int_0^{\infty } \Sigma_f (\vec{r} ,E^{\prime } ,t)\phi (\vec{r} ,E^{\prime } ,t)dE^{\prime }$ \end{flushleft}}
   \item{\begin{flushleft} $\displaystyle +\sum_{i=1}^6 \chi_{di} (E)\lambda_i C_i (\vec{r} ,t)$ \end{flushleft}}
\end{enumerate}

\matlabheadingtwo{Part A}

\begin{par}
\begin{flushleft}
Here are descriptions of the physical significance of each term:
\end{flushleft}
\end{par}

\begin{enumerate}
\setlength{\itemsep}{-1ex}
   \item{\begin{flushleft} This is the time rate of change of the flux $\phi$ divided by the average neutrons speed $v$. This is analogous to the time rate of change of the total neutron population $\dot{n} (\vec{r} ,E,t)$. \end{flushleft}}
   \item{\begin{flushleft} This is the neutron leakage term. This is a diffusion process that represents the physical movement of neutrons through space. Generally, it says neutrons are going to flow in the direction of the lowest concentration of neutrons. \end{flushleft}}
   \item{\begin{flushleft}  This is the neutron removal rate term. This includes absorption and neutrons used in fission events. \end{flushleft}}
   \item{\begin{flushleft} This is the neutron in-scattering term. This describes the rate at which neutrons scatter to energy $E$ from all other energies at position $\vec{r}$ at time $t$. \end{flushleft}}
   \item{\begin{flushleft} This is the external neutron source term.  \end{flushleft}}
   \item{\begin{flushleft} This is the production rate of prompt neutrons integrated over all energies $E$, multiplied by a distribution of $\chi_p$, which is the distribution of neutron energies after a fission event. The fraction $\nu (1-\beta )$ selects only prompt neutrons. These are fuel fission neutrons. \end{flushleft}}
   \item{\begin{flushleft} This is the production rate from the delayed precursors, with six precursor groups described here, and multiplied by $\chi_{di}$ to represent the distribution of neutron energy from each precursor. \end{flushleft}}
\end{enumerate}

\matlabheadingtwo{Part B}

\begin{par}
\begin{flushleft}
First we have the equation:
\end{flushleft}
\end{par}

\begin{par}
$$\phi (\vec{r} ,E,t)=F(\vec{r} )G(t)$$
\end{par}

\begin{par}
\begin{flushleft}
with the condition that:
\end{flushleft}
\end{par}

\begin{par}
$$\nabla^2 F(\vec{r} )=-B^2 F(\vec{r} )$$
\end{par}

\begin{par}
\begin{flushleft}
We can expand this second equation by using the definition of Laplace operator:
\end{flushleft}
\end{par}

\begin{par}
$$\frac{d^2 F}{dx^2 }+\frac{d^2 F}{dy^2 }+\frac{d^2 F}{dz^2 }+B^2 F=0$$
\end{par}

\begin{par}
\begin{flushleft}
This is the equation we need to solve for $F$ in order to find the spatial part of the flux $\phi$. To do so, we need to separate variables. Recognize that $F(\vec{r} )$ can be rewritten as $X(x)Y(y)Z(z)$, where each term makes up the behavior for that one component of $\vec{r}$. Doing this, we can rewrite our diffusion equation:
\end{flushleft}
\end{par}

\begin{par}
$$YZ\frac{d^2 X}{dx^2 }+XZ\frac{d^2 Y}{dy^2 }+XY\frac{d^2 Z}{dz^2 }+B^2 XYZ=0$$
\end{par}

\begin{par}
$$\frac{1}{X}\frac{d^2 X}{dx^2 }+\frac{1}{Y}\frac{d^2 Y}{dy^2 }+\frac{1}{Z}\frac{d^2 Z}{dz^2 }+B^2 =0$$

\begin{par}
\begin{flushleft}
where $\boxed{B^2 = B_x^2 + B_y^2 + B_z^2}$ (geometric buckling).
\end{flushleft}
\end{par}
\end{par}

\begin{par}
\begin{flushleft}
which we can then break into components:
\end{flushleft}
\end{par}

\begin{par}
$$\frac{1}{X}\frac{d^2 X}{dx^2 }+B_x^2 =0,~~\frac{1}{Y}\frac{d^2 Y}{dy^2 }+B_y^2 =0,~~\frac{1}{Z}\frac{d^2 Z}{dz^2 }+B_z^2 =0,~~$$
\end{par}

\begin{par}
\begin{flushleft}
which then become
\end{flushleft}
\end{par}

\begin{par}
$$\frac{d^2 X}{dx^2 }+B_x^2 X=0,~~\frac{d^2 Y}{dy^2 }+B_y^2 Y=0,~~\frac{d^2 Z}{dz^2 }+B_z^2 Z=0,~~$$
\end{par}

\begin{par}
\begin{flushleft}
Now, we've got a bunch of 2nd order differential equations to solve, which have a standard solution. Because $B^2$ is always positive, we know the roots of the characteristic equation are always complex, and then that the solutions to these equations are going to have oscillatory solutions and the general form:
\end{flushleft}
\end{par}

\begin{par}
$$X=A\cos (B_x x)+C\sin (B_x x)$$
\end{par}

\begin{par}
\begin{flushleft}
The other spacial coordinates are exactly the same, since we have a symmetric cube.
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
Our boundary conditions are also specified. We know just outside the boundary, there is zero flux. This means:
\end{flushleft}
\end{par}

\begin{par}
$$F(\hat{L} ,\hat{L} ,\hat{L} )=\vec{0}$$
\end{par}

\begin{par}
\begin{flushleft}
We can also abuse another fact of symmetry: because our cube is symmetric, and our origin is centered in the cube, the sine part of our solution \textit{must} be zero for our solution to be symmetric. So then, we know that:
\end{flushleft}
\end{par}

\begin{par}
$$X=A\cos (B_x x)$$
\end{par}

\begin{par}
\begin{flushleft}
and with our boundary condition:
\end{flushleft}
\end{par}

\begin{par}
$$X\left(\frac{\hat{L} }{2}\right)=A\cos \left(B_x \frac{\hat{L} }{2}\right)=0$$
\end{par}

\begin{par}
\begin{flushleft}
$A$ can be any scalar.
\end{flushleft}
\end{par}

\begin{par}
$$\cos \left(B_x \frac{\hat{L} }{2}\right)=0$$
\end{par}

\begin{par}
\begin{flushleft}
$B_x \frac{\hat{L} }{2}=\frac{n\pi }{2},~~n=1,3,5,...$ We'll just assume $n=1$, and that we don't care about far outside of the cube.
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
Then, $\boxed{B_x =\frac{\pi }{\hat{L} }}$, and by symmetry, $\boxed{B_y =\frac{\pi }{\hat{L} },~~B_z =\frac{\pi }{\hat{L} }}$.
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
Thus our solutions for $X$, $Y$, and $Z$ become:
\end{flushleft}
\end{par}

\begin{par}
$$X\left(x\right)=A\cos \left(\frac{\pi x}{\hat{L} }\right),~~Y\left(y\right)=A\cos \left(\frac{\pi y}{\hat{L} }\right),~~Z\left(z\right)=A\cos \left(\frac{\pi z}{\hat{L} }\right)$$
\end{par}

\begin{par}
\begin{flushleft}
and putting things back together,
\end{flushleft}
\end{par}

\begin{par}
$$\boxed{F(x,y,z)=A^3 \cos \left(\frac{\pi x}{\hat{L} }\right)\cos \left(\frac{\pi y}{\hat{L} }\right)\cos \left(\frac{\pi z}{\hat{L} }\right)}$$
\end{par}

\end{document}