Hi, is it possible to have variables in an ode45 that varies according to another function?

4 Ansichten (letzte 30 Tage)
This is my ode45 that I have to solve and, actually, the variables rho and speedsound (i have commented them out in the first script) is supposed to vary with altitude (x(1)). How do I allow this to work if I have another function called atmos and I want it to go back to atmos to retrieve the values of rho and speed sound for ever altitude it is before the next time step of integration?
%% Dynamic equations to do ODE45
function xdot = dynameqn(t,x,eledefl, thrustvec)
xdot = zeros(1,7);
g = 9.81;
c = 1.74;
S = 17.1;
initial_mass = 1248.5;
Ix = 1421;
Iy = 4067.5;
Iz = 4786;
Ixz = 200;
thrust_sfc = 200/60/60/1000; %kg/kN/h --> kg/N/s --> /60^2 & /1000
% rho = 1.225;
% speedsound = 340.26;
powersetting = 4;
height = x(1);
horz_displacement = x(2);
alpha = x(3);
u = x(4);
q = x(5);
theta = x(6);
mass = x(7);
CL = 6.44*alpha + 3.8*c*q/(2*u) + 0.355*eledefl; % lift coefficient
L = 0.5*rho*u^2 * CL * S; % lift force
CD = 0.03 + 0.05*CL^2; % drag coefficient
D = 0.5*rho*u^2 * CD* S; % drag force
Cm = 0.05 - 0.683*alpha - 9.96*c*q/(2*u) - 0.923*eledefl; % pitching moment coefficient
m = 0.5*rho*u^2*S*c*Cm; % pitching moment
thrust = powersetting * ( (7+u/speedsound)*200/3 + height/1000*(2*u/speedsound -11) ); % thrust force
xdot(1) = u*sin(theta-alpha); % h'= rate of vertical climb
xdot(2) = u*cos(theta-alpha); % x'= rate of change of horizontal displacement
xdot(3) = q - (thrust*sin(alpha+thrustvec) + L) / (mass* u) + g/u*cos(theta - alpha); % alpha'
xdot(4) = (thrust * cos(alpha+thrustvec) - D)/ mass - g*sin(theta-alpha); % u'
xdot(5) = m/Iy; % q'
xdot(6) = q; % theta'
xdot(7) = -thrust_sfc * thrust; % mass'
xdot=xdot';
end
atmos script below:
% Calculating T, p ,rho and speedsound for every altitude in the ISA atmosphere
function [rho, speedsound] = atmos_ode45(h)
h1 = 11; % Height of tropopause
h2 = 20; % End height of table
g = 9.81;
R = 287;
c = 6.51e-3; % temperature lapse dt/dh = - c = -6.51 degcelcius/km
T0 = 15+273.15; % Temperature sea level
p0 = 101325; % pressure sealevel
rho0 = 101325/R/T0; % density sealevel = pressure / R*T, R=287, T = 15 degcelcius
T1 = T0 - c*h1*1000; % Temperature at 11km
p1 = p0 * (T1/T0)^5.2506; % Pressure at 11km
rho1 = rho0 * (T1/T0)^4.2506; % Density at 11km
T2 = T1; % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1)*1000); % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1)*1000); % Density at 20km
if h <= h1
% disp('Troposphere');
T = T0 - c*h*1000
p = p0 * (T/T0)^5.2506
rho = rho0 * (T/T0)^4.2506
speedsound = (1.4*R*T)^0.5
elseif h <= h2
% disp('Tropopause');
T = T1
p = p1 * exp(-g/(R*T)*(h-h1)*1000)
rho = rho1 * exp(-g/(R*T)*(h-h1)*1000)
speedsound = (1.4*R*T)^0.5
end
return
end

Antworten (1)

Alan Stevens
Alan Stevens am 12 Jan. 2022
Put something like
[T, p, rho, speedsound] = atmos(height);
immediately after
height = x(1);
in function dynameqn.
Not sure why you return T and p from atmos, as you don't seem to use them in dynameqn. Perhaps you use them elsewhere.
  7 Kommentare
Torsten
Torsten am 12 Jan. 2022
If this is unphysical, you should check your equations for errors.
Stephen23
Stephen23 am 12 Jan. 2022
"Do you happen to have any idea how I can fix this? "
Make sure that continuous values are defined for all h.

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Dates and Time finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by