Error in @(t,y)dPress(t,y,params) 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);
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
clc, clear, close %clearing workspace, closing any open plots
params=struct; %this is a variable for passing parameters into the function which returns dP/dt
params.R=8.314/32; %getting R
params.T0=2930; %ignition temperature
params.a=(12)/((2.027*10^6)^0.45); %getting a from given information in question
params.rhoP=1920; %density of propellant
params.Astar=pi*0.25^2; %throat area
params.k=1.35; %ratio of specific heats
params.n=0.45; %burn exponent
P0=101325; %initial pressure=atmospheric pressure
[t,y]=ode45(@(t,y) dPress(t,y,params),[0,60],P0); %numerically integrates dP/dt and returns time step and the pressure at that time step
figure(1)
plot(t,y)
xlabel("Time (s)")
ylabel("Chamber Pressure (Pa)")
title("Chamber Pressure vs Time")
%redoing with timespan [0,0.1] to get a closer look at start up phase
[t,y]=ode45(@(t,y) dPress(t,y,params),[0,0.1],P0); %numerically integrates dP/dt and returns time step and the pressure at that time step
figure(2)
plot(t,y)
xlabel("Time (s)")
ylabel("Chamber Pressure (Pa)")
title("Chamber Pressure vs Time (Start-Up)")
function dP=dPress(t,P,params)
%this function calculates and returns dP/dt given time t, current pressure
%P and parameters in params
dP=0; %This variable stores dP/dt
persistent t1 rp %we need to remember these variables from last iteration and recalculate these in the current iteration
%v0 is the current free volume, t1 stores last time step, rp stores current
%port radius
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
%extracting parameters from params
R=params.R;
T0=params.T0;
a=params.a;
rhoP=params.rhoP;
Astar=params.Astar;
k=params.k;
n=params.n;
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); %recalculating port radius, port radius will increase by (regression rate)*(time elapsed from last iteration)
%rp cannot be more than 0.7, therefore minimum of two values used.
if rp>=0.7 %if grain gets exhausted
Ab=0; %burn area =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
dP=(Ab*a*P^n*(rhoP-rhoO)-P*Astar*sqrt(k/(R*T0))*(2/(k+1))^((k+1)/(2*(k-1))))*R*T0/v0; %calculating dp/dt
end
2 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu General Applications finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!