Using Ode45 to solve dynamics problem (ISA model)

2 次查看（过去 30 天）
JXT119
JXT119 2021-9-26
评论： JXT119 2021-9-26
How can I express a variable density in diferentialfuntion2? I have modeled athmospheric density in complex_atm_model in terms of height (which is position in the one dimensioinal motion model) but I am having trouble when trying to solve using ode45. Any ideas?
z0 = const.h0;
v0 = const.v0;
t0 = 0;
tf = 800;
N = 60000;
t = linspace(t0, tf, N);
X = [z0;v0];
[t,X] = ode45('diferentialfunction2', t, X);
------------------------------------------------------- Main code
function dXdt = diferentialfunction2(t, X)
c=2;
s=0.3;
g=9.81;
m = 100;
G = 6.67E-11;
R = 6371E3;
M = 5.9722E24;
h = 39045:-1:0;
h(39046) = 0;
[rho, T] = complex_atm_model(h,X);
dXdt(1,1) = X(2);
dXdt(2,1) = ((rho*c*s*((X(2))^2))/(2*m))-(G*(M/(R+X(1))^2));
end
---------------------------------------------------------------------------------
function [density, T] = complex_atm_model (h,X)
% Parameters to be used:
T_SL = 288.16; % Sea level temperature [K]
rho_SL = 1.225; % Sea level density [kg/m3]
Rg = 287.0; % Gas constant for the air [J/kg].
dT_dh = -0.0065; % Gradient of temperature in the troposphere [K/m]
g0 = 9.81; % Gravity acceleration at sea level [m/s2]
%Temperature model for all altitudes:
T = (T_SL * h ./ h) + (dT_dh*h); %Temperature vector modeled constant
inds = find(h) > 11000; %Find indexes for altitudes > 11000 m
T(inds) = T(11000); %Set constant T = T(11000) for h > 11000
%Density model for all altitudes:
inds_low = find(h < 11000); %Define inds_low as indexes for altitudes < 11000
rho = rho_SL * h./h; %Prelocate a density vector with uniform density
rho(inds_low) = (rho_SL * (T(inds_low)/T_SL).^4.26); %Density vector modeled constant
exp_athmos = (Rg * T(11000))/g0; %Measure how density decays with altitude
for i=11000:39045
rho(i) = rho(11000)*exp(-(h(i)-11000)/exp_athmos); %Set density values for h > 11000
end
density = rho(X(1));
end
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2021-9-26
What problem are you encountering?
z0 = const.h0;
v0 = const.v0;
we will need to know those in order to test.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Alan Stevens
Alan Stevens 2021-9-26
Ran in:
Like this
z0 = 39045; %const.h0;
v0 = 0; %const.v0;
t0 = 0;
tf = 800;
N = 60000;
tspan = linspace(t0, tf, N);
X = [z0;v0];
[t,X] = ode45(@diferentialfunction2, tspan, X);
plot(t,X(:,1)),grid
xlabel('time'), ylabel('height')
function dXdt = diferentialfunction2(~, X)
c=2;
s=0.3;
% g=9.81;
m = 100;
G = 6.67E-11;
R = 6371E3;
M = 5.9722E24;
h = 0:39045;
rho = complex_atm_model(h,X);
dXdt(1,1) = X(2);
dXdt(2,1) = ((rho*c*s*((X(2))^2))/(2*m))-(G*(M/(R+X(1))^2));
end
function density = complex_atm_model (h,X)
% Parameters to be used:
T_SL = 288.16; % Sea level temperature [K]
rho_SL = 1.225; % Sea level density [kg/m3]
Rg = 287.0; % Gas constant for the air [J/kg].
dT_dh = -0.0065; % Gradient of temperature in the troposphere [K/m]
g0 = 9.81; % Gravity acceleration at sea level [m/s2]
%Temperature model for all altitudes:
T = T_SL + dT_dh*h; %Temperature vector modeled constant
inds = 11000:max(h); %Find indexes for altitudes > 11000 m
T(inds) = T(11000); %Set constant T = T(11000) for h > 11000
%Density model for all altitudes:
inds_low = 1:11000; %Define inds_low as indexes for altitudes < 11000
%rho = rho_SL; %Prelocate a density vector with uniform density
rho(inds_low) = (rho_SL * (T(inds_low)/T_SL).^4.26); %Density vector modeled constant
exp_athmos = (Rg * T(11000))/g0; %Measure how density decays with altitude
for i=11000:max(h)+1
rho(i) = rho(11000)*exp(-(h(i)-11000)/exp_athmos); %Set density values for h > 11000
end
density = interp1(h,rho,X(1));
end
  1 个评论
JXT119
JXT119 2021-9-26
Thank you so much!! It worked perfectly. I see my main problem were the indexes in the complex_atm_model, thank you for the clarification!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by