# solving coupled ODE's by ode45

조회 수: 5(최근 30일)
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 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 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표시숨기기 이전 댓글 수: 1
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!