Using an ODE to solve for an equation

6 views (last 30 days)
Christina Reid
Christina Reid on 3 Feb 2021
Commented: darova on 4 Feb 2021
Hello all,
I am looking to solve a differential equation using ODE45. I have gotten pretty far in my script and having an issues that I imagine is really easy to solve. I have attached the equation I am trying to solve for. Labeled pC_eq. I was able to use the ode function to solve for dVc/dt and rho_c/dt ( I have also attached my function I used to calculate that. It is is called test_fnc.m. I am having a hard time calculating dpc/dt (which is the change in chamber pressure over time). dpc/dt is the last line in my main script below. For my assignment I need to calculate the equlibrium chamber pressure and the chamber pressure. I calculated the equlibrium chamber pressure, I labeled it as Pc_eq in my script)
What I am suppose to see is dpc/dt just slightly of a lesser value of Pc_eq, but that is not happening for me in my last equation.
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_0 = 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_p = 1790; % Density [kg/m^3]
rho_c = 1; % Initial chamber density [kg
% 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)
% User inputs - Pressurizing Gas Properties
% ========================================================================= %
p_exit = 1.3e+05; % [Pa]
T_initial = 300; % [K]
% 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]
cap_gamma = sqrt (gamma * (2/(gamma + 1))^((gamma+1)/(gamma-1))); % capital gamma
c_star = (1/cap_gamma)*sqrt((R*T_F)/M);
r_i = d_throat/2; % Need to come back to this value
%L = 10.6; % Need to come back to this value
%beta = a_0* (rho*a_0*c_star*(1/a_t))^(n/(1-n)); %coeff of constants
%t = 0:0.5:bt; % Time vector
%aa = r_i^((2*n-1)/(n-1)); % first part of y eq
%bb = ((2*n-1)/(n-1))* (2*pi*L)^(n/(1-n))*beta*t; % second part of y equation
%y = (aa + bb).^((n-1)/(2*n-1)) - r_i; %burnt web as a function of time
%s_b = 2*pi*(r_i+y)*L; %burnt suface area as a function of y
%K = s_b/a_t; % Klemmung
%Pc_eq = (a_0 *rho*c_star*K).^(1/(1-n));
tow = (V_c/a_t) * (1/(cap_gamma*c_star)); % characteristic time scale
y_0 = 0;
y_f = 1.075;
nn = 100; %sampling rate
%y_result = cell2mat(arrayfun(@(e) linspace(y_0, y_f, nn), y_f', 'UniformOutput', false));
%y_result = y_result';
time_interval = [0,bt];
y0 = 0;
[t_result,y_result] = ode45(@(t,y) diff_fnc(t,y,a_0,n,d_throat, rho_p, c_star), time_interval, y0);
[~,Pc_eq] = chamberPressure(t_result,y_result,a_0,n,d_throat, rho_p, c_star);
% Calculating the Force
x = 1:0.0001:7; % defining a reasonable range for Mach number @ the exit
num = (gamma + 1) /2;
exponent = (gamma+1)/(2*(gamma-1))
dem = (gamma -1)/2;
y = x.*(num./(1+dem*x.^2).^exponent)-0.02; %simplifed equation for my expansion ratio to M_exit
index = zerocross(y)
y(index)
plot(x,y)
hold on
plot(x(index),y(index),'o')
M_e = x(index); % Solving for Mach at the exit
ex_1 = (gamma-1)/gamma ; % I was having the hardest issue with the () so just wrote the exponent to C_f seperately
c_f = cap_gamma.*sqrt( 2*gamma./(gamma-1) * (1 - (p_exit./Pc_eq).^ex_1) ) + (p_exit*e)./Pc_eq;
F_1 = c_f .* a_t .* Pc_eq;
%% plots for motor chamber pressure asa function of both time and web
%%% assuming quassi-steady state (equilibrium)
subplot(1,2,1)
plot(t_result,Pc_eq)
xlabel('time [sec]')
ylabel('Pc_{eq} [Pa]')
title ('(Pc)_{eq} as a Function of Time')
subplot(1,2,2)
plot(y_result,Pc_eq)
xlabel('Burned Web Thickness [m]')
ylabel('Pc_{eq} [Pa]')
title ('(Pc)_{eq} as a Function of Burned Web Thickness')
%% plots for force as a function of both time and web
%%% assuming quassi-steady state (equilibrium)
figure(2)
subplot(1,2,1)
plot(t_result,F_1)
xlabel('time [sec]')
ylabel('Force_{eq} [N]')
title ('Force_{eq} as a Function of Time')
subplot(1,2,2)
plot(y_result,F_1)
xlabel('Burned Web Thickness [m]')
ylabel('Force_{eq{ [N]')
title ('Force_{eq} as a Function of Burned Web Thickness')
%% Step 1.2/1.4
% Find The motor chamber pressure and the vacuum thrust as a function of
% both time and web integrating the ODE representing the overall mass
% balance.
% The variation of the mass of the gases in the chamber over time,
% seperating the role of density variation and that of chamber volume
% variation
%the role of density variation is negligible all the time
% Repreat the calculation of point 2. Assuming a propellant initial
% temperature of 315 (+15 K wrt to baseline) and of 285 K (-15 K wrt
% baseline). Highlight and explain the difference in terms of pressure and
% thrust curves and their integrals over time
for T_ref = [285;300;315] % 3 different reference temperatures [K]
a = a_0*e.^(tow*(T_initial - T_ref));
%Y0 = [8.6; rho_c; 8.6; rho_c; 8.6; rho_c];
%[t_result, Y_result] = ode45(@(t,Y) diff_fnc_mb(t,Y,a,n,a_t, rho_p, c_star, T_ref,R, t_result, y_result), time_interval, Y0);
end
Y0 = [8.6; rho_c]; % initial condition of Vc and rho_c
[t_result_285, Y_result_285] = ode45(@(t,Y) diff_fnc_mb(t,Y,a(1),n,a_t, rho_p, c_star, T_ref(1),R, t_result, y_result), time_interval, Y0);
[t_result_300, Y_result_300] = ode45(@(t,Y) diff_fnc_mb(t,Y,a(2),n,a_t, rho_p, c_star, T_ref(2),R, t_result, y_result), time_interval, Y0);
[t_result_315, Y_result_315] = ode45(@(t,Y) diff_fnc_mb(t,Y,a(3),n,a_t, rho_p, c_star, T_ref(3),R, t_result, y_result), time_interval, Y0);
[dpcdt_300, ~] = ode45(@(t,Y) test_fnc(t,Y,a_0,n,a_t, rho_p, c_star, T_initial,R, t_result, y_result,cap_gamma, V_c), time_interval, Y0);
  1 Comment
darova
darova on 4 Feb 2021
It's too complicated. Hard to check where the mistake is

Sign in to comment.

Answers (0)

Categories

Find more on Quantum Mechanics in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!