MATLAB Answers

Using Ode45 to solve dynamics problem (ISA model)

31 views (last 30 days)
JXT119
JXT119 on 26 Sep 2021 at 9:20
Commented: JXT119 on 26 Sep 2021 at 14:44
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 Comments
Walter Roberson
Walter Roberson on 26 Sep 2021 at 12:24
What problem are you encountering?
z0 = const.h0;
v0 = const.v0;
we will need to know those in order to test.

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 26 Sep 2021 at 12:25
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 Comment
JXT119
JXT119 on 26 Sep 2021 at 14:44
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!

Sign in to comment.

More Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by