Using Ode45 to solve dynamics problem (ISA model)

1 Ansicht (letzte 30 Tage)
JXT119
JXT119 am 26 Sep. 2021
Kommentiert: JXT119 am 26 Sep. 2021
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 Kommentare
JXT119
JXT119 am 26 Sep. 2021
Thanks for your comment, just did it.
Walter Roberson
Walter Roberson am 26 Sep. 2021
What problem are you encountering?
z0 = const.h0;
v0 = const.v0;
we will need to know those in order to test.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Alan Stevens
Alan Stevens am 26 Sep. 2021
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 Kommentar
JXT119
JXT119 am 26 Sep. 2021
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!

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu MATLAB finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by