Using Ode45 to solve dynamics problem (ISA model)

조회 수: 1(최근 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