Ode15s , Plotting issue.

Question

MAYUR MARUTI DHIRDE 2023-12-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2059829-ode15s-plotting-issue

评论： Torsten 2023-12-14

在 MATLAB Online 中打开

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

%% 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

% 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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-12-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2059829-ode15s-plotting-issue#answer_1370329

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

MAYUR MARUTI DHIRDE 2023-12-14

I have a question regarding numerical stability which can be seen in the graphs, At node 25 I am using backward and forward discretization. For in between point from 2 to 49 I am using central discretization. These varying discretization is causing a stability issue. Is there any solution to solve this problem.