The graph is not showing up even though the ODE solver is working...

1 Ansicht (letzte 30 Tage)
StillANovice
StillANovice am 1 Aug. 2018
Erneut geöffnet: Walter Roberson am 22 Dez. 2018
Function:
function dAdT = polyadiabatic(T,A)
dAdT = zeros(2,1);
% Kinetics
Rp = (kp).*rho.*A(1).*((2.*ki.*((rho.*A(1).^3)))./ktc).^0.5;
% Formulae
ki = 2.019.*10.*exp(-13810./A(2)); % m6/kg.s
kp = 1.009.*(10.^5).*exp(-3557./A(2)); % m3/kg.s
%kp = 100900;
ktc = 2.205.*(10.^7).*exp(-844./A(2)).*exp(-2.*((A1.*wp)+(A2.*(wp.^2))+(A3.*(wp.^3)))); % m3/kg.s
A1 = 2.57-((5.05.*0.001).*A(2));
A2 = 9.56-((1.76.*0.01).*A(2));
A3 = -3.032+((7.85.*0.001).*A(2));
rho = 845-(A(2)-353)+((200+(A(2)-353)).*wp); % kg/m3
wp = 1-A(1);
% Pressure
Pv_s = exp(135.23395+((-9226.418).*(A(2).^(-1)))+((-17.99862).*log(A(2)))+(0.01669633.*A(2))); % Vapour pressure, Pa
% Parameters
deltaH = -6.7.*(10.^5); % J/kg
cp = 1.884.*1000; % J/kg.K
% Mass balance
dAdT(1) = -(Rp./rho);
% Energy balance
dAdT(2) = -(Rp.*deltaH)./(cp.*rho);
end
Calling the function:
function [T,A] = call_polyadiabatic()
tspan = [0 452100];
% Initial conditions
A1_0 = 1;
A2_0 = 318;
[T,A] = ode15s(@polyadiabatic,tspan,[A1_0 A2_0]);
yyaxis left
ylabel("Temperature, K")
hold on
grid on
plot(T,A(:,2),'linewidth',1)
yyaxis right
ylabel("Vapour pressure, Pa")
plot(T,Pv_s,'linewidth',1)
hold off
xlabel("Time, s")
title("Adiabatic Temperature and Pressure Profiles")
end
Error:
Undefined function or variable 'Pv_s'.
Error in call_polyadiabatic (line 31)
plot(T,Pv_s,'linewidth',1)
How should I modify my code in order for Pv_s to be plotted?
  5 Kommentare
StillANovice
StillANovice am 30 Aug. 2018
Appreciate your consideration :)

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Star Strider
Star Strider am 1 Aug. 2018
First, I had to rearrange ‘polyadiabatic’ to get it to work, because in the version you posted, several statements were out of order.
This version works:
function dAdT = polyadiabatic(T,A)
dAdT = zeros(2,1);
% Formulae
ki = 2.019.*10.*exp(-13810./A(2)); % m6/kg.s
kp = 1.009.*(10.^5).*exp(-3557./A(2)); % m3/kg.s
%kp = 100900;
A1 = 2.57-((5.05.*0.001).*A(2));
A2 = 9.56-((1.76.*0.01).*A(2));
A3 = -3.032+((7.85.*0.001).*A(2));
wp = 1-A(1);
rho = 845-(A(2)-353)+((200+(A(2)-353)).*wp); % kg/m3
ktc = 2.205.*(10.^7).*exp(-844./A(2)).*exp(-2.*((A1.*wp)+(A2.*(wp.^2))+(A3.*(wp.^3)))); % m3/kg.s
% Kinetics
Rp = (kp).*rho.*A(1).*((2.*ki.*((rho.*A(1).^3)))./ktc).^0.5;
% Pressure
Pv_s = exp(135.23395+((-9226.418).*(A(2).^(-1)))+((-17.99862).*log(A(2)))+(0.01669633.*A(2))); % Vapour pressure, Pa
% Parameters
deltaH = -6.7.*(10.^5); % J/kg
cp = 1.884.*1000; % J/kg.K
% Mass balance
dAdT(1) = -(Rp./rho);
% Energy balance
dAdT(2) = -(Rp.*deltaH)./(cp.*rho);
end
Second, you need to calculate ‘Pv_s’ in ‘call_polyadiabatic’ in order to plot it. Add it (with a minor edit to use the second column of ‘A’) just after your ode15s call:
[T,A] = ode15s(@polyadiabatic,tspan,[A1_0 A2_0]);
Pv_s = exp(135.23395+((-9226.418).*(A(:,2).^(-1)))+((-17.99862).*log(A(:,2)))+(0.01669633.*A(:,2))); % Vapour pressure, Pa
Then the plot works.
Also, it might be more efficient to replace:
exp(135.23395+((-9226.418).*(A(2).^(-1)))
with:
exp(135.23395+((-9226.418)./A(2))
I will let you experiment with that.

Weitere Antworten (0)

Kategorien

Mehr zu Mathematics finden Sie in Help Center und File Exchange

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by