Info

此问题已关闭。 请重新打开它进行编辑或回答。

Can anyone help me fix that error please

1 次查看（过去 30 天）
Mohamed El kholy
Mohamed El kholy 2020-11-19
关闭： MATLAB Answer Bot 2021-8-20
function dydt = vdp1(t,y)
% input parameters
Dr = 0.05;
Dc = 0.05;
Dh = 0.05;
Dt = 0.09;
Db = 0.099;
Lr = 0.5;
Lc = 0.5;
Lh = 0.5;
Li = 0.5;
g = 9.81;
rho = 1000;
mu = 1.002e-3;
Lp = 1;
Dp = 0.0414;
H0 = 1.5;
Hs = 0.025;%
C0 = 10;%
Cs = 10;%
Ct = 10;
Pa = 10^5;
D0 = 0.032;
Ds = 0.025;%
Kv0 = 2;%
Kvs = 2;%
Kp = 2992;
mI2 = 0.3;%
Tr = 30;%
Te = 100;
Tc = 20;
%%
Vr = (pi*(Dr/2)^2)*Lr;
Ac = pi*(Dc/2)^2;
Ah = pi*(Dh/2)^2;
Ap = pi*(Dp/2)^2;
A0 = pi*(D0/2)^2;
As = pi*(Ds/2)^2;
At = pi*(Dt/2)^2;
Ab = pi*(Db/2)^2;
%%
P = y(1);
lh = y(2);
dlh = y(3);
lc = y(4);
dlc = y(5);
lt = y(6);
dlt = y(7);
% Pm equation 59
Pm = ((Ab^(2)/(mI2+mf))+(At/(rho*(lt+Ct)))+(Ac/(rho*lc))+(Ah/(rho*lh)))^(-1)...
*((-Ab/(mI2+mf))*((rho*g*(Li+xp)-Pa)*Ab + Cf*dxp+Hf+Cfd*dxp^(2)-Kp*xp)-...
(At/(rho*lt))*(-Pa-rho*g*lt-8*mu*pi*(lt-Ct)*dlt/At)-(Ac/(rho*lc))*(-P-rho*g*lc-8*mu*pi*lc*dlc/Ac)...
-(Ah/(rho*lh))*(P-rho*g*lh-8*pi*mu*lh*dlh/Ah));
%%
dydt = zeros(9,1);
Vc = Ac*(Lc - lc);
Ve = Ah*(Lh - lh);
dydt(1) = (-P/((Vr/Tr)+(Vc/Tc)+(Ve/Te)))*((1/Tc)*(-Ac*dlc)+(1/Te)*(-Ah*dlh));
dydt(2) = y(3);
dydt(3) = (1/(rho*lh))*((Pm - P) - rho*g*lh - 8*pi*mu*(lh/Ah)*dlh);
dydt(4) = y(5);
dydt(5) = (1/(rho*lc))*((Pm-P)-rho*g*lc-8*pi*mu*lc*dlc/Ac);
dydt(6) = y(7);
dydt(7) = (1/(rho*(lt+Ct)))*((Pm-Pa)-rho*g*lt-8*mu*pi*(lt+Ct)*dlt/At);
end
close all;
clear ;
tspan = 0:0.01:20;
y0 = [1.1025e05 0.5 0 0.4 0 1.5 0 0.04 0]' ;
[t,y] = ode45(@vdp1, tspan, y0);
%%
figure(1)
plot(tspan,y(:,1),'b','LineWidth',2)
title('Variations of working gas pressure')
xlabel('time (Sec)')
ylabel('P (pa)')
%%
figure(2)
plot(tspan,y(:,2),'b','LineWidth',2)
title('Hot liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lh')
%%
figure(3)
plot(tspan,y(:,4),'b','LineWidth',2)
title('Cold liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lc')
%%
figure(4)
plot(tspan,y(:,6),'b','LineWidth',2)
title('Tuning column oscillations')
xlabel('time (Sec)')
ylabel('lc')

此问题已关闭。

回答（2 个）

Stephan
Stephan 2020-11-20
You missed to define several variables, which are used in your equations:
% several used variables are not defined! - i set them all = 1:
mf = 1;
xp = 1;
Cf = 1;
dxp = 1;
Hf = 1;
Cfd = 1;
The complete and working code is here - just change the values for the variables as needed:
tspan = 0:0.01:20;
y0 = [1.1025e05 0.5 0 0.4 0 1.5 0 0.04 0]' ;
[t,y] = ode45(@vdp1, tspan, y0);
%%
figure(1)
plot(tspan,y(:,1),'b','LineWidth',2)
title('Variations of working gas pressure')
xlabel('time (Sec)')
ylabel('P (pa)')
%%
figure(2)
plot(tspan,y(:,2),'b','LineWidth',2)
title('Hot liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lh')
%%
figure(3)
plot(tspan,y(:,4),'b','LineWidth',2)
title('Cold liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lc')
%%
figure(4)
plot(tspan,y(:,6),'b','LineWidth',2)
title('Tuning column oscillations')
xlabel('time (Sec)')
ylabel('lc')
function dydt = vdp1(~,y)
% several used variables are not defined! - i set them all = 1:
mf = 1;
xp = 1;
Cf = 1;
dxp = 1;
Hf = 1;
Cfd = 1;
% input parameters
Dr = 0.05;
Dc = 0.05;
Dh = 0.05;
Dt = 0.09;
Db = 0.099;
Lr = 0.5;
Lc = 0.5;
Lh = 0.5;
Li = 0.5;
g = 9.81;
rho = 1000;
mu = 1.002e-3;
Lp = 1;
Dp = 0.0414;
H0 = 1.5;
Hs = 0.025;%
C0 = 10;%
Cs = 10;%
Ct = 10;
Pa = 10^5;
D0 = 0.032;
Ds = 0.025;%
Kv0 = 2;%
Kvs = 2;%
Kp = 2992;
mI2 = 0.3;%
Tr = 30;%
Te = 100;
Tc = 20;
%%
Vr = (pi*(Dr/2)^2)*Lr;
Ac = pi*(Dc/2)^2;
Ah = pi*(Dh/2)^2;
Ap = pi*(Dp/2)^2;
A0 = pi*(D0/2)^2;
As = pi*(Ds/2)^2;
At = pi*(Dt/2)^2;
Ab = pi*(Db/2)^2;
%%
P = y(1);
lh = y(2);
dlh = y(3);
lc = y(4);
dlc = y(5);
lt = y(6);
dlt = y(7);
% Pm equation 59
Pm = ((Ab^(2)/(mI2+mf))+(At/(rho*(lt+Ct)))+(Ac/(rho*lc))+(Ah/(rho*lh)))^(-1)...
*((-Ab/(mI2+mf))*((rho*g*(Li+xp)-Pa)*Ab + Cf*dxp+Hf+Cfd*dxp^(2)-Kp*xp)-...
(At/(rho*lt))*(-Pa-rho*g*lt-8*mu*pi*(lt-Ct)*dlt/At)-(Ac/(rho*lc))*(-P-rho*g*lc-8*mu*pi*lc*dlc/Ac)...
-(Ah/(rho*lh))*(P-rho*g*lh-8*pi*mu*lh*dlh/Ah));
%%
dydt = zeros(9,1);
Vc = Ac*(Lc - lc);
Ve = Ah*(Lh - lh);
dydt(1) = (-P/((Vr/Tr)+(Vc/Tc)+(Ve/Te)))*((1/Tc)*(-Ac*dlc)+(1/Te)*(-Ah*dlh));
dydt(2) = y(3);
dydt(3) = (1/(rho*lh))*((Pm - P) - rho*g*lh - 8*pi*mu*(lh/Ah)*dlh);
dydt(4) = y(5);
dydt(5) = (1/(rho*lc))*((Pm-P)-rho*g*lc-8*pi*mu*lc*dlc/Ac);
dydt(6) = y(7);
dydt(7) = (1/(rho*(lt+Ct)))*((Pm-Pa)-rho*g*lt-8*mu*pi*(lt+Ct)*dlt/At);
end
  3 个评论
显示 1更早的评论
Mohamed El kholy
Mohamed El kholy 2020-11-20
Also, please check this corrected one after remving the variables that not used.
Your help is greatly apprecaietd regadring that matter !
tspan = 0:0.01:20;
y0 = [1.1025e05 0.5 0 0.4 0 1.5 0 0.04 0]' ;
[t,y] = ode45(@vdp1, tspan, y0);
%%
figure(1)
plot(tspan,y(:,1),'b','LineWidth',2)
title('Variations of working gas pressure')
xlabel('time (Sec)')
ylabel('P (pa)')
%%
figure(2)
plot(tspan,y(:,2),'b','LineWidth',2)
title('Hot liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lh')
%%
figure(3)
plot(tspan,y(:,4),'b','LineWidth',2)
title('Cold liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lc')
%%
figure(4)
plot(tspan,y(:,6),'b','LineWidth',2)
title('Tuning column oscillations')
xlabel('time (Sec)')
ylabel('lc')
function dydt = vdp1(~,y)
% input parameters
Dr = 0.05;
Dc = 0.05;
Dh = 0.05;
Dt = 0.09;
Lr = 0.5;
Lc = 0.5;
Lh = 0.5;
g = 9.81;
rho = 1000;
mu = 1.002e-3;
Pa = 10^5;
Tr = 30;%
Te = 100;
Tc = 20;
%%
Vr = (pi*(Dr/2)^2)*Lr;
Ac = pi*(Dc/2)^2;
Ah = pi*(Dh/2)^2;
At = pi*(Dt/2)^2;
%%
P = y(1);
lh = y(2);
dlh = y(3);
lc = y(4);
dlc = y(5);
lt = y(6);
dlt = y(7);
% Pm equation 59
Pm = (-(At/(rho*lt))*(-Pa-rho*g*lt-8*mu*pi*(lt-Ct)*dlt/At)-(Ac/(rho*lc))*(-P-rho*g*lc-8*mu*pi*lc*dlc/Ac)...
-(Ah/(rho*lh))*(P-rho*g*lh-8*pi*mu*lh*dlh/Ah));
%%
dydt = zeros(9,1);
Vc = Ac*(Lc - lc);
Ve = Ah*(Lh - lh);
dydt(1) = (-P/((Vr/Tr)+(Vc/Tc)+(Ve/Te)))*((1/Tc)*(-Ac*dlc)+(1/Te)*(-Ah*dlh));
dydt(2) = y(3);
dydt(3) = (1/(rho*lh))*((Pm - P) - rho*g*lh - 8*pi*mu*(lh/Ah)*dlh);
dydt(4) = y(5);
dydt(5) = (1/(rho*lc))*((Pm-P)-rho*g*lc-8*pi*mu*lc*dlc/Ac);
dydt(6) = y(7);
dydt(7) = (1/(rho*(lt+Ct)))*((Pm-Pa)-rho*g*lt-8*mu*pi*(lt+Ct)*dlt/At);
end
Stephan
Stephan 2020-11-20
Your system appears to be highly stiff. I rechanged the formula for Pm to the long version, because i was not able to get a solution with the one you used. Note that this means you need parameter values for the additional parameters. Maybe the problem is more easy to solve when meaningful parameters are used. Also i had to reduce the 'RelTol' Parameter extremly and used a stiff solver. Also note that i have changed tspan to a small value - have a look at what happens when you set the upper limit to 20. Since i have no insight in your problem, i guess there is not much more i can do for you.
tspan = [0 0.5];
y0 = [1.1025e05, 0.5, 0, 0.4, 0, 1.5, 0, 0.04, 0]';
opts = odeset('RelTol',1e-2);
[t,y] = ode15s(@vdp1, tspan, y0, opts);
%%
figure(1)
plot(t,y(:,1),'b','LineWidth',2)
title('Variations of working gas pressure')
xlabel('time (Sec)')
ylabel('P (pa)')
%%
figure(2)
plot(t,y(:,2),'b','LineWidth',2)
title('Hot liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lh')
%%
figure(3)
plot(t,y(:,4),'b','LineWidth',2)
title('Cold liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lc')
%%
figure(4)
plot(t,y(:,6),'b','LineWidth',2)
title('Tuning column oscillations')
xlabel('time (Sec)')
ylabel('lc')
function dydt = vdp1(~,y)
% Missing parameters:
mf = 1;
% input parameters
Dr = 0.05;
Dc = 0.05;
Dh = 0.05;
Dt = 0.09;
Lr = 0.5;
Lc = 0.5;
Lh = 0.5;
g = 9.81;
rho = 1000;
mu = 1.002e-3;
Pa = 10^5;
Tr = 30;%
Te = 100;
Tc = 20;
%%
% Dp = 0.0414;
% D0 = 0.032;
% Ds = 0.025;%
Ct = 10;
Db = 0.099;
mI2 = 0.3;%
Li = 0.5;
xp = 1;
Cf = 1;
dxp = 1;
Hf = 1;
Cfd = 1;
Kp = 2992;
Vr = (pi*(Dr/2)^2)*Lr;
Ac = pi*(Dc/2)^2;
Ah = pi*(Dh/2)^2;
% Ap = pi*(Dp/2)^2;
% A0 = pi*(D0/2)^2;
% As = pi*(Ds/2)^2;
% At = pi*(Dt/2)^2;
Ab = pi*(Db/2)^2;
At = pi*(Dt/2)^2;
%%
P = y(1);
lh = y(2);
dlh = y(3);
lc = y(4);
dlc = y(5);
lt = y(6);
dlt = y(7);
% Pm equation 59
Pm = ((Ab^(2)/(mI2+mf))+(At/(rho*(lt+Ct)))+(Ac/(rho*lc))+(Ah/(rho*lh)))^(-1)...
*((-Ab/(mI2+mf))*((rho*g*(Li+xp)-Pa)*Ab + Cf*dxp+Hf+Cfd*dxp^(2)-Kp*xp)-...
(At/(rho*lt))*(-Pa-rho*g*lt-8*mu*pi*(lt-Ct)*dlt/At)-(Ac/(rho*lc))*(-P-rho*g*lc-8*mu*pi*lc*dlc/Ac)...
-(Ah/(rho*lh))*(P-rho*g*lh-8*pi*mu*lh*dlh/Ah));
%%
dydt = zeros(9,1);
Vc = Ac*(Lc - lc);
Ve = Ah*(Lh - lh);
dydt(1) = (-P/((Vr/Tr)+(Vc/Tc)+(Ve/Te)))*((1/Tc)*(-Ac*dlc)+(1/Te)*(-Ah*dlh));
dydt(2) = y(3);
dydt(3) = (1/(rho*lh))*((Pm - P) - rho*g*lh - 8*pi*mu*(lh/Ah)*dlh);
dydt(4) = y(5);
dydt(5) = (1/(rho*lc))*((Pm-P)-rho*g*lc-8*pi*mu*lc*dlc/Ac);
dydt(6) = y(7);
dydt(7) = (1/(rho*(lt+Ct)))*((Pm-Pa)-rho*g*lt-8*mu*pi*(lt+Ct)*dlt/At);
end
Mohamed El kholy
Mohamed El kholy 2020-11-20
Thanks so much for your help. Highly appreciated!
I will run it and see.

此问题已关闭。

标签

Community Treasure Hunt

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

Start Hunting!

Translated by