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) odeargumen​ts(FcnHand​lesUsed, solver_name, ode, tspan, y0, options, varargin);

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

[Walter Roberson]

The code runs without error for me.

[Ous Chkiri]

Yes. Thank you