Azzera filtri
Azzera filtri

solving coupled ODE's by ode45

1 visualizzazione (ultimi 30 giorni)
Christina Reid
Christina Reid il 9 Mar 2021
Commentato: darova il 10 Mar 2021
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 Commento
Christina Reid
Christina Reid il 9 Mar 2021
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);

Accedi per commentare.

Risposte (1)

darova
darova il 9 Mar 2021
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 Commenti
Christina Reid
Christina Reid il 10 Mar 2021
Thank you! So I am assuming when I write the function, the out put of the function would be dY?
darova
darova il 10 Mar 2021
Yes, exactly

Accedi per commentare.

Categorie

Scopri di più su Statics and Dynamics in Help Center e File Exchange

Tag

Community Treasure Hunt

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

Start Hunting!

Translated by