How to implement this coupled ODE

3 vues (au cours des 30 derniers jours)
Romio
Romio le 8 Déc 2022
Commenté : Sam Chak le 8 Déc 2022
I want to solve this ODE but I don't know how to implement it. Basically I want du/dt to be some input function that I determine (e.g. heviside or constant function) and use its integral elsewhere, rho is obtained from some data fitting that is dependent on h = u * t.
Here is my try. I don't know how to define du/dt as an input, and at the same time use its solution to calcualte h and dm/dt
clc,clear,close all
global Ve Cd A_r g
g = 9.8; % m/s^2
Cd = 0.47;
r = 3.7; % m
A_r = pi * r^2;
rho_fuel = 1000; % Kg/m^3
Ve = 3000; % m/s
Mbody = 60684; % Kg
Mfuel = 488370; % Kg
% Mfuel = Mfuel + 450000;
Mpayload = 31184; % Kg
t0 = 0;
Tf = 540;
m0 = Mbody + Mfuel + Mpayload;
a0 = 14;
% constant acceleration input
a_in = 14*ones(1,Tf);
tspan = [t0 Tf];
[t,m] = ode45(@(t,m) myode(t,m), tspan, m0);
function dydt = myode(t,y,a_in)
global Ve Cd A_r g
% dydt(2) = y(2);
dydt(2) = interp1(a_in,t);
H=[-1;0;1;2;3;4;5;6;7;8;9;10;11;13;15;17;20;25;30;32;35;40;45;47;50;51;60;70;71;80;84.9;89.7;100.4;105;110;120];
H=H*1000;
rho=[1.347;1.225;1.1116;1.0065;0.9091;0.8191;0.7361;0.6597;0.5895;0.5252;0.4664;0.4127;0.3639;0.2655;0.1937;0.1423;0.0880;0.0395;0.0180;0.0132;0.0082;0.0039;0.0019;0.0014;0.0010;0.00086;0.000288;0.000074;0.000064;0.000015;0.000007;0.000003;0.0000005;0.0000002;0.0000001;0.0000001];
f_rho=fit(H, rho, 'linearinterp');
h = y(2) * t; % what is y(2)??
rho_air = f_rho(h);
Fd = 0.5 * rho_air * (y(2))^2 * A_r * Cd; % what is y(2)??
dydt(1) = (1/Ve) * (-m * g - Fd - m * dudt(2));
end

Réponses (2)

Torsten
Torsten le 8 Déc 2022
Modifié(e) : Torsten le 8 Déc 2022
clc,clear,close all
global Ve Cd A_r g
g = 9.8; % m/s^2
Cd = 0.47;
r = 3.7; % m
A_r = pi * r^2;
rho_fuel = 1000; % Kg/m^3
Ve = 3000; % m/s
Mbody = 60684; % Kg
Mfuel = 488370; % Kg
% Mfuel = Mfuel + 450000;
Mpayload = 31184; % Kg
t0 = 0;
Tf = 540;
m0 = Mbody + Mfuel + Mpayload;
a0 = 14;
% constant acceleration input
a_in = 14*ones(1,Tf+1);
tspan = [t0 Tf];
H=[-1;0;1;2;3;4;5;6;7;8;9;10;11;13;15;17;20;25;30;32;35;40;45;47;50;51;60;70;71;80;84.9;89.7;100.4;105;110;120];
H=H*1000;
rho=[1.347;1.225;1.1116;1.0065;0.9091;0.8191;0.7361;0.6597;0.5895;0.5252;0.4664;0.4127;0.3639;0.2655;0.1937;0.1423;0.0880;0.0395;0.0180;0.0132;0.0082;0.0039;0.0019;0.0014;0.0010;0.00086;0.000288;0.000074;0.000064;0.000015;0.000007;0.000003;0.0000005;0.0000002;0.0000001;0.0000001];
f_rho=fit(H, rho, 'linearinterp');
[t,y] = ode45(@(t,m) myode(t,m,a_in,f_rho), tspan, [m0,a0]);
figure(1)
plot(t,y(:,1))
figure(2)
plot(t,y(:,2))
function dydt = myode(t,y,a_in,f_rho)
global Ve Cd A_r g
m = y(1);
u = y(2);
dudt = interp1(0:540,a_in,t);
h = u * t;
rho_air = f_rho(h);
Fd = 0.5 * rho_air * u^2 * A_r * Cd;
dmdt = (1/Ve) * (-m * g - Fd - m * dudt);
dydt = [dmdt;dudt];
end
  1 commentaire
Sam Chak
Sam Chak le 8 Déc 2022
I think the final mass m should be Mbody + Mpayload, after the fuel is exhausted.
The velocity of the rocket when it runs out of fuel should be , not increasing continuously.
The altitude of the rocket is also unrealistically high.
I have not fixed the velocity. But I think it has something to do with the acceleration.
global Ve Cd A_r g Tf Mbody Mpayload
g = 9.8; % m/s^2
Cd = 0.47;
r = 3.7; % m
A_r = pi * r^2;
rho_fuel = 1000; % Kg/m^3
Ve = 3000; % m/s
Mbody = 60684; % Kg
Mfuel = 488370; % Kg
% Mfuel = Mfuel + 450000;
Mpayload = 31184; % Kg
t0 = 0;
Tf = 1000;
m0 = Mbody + Mfuel + Mpayload;
a0 = 10;
% constant acceleration input
a_in = a0*ones(1, Tf+1);
tspan = [t0 Tf];
H = [-1;0;1;2;3;4;5;6;7;8;9;10;11;13;15;17;20;25;30;32;35;40;45;47;50;51;60;70;71;80;84.9;89.7;100.4;105;110;120];
H = H*1000;
rho = [1.347;1.225;1.1116;1.0065;0.9091;0.8191;0.7361;0.6597;0.5895;0.5252;0.4664;0.4127;0.3639;0.2655;0.1937;0.1423;0.0880;0.0395;0.0180;0.0132;0.0082;0.0039;0.0019;0.0014;0.0010;0.00086;0.000288;0.000074;0.000064;0.000015;0.000007;0.000003;0.0000005;0.0000002;0.0000001;0.0000001];
f_rho = fit(H, rho, 'linearinterp');
[t, y] = ode45(@(t, m) myode(t, m, a_in, f_rho), tspan, [m0, a0]);
figure(1)
plot(t, y(:,1)), grid on, xlabel('t'), title('Rocket Mass (kg)')
figure(2)
plot(t, y(:,2)/1e3), grid on, xlabel('t'), title('Vertical Velocity (km/s)')
figure(3)
plot(t, y(:,2).*t/1e3), grid on, xlabel('t'), title('Rocket Altitude (km)')
function dydt = myode(t, y, a_in, f_rho)
global Ve Cd A_r g Tf Mbody Mpayload
m = y(1) - Mbody - Mpayload; % Final mass = Mbody + Mpayload
u = y(2);
dudt = interp1(0:Tf, a_in, t);
h = u * t;
rho_air = f_rho(h);
Fd = 0.5 * rho_air * u^2 * A_r * Cd;
dmdt = (1/Ve) * (- m * g - Fd - m * dudt);
dydt = [dmdt; dudt];
end

Connectez-vous pour commenter.


Jan
Jan le 8 Déc 2022
[t,m] = ode45(@(t,m) myode(t,m), tspan, m0);
Here myode has 2 inputs only.
function dydt = myode(t,y,a_in)
Here it has 2 inputs. I assume you want to provide a_in also:
[t,m] = ode45(@(t,m) myode(t, m, a_in), tspan, m0);
Your equation has two variables. y(1) is related to m and y(2) to u, as far as I can see. But you do not integrate u with ODE45? Then it is not a part of the equation to be integrated.
Note, that ODE45 is designed to integrate smooth functions only. Linear interpolations are not smooth. If you are lucky the integration stops with an error message. In many cases ODE45 reduces the step size until the discontinuity is hidden by rounding errors, but the accumulated errors due to the huge number of steps can cause very random final values. Do not call this "result".

Tags

Produits


Version

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by