solving coupled ODE's by ode45

1 次查看(过去 30 天)
Christina Reid
Christina Reid 2021-3-9
评论: darova 2021-3-10
I am trying to write a MATLAB script which solves the Reynold's transport theorem. I have two equations, attached here. I know that these two equations form a coupled ODE system, but I have no idea how to approach it. The two screenshots are shows in the .png files called 'Screen Shot 2021-03-08 at 8.45.43 PM.png' and 'Screen Shot 2021-03-08 at 8.46.58 PM.png'
I wrote a function for a similar problen, called diff_fnc_mb, and this function is used to solve an ODE for one equation. The equation for this function is also attached 'Screen Shot 2021-03-08 at 8.57.48 PM.png' The difference between this equation and the equations that need to be used for the coupled ODE, is the variable T_c is constant in 'Screen Shot 2021-03-08 at 8.57.48 PM.png' while T_c varys for the coupled ODE. Could someone please help me with this
  1 个评论
Christina Reid
Christina Reid 2021-3-9
Below are my defined variables:
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]
temp_sens = 0.0015 % K^-1
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]
T_1 = 285; %[ K ]
T_2 = 315; %[ K ]
% constants
% ========================================================================= %
R = 8314.5; % gas constant [J/kmol)
R_gas = 287; % J/kg K
% 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);

请先登录,再进行评论。

回答(1 个)

darova
darova 2021-3-9
Here is how i see it
dTc = R*Tc/Pc/Vc*(gamma*(mb*Tf-mt*Tc)-(mb-mt)*Tc);
dPc = R*Tc/Vc*(rho*Sb*a*Pn - Pc*At/cstar)
dVc = gamma*R/Vc*(mb*Tf-mt*Tc)-dPc;
dVc = Vc/Pc*dVc;
dY = [dTc dPc dVc]';
  2 个评论
Christina Reid
Christina Reid 2021-3-10
Thank you! So I am assuming when I write the function, the out put of the function would be dY?

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Statics and Dynamics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by