Ode15s , Plotting issue.
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
Torsten
le 12 Déc 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 commentaires
Torsten
le 14 Déc 2023
I told you my opinion about your discretization. I won't try to heal something that I think is not curable.
Voir également
Catégories
En savoir plus sur Mathematics and Optimization dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!