Simulation of point kinetics reactor equations
41 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Hadeer Abdullah
el 7 de Oct. de 2021
Respondida: Swu
el 7 de Mzo. de 2023
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.
0 comentarios
Respuesta aceptada
Alan Stevens
el 7 de Oct. de 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
6 comentarios
Cody James
el 31 de Oct. de 2021
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
el 1 de Nov. de 2021
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
Más respuestas (1)
Swu
el 7 de Mzo. de 2023
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);
0 comentarios
Ver también
Categorías
Más información sobre Assembly en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!