# Could someone please help me with my code I have spent a lot of time working on it and I still didn4t figure out what is the problem

1 view (last 30 days)
Amal Fennich on 3 Nov 2019
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)
%at beginning of the integration set initial values for the persistent variables
v0=pi*rp^2*8;
t1=0; %initial time step
dP=@(t,P)Fun(t,P,R,T0,rp,a,n,Ab,P0,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.001], P0);
figure(1)
plot(t,y)
xlabel("Time (s)")
ylabel("Chamber Pressure (Pa)")
title("Chamber Pressure vs Time (Start-Up)")
dP=@(t1,P)Fun(t,P,R,T0,rp,a,n,t1,Ab,P0,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,60], P0);
hold on
figure(2)
plot(t,y)
xlabel("Time (s)")
ylabel("Chamber Pressure (Pa)")
title("Chamber Pressure vs Time ")
hold off
function dP = Fun(t,P,Ab,R,T0,rp,a,n,rhoP,Astar,k)
dP=0;
if t==0
rp=0.35;
end
Ab=2*pi*rp*8;
rhoO=P/(R*T0);
rp>=0.7
Ab=0;
v0=pi*rp^2*8;
t1=t;
rp=min(rp+((a*P^n)*10^-3)*(t-t1),0.7);
Ab=2*pi*rp*8;
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

John D'Errico on 3 Nov 2019
There is no problem that I see. It does exactly what it does. Of course, you give no hint at all as to where the problem may lie, or even what the code should, in your opinion do, or why you think there may be a problem.
To this end, I am now using my MATLAB tarot cards, so see if they identify what is the problem, but all they want to tell me is you may come into some money soon.
Amal Fennich on 3 Nov 2019
my code gives me this error I forgot to include it >> rocket3
Error using rocket3>Fun
Too many input arguments.
Error in rocket3>@(t,P)Fun(t,P,R,T0,rp,a,n,Ab,P0,rhoP,Astar,k,v0) (line 14)
dP=@(t,P)Fun(t,P,R,T0,rp,a,n,Ab,P0,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;
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in rocket3 (line 15)
[t,P]=ode45(dP, [0,0.001], P0);
Stephen Cobeldick on 3 Nov 2019
The function input arguments are inconsistent:
Fun(t,P,R,T0,rp,a,n,Ab,P0,rhoP,Astar,k,v0) % function call
Fun(t,P,Ab,R,T0,rp,a,n,rhoP,Astar,k) % function definition
Aligned on the matching names:
Fun(t,P, R,T0,rp,a,n,Ab,P0,rhoP,Astar,k,v0) % function call
Fun(t,P,Ab,R,T0,rp,a,n, rhoP,Astar,k ) % function definition

Subhadeep Koley on 6 Nov 2019
As rightly pointed by Stephen Cobeldick, the function input arguments are inconsistent. I have made some changes in the script and the function. Check whether it is providing your expected output or not.
SCRIPT
clear;close all;clc;
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)
%at beginning of the integration set initial values for the persistent variables
v0=pi*rp^2*8;
t1=0; %initial time step
dP=@(t,P)Fun(t,P,R,T0,rp,a,n,rhoP,Astar,k);%@(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.001], P0);
figure(1);
plot(t,P);
xlabel("Time (s)");
ylabel("Chamber Pressure (Pa)");
title("Chamber Pressure vs Time (Start-Up)");
dP=@(t,P)Fun(t,P,R,T0,rp,a,n,rhoP,Astar,k);%@(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,60], P0);
hold on;
figure(2);
plot(t,P);
xlabel("Time (s)");
ylabel("Chamber Pressure (Pa)");
title("Chamber Pressure vs Time ");
hold off;
FUNCTION
function dP = Fun(t,P,R,T0,rp,a,n,rhoP,Astar,k)
if t==0
rp=0.35;
end
rhoO=P/(R*T0);
rp>=0.7;
v0=pi*rp^2*8;
t1=t;
rp=min(rp+((a*P^n)*10^-3)*(t-t1),0.7);
Ab=2*pi*rp*8;
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
Hope this helps!