solving coupled ODE's by ode45

조회 수: 5(최근 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
darova
darova 2021년 3월 10일
Yes, exactly

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

태그

Community Treasure Hunt

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

Start Hunting!

Translated by