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 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で General Applications についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!