how to plot the motor chamber pressure and the vacuum thrust as a function of both time and web assuming quasi-steady state (equilibrium)

I have computed the equilibrium pressure assuming a quasi steady state condition: using the following eq attached to this question and below is my code. How do I plot Pc_eq as a function of time?

clear all
%% Step 1
% ========================================================================= %
% Basing on the motor and propellant data indicated in Table 1 and on the 
% burning surface evolution shown in Figs. 5 and 6 and assuming a constant
% chambre temperature profile compute the following curves without taking
% into account the non-ideal parameters.
%% Step 1.1 
%Find: The motor chamber pressure and the vaccume thrust as a
% function of  both time and web assuming quasi-steady state (equilibrium)
% User inputs - First stage P80 SRM
% ========================================================================= %
mp = 88000; % Propellant Mass [kg]
ms = 7330;  % Structural mass [kg]
l = 10.6;   % length [m]
d = 3.0;    % Diameter [m]
f = 3015;   % Max thrust(vaccum) [kN]
bt = 110;   % Buring time [sec]
isp = 280;  % Specifi Impulse (vaccuum) [sec]


% User inputs - Propellant Ballistic Properties 
% ========================================================================= %
a = 1.847e-05; % Temperature coefficient @ 300 K [m/s * Pae-0.4]
n = 0.4;       % Combustion index
tau = 0.0015;  % Temperature sensitivity [k^-1]
rho = 1790;    % Density [kg/m^3]

% User inputs - Propellant thermochemocal properties 
% ========================================================================= %
T_F = 3550;   % Flame temperature [K]
M = 29;       % Molecular mass [kg/kmole]
gamma = 1.13; % Specific heat ratio

% User inputs - Motor geometrical properties 
% ========================================================================= %
d_throat = 0.496; % Throat diameter [m]
e = 16;           % expansion ratio
v_c = 8.6;        % Initial chamber volume [m^3]
v_frac = 0.85;    % volumetric loading fraction (V_c/(V_c+V_p)
% constants 
% ========================================================================= %
r = 8314.5       % gas constant [J/kmol)
% Calculations 
% Find: The motor chamber pressure and the vaccum thrust as a
% ========================================================================= %
a_t = pi* (d_throat/2)^2; % Thrat area [m^2]
m_dot = mp/bt;                 % mass flow rate [kg/s]
s_b = 48;                      % burning surface area - pulled from the Sb vs Y plot [m^2]
cap_gamma = sqrt (gamma * (2/(gamma + 1))^((gamma+1)/(gamma-1))); % capital gamma
c_star = (1/cap_gamma)*sqrt((r*T_F)/M);
K = s_b/a_t;              % Klemmung
Pc_eq = (a *rho*c_star*K)^(1/(1-n));


 

[Mischa Kim]

Hi Christina, I do not see any time dependency of Pc_eq in your code or the equation (screen shot). But I guess at least one of the variables in Pc_eq is time dependent. Can you share some more info about this?

[Christina Reid]

Hi Mischa,

The source term in the Pc_eq equation is the s_b (burning surface area) which lives in the K (the klemmung). So s_b is a fucntion of time and decreasing as as the burn time increases

[Christina Reid]

Hi again Mishca,

I actually added an attachment which explains PC_eq dependecy. The variable,y, being the web thickness = 1.09.

Answer 1

[Gaurav Garg]

Hi Christina,

The variable Pc_eq is a scalar value, in your case. You need to have a vector in order to plot it.

Moreover, your equation doesn't seem to have any time dependency.

So, I would suggest you to rewrite the code and give the variables you want to plot, some dependency over time and then use plot funciton.