Simulation of point kinetics reactor equations

Hello!Those two equations are needed to be solved (the attached picture)

The initial conditions n(0)=0.1, c(0)=0

The required: find the required time to increase n from 0.1 to 1

I got errors regarding syms functions. I am not sure if I do that right.I attached to this a matlab file which contains all the parameters and what I tried to do.

Answer 1

[Alan Stevens]

There are seven equations if you are using all six delayed neutron groups. You don't give your reactivity, nor the individual beta values. The program below uses arbitrary data for rho and the timespan, and Glasstone and Sesonske values for beta_i. If you are only interested in the case of a single group of delayed neutrons you should be able to modify the following appropriately:

% Point reactor kinetics 6 groups of delayed neutrons
beta = [0.00021; 0.00141; 0.00127; 0.00255; 0.00074; 0.00027]; % Taken from Glasstone and Sesonske
betasum = sum(beta);

rho = 1.1*betasum; % reactivity

% Initial consitions
c0 = zeros(6,1);
n0 = 0.1;

tspan = [0, 1];

nc0 = [n0; c0];
[t, nc] = ode45(@(t,nc) kinetics(t,nc,rho,beta,betasum), tspan, nc0);
n = nc(:,1);
c = nc(:, 2:7);

plot(t,n),grid,
xlabel('time'), ylabel('n')

% figure
% plot(t,c),grid
% xlabel('time'), ylabel('c')

function dncdt = kinetics(~,nc,rho,beta,betasum)
        
        L = 0.0001;
        lam = [0.0126; 0.0337; 0.111; 0.301; 1.14; 3.01];
        n = nc(1);
        c = nc(2:7);
        dndt = (rho - betasum)/L + sum(lam.*c);
        dcdt = beta*n/L - lam.*c;
        
        dncdt = [dndt; dcdt];        

end

[Hadeer Abdullah]

Thank you for your answer. I need to simulate all the six groups. The reactivity is 0.002

How would that change if I need to make the simulation for all six groups?

[Alan Stevens]

The program already deals with all 6 groups. Just change rho to 0.002.

[Hadeer Abdullah]

Thank you! that's very helpful.

Another point, If i want to change the reactivity from constant 0.002 to be a 0.002 increase over 1 second. I tried that:

function [rho] = InReactivity(t)
if 0<t<1
    rho = (0.002/betasum)+t
else
    rho = 0.002/betasum
end
end

But it gives me error with the ode45

beta = [0.00021; 0.00141; 0.00127; 0.00255; 0.00074; 0.00027];
>> betasum = sum(beta);
>> [rho] = (@(rho) InReactivity(t));
>> c0 = zeros(6,1);
>> tspan = [0, 0.1];
>> nc0 = [n0; c0];
>> [t, nc] = ode45(@(t,nc) kinetics(t,nc,rho,beta,betasum), tspan, nc0);

[Alan Stevens]

The following works

% Point reactor kinetics 6 groups of delayed neutrons
beta = [0.00021; 0.00141; 0.00127; 0.00255; 0.00074; 0.00027]; % Taken from Glasstone and Sesonske
betasum = sum(beta);

% Initial consitions
c0 = zeros(6,1);
n0 = 0.1;

tspan = [0, 1];

nc0 = [n0; c0];
[t, nc] = ode45(@(t,nc) kinetics(t,nc,beta,betasum), tspan, nc0);
n = nc(:,1);
c = nc(:, 2:7);

plot(t,n),grid,
xlabel('time'), ylabel('n')
% figure
% plot(t,c),grid
% xlabel('time'), ylabel('c')

function dncdt = kinetics(t,nc,beta,betasum)
        rho = InReactivity(t,betasum);
        L = 0.0001;
        lam = [0.0126; 0.0337; 0.111; 0.301; 1.14; 3.01];
        n = nc(1);
        c = nc(2:7);
        dndt = (rho - betasum)/L + sum(lam.*c);
        dcdt = beta*n/L - lam.*c;
        
        dncdt = [dndt; dcdt];        

end

function [rho] = InReactivity(t,betasum)
    if (0<t)&&(t<1)
        rho = (0.002/betasum)+t;
    else
        rho = 0.002/betasum;
    end
end

[Cody James]

Is it possible to simultaniously solve for multiple values of rho and plot them? For example: rho = [0.01, 0.02, 0.03]

[Alan Stevens]

Something like this?

% Point reactor kinetics 6 groups of delayed neutrons
beta = [0.00021; 0.00141; 0.00127; 0.00255; 0.00074; 0.00027]; % Taken from Glasstone and Sesonske
betasum = sum(beta);
rho = [0.01 0.02 0.03];
% Initial consitions
c0 = zeros(6,1);
n0 = 0.1;

tspan = [0, 1];

nc0 = [n0; c0];
for i = 1:numel(rho)
    [t, nc] = ode45(@(t,nc) kinetics(t,nc,beta,betasum,rho(i)), tspan, nc0);
    n = nc(:,1);
    c = nc(:, 2:7);
    figure
    plot(t,n),grid,
    xlabel('time'), ylabel('n')
    legend(['rho = ' num2str(rho(i))])
end

function dncdt = kinetics(~,nc,beta,betasum,rho)
        
        L = 0.0001;
        lam = [0.0126; 0.0337; 0.111; 0.301; 1.14; 3.01];
        n = nc(1);
        c = nc(2:7);
        dndt = (rho - betasum)/L + sum(lam.*c);
        dcdt = beta*n/L - lam.*c;
        
        dncdt = [dndt; dcdt];        

end

Answer 2

[Swu]

in "function dncdt = kinetics(~,nc,rho,beta,betasum)"

"dndt = (rho - betasum)/L + sum(lam.*c);"

should be ""dndt = (rho - betasum)* n/L + sum(lam.*c);

[Yuliang Fang]

You are right！