Error. I want to plot this ode45 pressure versus time

1 Ansicht (letzte 30 Tage)
R=8.314/32;
T0=2930;
a=(12)/((2.027*10^6)^0.45);
rhoP=1920;
Astar=pi*0.25^2;
k=1.35;
n=0.45;
P0=101325;
syms P(t)
for t = [0,0.1]
dP=@(P,t)(Ab*a*P^n*(rhoP-rhoO)-P*Astar*sqrt(k/(R*T0))*(2/(k+1))^((k+1)/(2*(k-1))))*R*T0/v0;
[P,t]=ode45(dP, [0,0,1], P0);
end
if t==0 %at beginning of the integration set initial values for the persistent variables
rp=0.35; %initial port radius
t1=0; %initial time step
end
Ab=2*pi*rp*8;%burn area
rhoO=P/(R*T0); %gas density
rp=min(rp+((a*P^n)*10^-3)*(t-t1),0.7);
figure(2)
plot(t,y)
xlabel("Time (s)")
ylabel("Chamber Pressure (Pa)")
title("Chamber Pressure vs Time (Start-Up)")

Akzeptierte Antwort

Thiago Henrique Gomes Lobato
Bearbeitet: Thiago Henrique Gomes Lobato am 3 Nov. 2019
There were some errors in your code:
  • For ode45 you give numerical functions, not symbolic ones
  • There were some variables that were used in the function but initialized after the loop
  • If you have time dependent variables, you ideally should define them inside the function
  • Your variable order in the function was inverted
  • You didn't declared v0 anywhere
I choosed a random value for v0 and fix the other issues of your code, you just have to make sure that the considerations you made for your function and variables are right to fully trust the results.
R=8.314/32;
T0=2930;
a=(12)/((2.027*10^6)^0.45);
rhoP=1920;
Astar=pi*0.25^2;
k=1.35;
n=0.45;
P0=101325;
P = P0;
%syms P(t)
%at beginning of the integration set initial values for the persistent variables
rp=0.35; %initial port radius
t1=0; %initial time step
rhoO=P/(R*T0); %gas density
rp=min(rp+((a*P^n)*10^-3)*(t1-t1),0.7);
Ab=2*pi*rp*8;%burn area
v0 = 0.1;
dP=@(t,P)Fun(t,P,R,T0,rp,a,n,t1,Ab,rhoP,Astar,k,v0);%@(t,P)(Ab.*a.*P.^n.*(rhoP-rhoO)-P.*Astar.*sqrt(k/(R.*T0)).*(2/(k+1)).^((k+1)/(2.*(k-1)))).*R.*T0./v0;
[t,P]=ode45(dP, [0,0.1], P);
y = P;
figure(2)
plot(t,y)
xlabel("Time (s)")
ylabel("Chamber Pressure (Pa)")
title("Chamber Pressure vs Time (Start-Up)")
function dP = Fun(t,P,R,T0,rp,a,n,t1,Ab,rhoP,Astar,k,v0)
rhoO=P/(R*T0); %gas density
rp=min(rp+((a*P^n)*10^-3)*(t-t1),0.7);
Ab=2*pi*rp*8;%burn area
dP = (Ab.*a.*P.^n.*(rhoP-rhoO)-P.*Astar.*sqrt(k/(R.*T0)).*(2/(k+1)).^((k+1)/(2.*(k-1)))).*R.*T0./v0;
end

Weitere Antworten (1)

Ous Chkiri
Ous Chkiri am 3 Nov. 2019
if i want to add this condition if rp>=0.7 %if grain gets exhausted
Ab=0; %burn area =0 and plotting from 0 to 60 but from 0>>0.1 we apply rp 0.3 after we have 0.7
[0,0.1,60]
  2 Kommentare
Thiago Henrique Gomes Lobato
If you have conditions that make a discontinuity you will have to perform the integration in a piece-wise matter and add the conditions in your integration function ("Fun" in the example I gave you). This answer and comments shows an example about how you could do it and also possible pitfalls https://de.mathworks.com/matlabcentral/answers/487643-adding-a-piecewise-time-dependent-term-in-system-of-differential-equation#answer_398394?s_tid=prof_contriblnk
Ous Chkiri
Ous Chkiri am 3 Nov. 2019
I did not understand to do it and i already see the link

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by