Solving Equations of motion with ode15s

Question

Brendan Görres 2019-5-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/460024-solving-equations-of-motion-with-ode15s

评论： Torsten 2019-5-6

Recently I am trying to solve a system of equations with ode15s. The goal is to calculate the Ascent trajectory of a Rocket with a control variable alpha, which is the angle between the velocity vector v and the thrust vector, so that by changing alpha manually I should change the flight direction.

%Thrust vector control variable
alpha=(-2*pi/180)
%alpha=(-2*pi/180)*t+10
%alpha=0.5*(-2*pi/180)*t^2+10*t+0.01...assumed to be konstant/linear/quadratic[rad]
%functions
gh=g0*((re/h)^2);%Gravi*tationsbeschleunigung abh. von h
rhoh=rho0*exp(-h/sh);%Luftdichte abh. von h
%EoM
dydt=(thrust*sin(alpha)/m*v)-(v*cos(y)/h)-gh*cos(y);
dvdt=(thrust*cos(alpha)/m)-gh*sin(y)-(rhoh*Adrag*(v^2)/2*m);
dhdt=v*sin(y);
dxddt=(re/h)*v*cos(y);
dmdt=-propflow;
dwdt=(1/Ix)*thrust*sin(alpha)*rs;

in a first run I assume alpha to be konstant and I experimented a little bit with the value (so for example alpha=-2*... and in the next run alpha= +10*...). In a next step I assumed alpha to be linear so for example

alpha=(-2*pi/180)*t+10

and in a last step I assumed it to be quadratic:

alpha=0.5*(-2*pi/180)*t^2+10*t+0.01

I am using ode15s to solve the equations, because my code needs to include the algebraic variable t.

I plotted the results over t and I receive the exact same plots for every single run, regardless of alpha being konstant linear or quadratic.

How is that possible?

11 个评论
显示 9更早的评论隐藏 9更早的评论

Brendan Görres 2019-5-3

在 MATLAB Online 中打开

%earth parameterssni
re=6378e3;%[m]
g0=9.81;%Erdbeschleunigung[m/s^2]
%vehicleparameters
m0=68e3;%Startmasse[kg]
Isp=390;% spezifischer Impuls[s]
ceff=Isp*g0;%effektive Austrittsgeschwindigkeit[ms/]
thrust=933910;%Schub[N]
propflow=thrust/ceff;%Massenstrom Treibstoff[kg/s]
R=15;%1/MR
mp=m0*(1-(1/R));%Treibstoffmasse[kg]
mf=m0-mp;%finale Masse[kg]
tburn=mp/propflow;%Brennzeit[s]
%pitchover values
hpo=re+130;%pitchover altitude
ghpo=g0*((re/(hpo))^2);%Erdbeschleunigung bei hpo[m/s^2]
syms x
f=(-ghpo/2)*(((m0-x)/propflow)^2)+(ceff/propflow)*(x*log(x/m0)+m0-x)==130;
mpo=double(vpasolve(f,x));
vpo=-ghpo*((m0-mpo)/propflow)-ceff*log(mpo/m0);
tpo=((m0-mpo)/propflow);
%Integration time
tmax=tburn-tpo;
tRange=[0 tmax];
%Initial values
y0=89.85*pi/180;%[rad]
v0=vpo;%Anfangsgeschwindigkeit[m/s]
h0=hpo;%Anfangshöhe[m]
x0=0;%Downrange zu Beginn[m]
w0=0;%[rad]
Y0=[y0;v0;h0;x0;m0;w0];
%Solve System of Equations
options = odeset('RelTol',1e-3,'AbsTol',1e-6);
[tSol,YSol]=ode15s(@APalpha,tRange,Y0,options);
ySol=YSol(:,1)*180/pi;%[deg]
vSol=YSol(:,2);
hSol=YSol(:,3)/1000;%[km]
xdSol=YSol(:,4)/1000;%[km]
mSol=YSol(:,5);
wSol=YSol(:,6)*180/pi;%[deg]
%plot functions
plot(tSol,ySol)
xlim([0 tmax]);
xlabel('time[s]');
ylabel('flight path angle[deg]');
figure
plot(tSol,vSol)
xlim([0 tmax]);
xlabel('time[s]');
ylabel('velocity[m/s]');
figure
plot(tSol,hSol)
xlim([0 tmax]);
xlabel('time[s]');
ylabel('height[km]');
figure
plot(tSol,xdSol)
xlim([0 tmax]);
xlabel('time[s]');
ylabel('downrange[km]');
figure
plot(tSol,mSol)
xlim([0 tmax]);
xlabel('time[s]');
ylabel('vehicle mass[kg]');
figure
plot(xdSol,hSol)
xlabel('downrange[km]');
ylabel('height[km]');
plot(tSol,wSol)
xlabel('time[s]');
ylabel('Rotation angle velocity');

and the function

function dYdt= APalpha(t,Y)
format short
y=Y(1);
v=Y(2);
h=Y(3);
xd=Y(4);
m=Y(5);
w=Y(6);
%earth parameters
re=6378e3;%Erdradius[m]
g0=9.81;%Erdbeschleunigung[m/s^2]
rho0=1.225;%Luftdichte[kg/m^3] 
sh=7500;%skalare Höhe[m]
%vehicle parameters
m0=68e3;%Startmasse[kg]
thrust=933910;%Schub[N]
diameter=5;%Durchmesser der Rakete[m]
R=diameter/2;
Cdrag=0.5; %Strömungswiderstandskoeffizient
Isp=390;%spezifischer Impuls[s]
H=20;%Height of the Rocket[m]
area=pi*(R^2);%Raketenfläche[m^2]
Adrag=Cdrag*area;%angeströmte Fläche[m^2]
ceff=Isp*g0;%[m/s]
propflow=thrust/ceff;% Treibstoffmassenstrom[kg/s]
V0=pi*R^2*H;%initial Volume of the rocket[m^3]
rhop=m0/V0;%Dichte Treibstoff[kg/m^3]=const.
rs=10;%Hebelarm zum Schwerpunkt[m]
tburn=m0/propflow;%Brenndauer
rt=(2.5/tburn)*t;%Innenradius der Treibstoffsäule[m]
Ix=(pi*rhop*H/12)*(3*(R^4-rt^4)+H^2*(R^2-rt^2));%Trägheitsmoment um x-Achse[kg*m^2]
dYdt=[dydt;dvdt;dhdt;dxddt;dmdt;dwdt];
end

Brendan Görres 2019-5-3

在 MATLAB Online 中打开

oh, I forgot to implement them into my question. This is the whole code I wrote:

%earth parameterssni
re=6378e3;%[m]
g0=9.81;%Erdbeschleunigung[m/s^2]
%vehicleparameters
m0=68e3;%Startmasse[kg]
Isp=390;% spezifischer Impuls[s]
ceff=Isp*g0;%effektive Austrittsgeschwindigkeit[ms/]
thrust=933910;%Schub[N]
propflow=thrust/ceff;%Massenstrom Treibstoff[kg/s]
R=15;%1/MR
mp=m0*(1-(1/R));%Treibstoffmasse[kg]
mf=m0-mp;%finale Masse[kg]
tburn=mp/propflow;%Brennzeit[s]
%pitchover values
hpo=re+130;%pitchover altitude
ghpo=g0*((re/(hpo))^2);%Erdbeschleunigung bei hpo[m/s^2]
syms x
f=(-ghpo/2)*(((m0-x)/propflow)^2)+(ceff/propflow)*(x*log(x/m0)+m0-x)==130;
mpo=double(vpasolve(f,x));
vpo=-ghpo*((m0-mpo)/propflow)-ceff*log(mpo/m0);
tpo=((m0-mpo)/propflow);
%Integration time
tmax=tburn-tpo;
tRange=[0 tmax];
%Initial values
y0=89.85*pi/180;%[rad]
v0=vpo;%Anfangsgeschwindigkeit[m/s]
h0=hpo;%Anfangshöhe[m]
x0=0;%Downrange zu Beginn[m]
w0=0;%[rad]
Y0=[y0;v0;h0;x0;m0;w0];
%Solve System of Equations
options = odeset('RelTol',1e-3,'AbsTol',1e-6);
[tSol,YSol]=ode15s(@APalpha,tRange,Y0,options);
ySol=YSol(:,1)*180/pi;%[deg]
vSol=YSol(:,2);
hSol=YSol(:,3)/1000;%[km]
xdSol=YSol(:,4)/1000;%[km]
mSol=YSol(:,5);
wSol=YSol(:,6)*180/pi;%[deg]
%plot functions
plot(tSol,ySol)
xlim([0 tmax]);
xlabel('time[s]');
ylabel('flight path angle[deg]');
figure
plot(tSol,vSol)
xlim([0 tmax]);
xlabel('time[s]');
ylabel('velocity[m/s]');
figure
plot(tSol,hSol)
xlim([0 tmax]);
xlabel('time[s]');
ylabel('height[km]');
figure
plot(tSol,xdSol)
xlim([0 tmax]);
xlabel('time[s]');
ylabel('downrange[km]');
figure
plot(tSol,mSol)
xlim([0 tmax]);
xlabel('time[s]');
ylabel('vehicle mass[kg]');
figure
plot(xdSol,hSol)
xlabel('downrange[km]');
ylabel('height[km]');
plot(tSol,wSol)
xlabel('time[s]');
ylabel('Rotation angle velocity');
function dYdt= APalpha(t,Y)
format short
y=Y(1);
v=Y(2);
h=Y(3);
xd=Y(4);
m=Y(5);
w=Y(6);
%earth parameters
re=6378e3;%Erdradius[m]
g0=9.81;%Erdbeschleunigung[m/s^2]
rho0=1.225;%Luftdichte[kg/m^3] 
sh=7500;%skalare Höhe[m]
%vehicle parameters
m0=68e3;%Startmasse[kg]
thrust=933910;%Schub[N]
diameter=5;%Durchmesser der Rakete[m]
R=diameter/2;
Cdrag=0.5; %Strömungswiderstandskoeffizient
Isp=390;%spezifischer Impuls[s]
H=20;%Height of the Rocket[m]
area=pi*(R^2);%Raketenfläche[m^2]
Adrag=Cdrag*area;%angeströmte Fläche[m^2]
ceff=Isp*g0;%[m/s]
propflow=thrust/ceff;% Treibstoffmassenstrom[kg/s]
V0=pi*R^2*H;%initial Volume of the rocket[m^3]
rhop=m0/V0;%Dichte Treibstoff[kg/m^3]=const.
rs=10;%Hebelarm zum Schwerpunkt[m]
tburn=m0/propflow;%Brenndauer
rt=(2.5/tburn)*t;%Innenradius der Treibstoffsäule[m]
Ix=(pi*rhop*H/12)*(3*(R^4-rt^4)+H^2*(R^2-rt^2));%Trägheitsmoment um x-Achse[kg*m^2]
%Thrust vector control variable
alpha=-2*pi/180;%[rad]
%%alpha=(-2*pi/180)*t+10
%alpha=0.5*(-2*pi/180)*t^2+10*t+0.01....assumed to be konstant/linear/quadratic[rad]
%functions
gh=g0*((re/h)^2);%Gravitationsbeschleunigung abh. von h
rhoh=rho0*exp(-h/sh);%Luftdichte abh. von h
%EoM
dydt=(thrust*sin(alpha)/m*v)-(v*cos(y)/h)-gh*cos(y);
dvdt=(thrust*cos(alpha)/m)-gh*sin(y)-(rhoh*Adrag*(v^2)/2*m);
dhdt=v*sin(y);
dxddt=(re/h)*v*cos(y);
dmdt=-propflow;
dwdt=(1/Ix)*thrust*sin(alpha)*rs;
dYdt=[dydt;dvdt;dhdt;dxddt;dmdt;dwdt];
end