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(1)
plot(t,y)
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,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.1,60], P);
y = P;
figure(2)
plot(t,y)
xlabel("Time (s)")
ylabel("Chamber Pressure (Pa)")
title("Chamber Pressure vs Time ")
function dP = Fun(t,P,R,T0,rp,a,n,t1,Ab,rhoP,Astar,k,v0)
char t1;char rp;
if t==0
rp=0.35;
t1=0;
end
if t>=0.1
rp>=0.7;
Ab=0;
end
v0=pi*rp^2*8; %recalculating free volume available in chamber
t1=t; %storing current time step to be used in next iteration
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
