# Using an ODE to solve for an equation

조회 수: 2 (최근 30일)
Christina Reid . 2021년 2월 3일
댓글: darova . 2021년 2월 4일
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]
% 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
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
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표시 없음숨기기 없음
darova 2021년 2월 4일
It's too complicated. Hard to check where the mistake is

댓글을 달려면 로그인하십시오.

### 카테고리

Help CenterFile Exchange에서 General Physics에 대해 자세히 알아보기

### Community Treasure Hunt

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

Start Hunting!