Solving Equations of motion with ode15s

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?

[Torsten]

We can't tell because we don't have executable code to run.

[Brendan Görres]

%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

[Torsten]

dydt;dvdt;dhdt;dxddt;dmdt;dwdt are all undefined.

[Torsten]

Consider,e.g., the units in

dydt=(thrust*sin(alpha)/m*v)-(v*cos(y)/h)-gh*cos(y);

All three contributions to dydt have a different unit. There must be an error here.

(I did not check the other equations).

[Brendan Görres]

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

[Brendan Görres]

Yeah you are right. Had a few typos in my code but changed it to

dydt=(thrust*sin(alpha)/m*v)-(v*cos(y)/h)-(gh*cos(y))/v;
dvdt=(thrust*cos(alpha)/m)-gh*sin(y)-(rhoh*Adrag*(v^2)/2*m)*cos(alpha);

But that doesn't change my original problem. How is it possible that the function of alpha has no influence on the outcome of my EoM?

[Torsten]

Must still be wrong.

thrust*sin(alpha)/m*v has unit m^2/s^3, the other contributions have unit 1/s.

thrust*cos(alpha)/m has unit m/s^2, (rhoh*Adrag*(v^2)/2*m)*cos(alpha) has unit m/s^2 * kg^2.

[Brendan Görres]

no they all have rad/s which is correct because dydt is an angle velocity.

dydt[unit]=(kg*m*rad*s/kg*m)-(m*rad/s*m)-(m*rad*s/s^2*m)

which is all rad/s.

[Torsten]

unit(thrust*sin(alpha)/m*v) = kg*m/s^2 * rad / kg * m/s = rad * m^2 / s^3

I think you mean

unit(thrust*sin(alpha)/(m*v)) = kg*m/s^2 * rad / (kg * m/s) = rad/s

Same for the other expression.

[Brendan Görres]

ok I see I forgot the brackets, but again can you help me with my question?

[Torsten]

And the values you get for the variables look realistic to you ? E.g. a velocity of 12000 m/s ?

I suggest you reinspect your equations first before going into details about the setting of "alpha".