76 views (last 30 days)

Show older comments

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.

Alan Stevens
on 7 Oct 2021

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

Alan Stevens
on 1 Nov 2021 at 8:05

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

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!