Using Ode45 to solve dynamics problem (ISA model)

3 ビュー (過去 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 件のコメント
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 日
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 件)

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by