Ode15s , Plotting issue.
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
Nx = 50; % Number of grid points
L = 5000; % Length of the spatial domain
dx = L / (Nx-1); % Spatial step size
x = linspace(0, L, Nx); % Spatial points
A=0.1963;
f = 0.008;
D = 500e-3;
rho=1.2;
c = 348.5;
P_atm = 1.013e5; % Atmospheric pressure in Pa
tspan = [0, 3600]; % Time span
tcut = [tspan(1) 600 1800 tspan(end)];
mcut = [0 509.28 0];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
t_solution = [];
P_solution = [];
m_solution = [];
for i = 1:numel(mcut)
tspan = [tcut(i) tcut(i+1)] ;
if i == 1
% Initial conditions
P0 = 5e6*ones(Nx,1); % Initial condition for pressure (P)
m0 = zeros(Nx, 1); % Initial condition for mass flow rate (m)
m0(end) = mcut(1);
else
P0 = u(end, 1:Nx).';
m0 = u(end, Nx + 1:end).';
m0(end) = mcut(i);
end
% Initial conditions vector
u0 = [P0;m0];
% Solve the system of ODEs using ode15s
[t, u] = ode15s(@(t, u) MyODEs(t, u, Nx, dx, f, c, D), tspan, u0,options);
% Extract the solutions for P and m
t_solution = [t_solution;t];
P_solution = [P_solution;u(:, 1:Nx)];
m_solution = [m_solution;u(:, Nx + 1:end)];
if t(end) < tspan(2)
break
end
end
%% when m =kg/s*m^2
Q_n= (m_solution./rho);
Q=(Q_n .*A);
nodes_to_plot = 1:7:Nx;
% Plotting P_solution for selected nodes with transparency
figure;
hold on;
colors = rand(length(nodes_to_plot), 3);
for idx = 1:length(nodes_to_plot)
i = nodes_to_plot(idx);
plot(t_solution , P_solution(:, i), 'LineWidth', 0.3, 'DisplayName', ['Node ' num2str(i)], 'Color', colors(idx,:), 'LineStyle', '-');
end
hold off;
xlabel('Time (t)');
ylabel('Pressure (P)');
title('Pressure (P) at Selected Nodes over Time');
legend('show');
% % Choose Node 1 (x = 1) for plotting
node1_idx = 25;
% % Extract the pressure and Q values at Node 1
Q_node25 = Q(:, node1_idx);
figure;
plot(t_solution, Q_node25 , 'LineWidth',2);
xlabel('Time (t)');
ylabel('Q at Node 25');
title('Q at Node 25 over Time with covective term');
function dudt = MyODEs(t, u, Nx, dx, f, c, D)
P = u(1:Nx);
m = u(Nx + 1:end);
dPdt = zeros(Nx, 1);
dmdt = zeros(Nx, 1);
% Boundary Condition
P(1) = 5e6; % P(1) boundary condition
% Nozzle at i=25 (Algebraic relationships)
p_hst = 503.51e3; % p_hst = rho * g * d (where d = 40m)
Cd = 0.55;
gamma = 1.4;
rho1 = 1.2;
% Mass flowrate of nozzle
m_nz = Cd *sqrt((2 * gamma / (gamma - 1)) * P(25) * rho1 * (1 - (p_hst / P(25))^((gamma - 1) / gamma)) * (p_hst / P(25))^(2 / gamma));
%fprintf(' m_nz: %f kg/s*m^2\n', m_nz);
% Equation
dmdt(25) =-(2*m(25)*c^2/P(25)*((- 3*m(25)+4*m(26)-m(27)) / (2 * dx)))+(c^2*m(25)^2/P(25)^2)*((- 3*P(25)+ 4*P(26)-P(27)) / (2*dx)) - ((-3*P(25)+4*P(26)-P(27)) / (2*dx)) - (f * m(25) * abs(m(25)) * c^2) / (2 *D* (P(25))); % Forward discretization for dm/dt at i=25
boundary_cond_m(25)=m(25)+m_nz;
dPdt(25) =-c^2 * (3*boundary_cond_m(25) - 4*m(24) + m(23)) / (2*dx); % Backward discretization of dP/dt at i=25
% Equation
dmdt(1) =-(2*m(1)*c^2/P(1)*((- 3*m(1)+4*m(2)-m(3)) / (2 * dx)))+(c^2*m(1)^2/P(1)^2)*((- 3*P(1)+ 4*P(2)-P(3)) / (2*dx)) - ((-3*P(1)+4*P(2)-P(3)) / (2*dx)) - (f * m(1) * abs(m(1)) * c^2) / (2 *D* (P(1))); % Forward discretization for dm/dt at i=1
dPdt(Nx) =-c^2 * (3*m(Nx) - 4*m(Nx - 1) + m(Nx-2)) / (2*dx); % Backward discretization of dP/dt at i=Nx
% Central discretization for in-between grid points (i=2 to 49)
for i = 2:49
if i == 25
% Skip the calculation for dPdt(25) since it's already calculated outside the loop
%dPdt(25) = ...;
%dmdt(25) =...;
else
% Calculate dP/dt and dm/dt for in-between grid points using central discretization
dPdt(i) = -c^2 * (m(i+1) - m(i-1)) / (2 * dx); % Central discretization of dP/dt
dmdt(i) =-2*m(i)*c^2/P(i)*((m(i + 1) - m(i - 1)) / (2 * dx))+(c^2*m(i)^2/P(i)^2)*(P(i+1) - P(i-1)) / (2*dx) - (P(i+1) - P(i-1)) / (2*dx) - (f * m(i) * abs(m(i))* c^2 ) /( 2*D*P(i)); % Central discretization of dm/dt
end
end
% Combine dPdt and dmdt into a single vector
dudt = [dPdt; dmdt];
end
%% why the graph is bulge at Node 25, is it correct or wrong?
0 Commenti
Risposte (1)
Torsten
il 12 Dic 2023
Choose less output times from the solver, e.g. by setting
tspan = linspace(tcut(i), tcut(i+1),(tcut(i+1)-tcut(i))/10)
instead of
tspan = [tcut(i) tcut(i+1)] ;
2 Commenti
Torsten
il 14 Dic 2023
I told you my opinion about your discretization. I won't try to heal something that I think is not curable.
Vedere anche
Categorie
Scopri di più su Mathematics and Optimization in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!