Simulation eines Raketenstarts

Simulation eines Raketenstarts#

Sie sind Wissenschaftler:in im Head-Bureau for Rocket Simulation (H-BRS) und wurden beauftragt einen geplanten Raketenstart zu simulieren. Von der Raketenplattform in Sankt Augustin sollen vier Satelliten mit einem Gesamtgewicht von 2952 kg in MEO (Medium Earth Orbit) auf eine leicht elliptische Laufbahn in 22922 km Höhe gebracht werden, also etwas unterhalb einer geostationären Umlaufbahn.

Zu Beginn der Planungsphase sollen eine Reihe Simulationen durchgeführt werden, um den Flugplan zu optimieren. Für Sie eigentlich eine Routineuntersuchung, aber da Ihre Mitarbeiter:innen noch nie etwas von Versionskontrolle gehört haben, sind leider Teile Ihres Simulationsprogrammes abhandengekommen.

Hinweis

Das Ziel dieser Übung ist, Ihnen 1.) das Aufstellen und Programmieren mehrdimensionaler Differentialgleichungen und 2.) das Arbeiten mit größtenteils bereits fertigem Code näher zu bringen. Sie müssen nur rocket_ode.m selber schreiben, das Meiste ist bereits für Sie programmiert!

Sie finden unten die sechs .m-Dateien mit jeweils einer kurzen Beschreibung. Sie können sich den Inhalt der Dateien durch Aufklappen der Codeboxen anschauen.

Ihre Aufgaben finden Sie weiter unten unterhalb der Beschreibung der Simulationssoftware.

Ihre Simulationssoftware#

Als Matlab-Expert:in haben Sie eine einfache Simulationssoftware für die frühe Planungsphase entworfen. Ihr Programm simuliert ein “gravity turn maneuver” in 2D und besteht aus den folgenden sechs Dateien.

`rocket_launch.m`#

Dies ist das Hauptskript Ihrer Simulation. Hier werden die Simulationsparameter festgelegt, die zugrunde liegende gewöhnliche Differentialgleichung gelöst und, je nach Bedarf, eine Animation der Simulation gestartet.

Hinweis

Die Simulation beinhaltet eine Animation und kann deswegen nicht interaktiv auf dieser Seite laufen. Laden Sie sich die Dateien herunter und rechnen Sie lokal auf Ihrem Computer.

Show code cell source Hide code cell source

%%file rocket_launch.m
%% 2D rocket gravity turn simulation

% this script simulates a gravity turn rocket launch in 2D cartesian
% coordinates (3D earth is just for visualization)

% The physical parameters of the rocket are very unlikely to be realistic...
% Please do not build a space rocket based on this code

% dependencies:
%   rocket_ode.m
%   mass_flow_rate.m
%   F_thrust.m
%   rocket_status.m
%   rocket_visualization.m

clear;
close all;

%% define parameters inside global struct

global params;

params.G = 6.67408e-11;        % gravitational constant [m3 kg-1 s-2]

params.earth.r  = 6378137;     % radius of earth at the equator [m]
params.earth.m  = 5.972e24;    % mass of the earth [kg]

% inspired by Ariane flight VA244
params.rocket.m0  = 84852;          % rocket dry mass (inc. payload, without fuel) [kg]
params.rocket.mf  = 660000;         % rocket fuel mass [kg]
params.rocket.mfr = 968;            % rocket ideal mass flow rate [kg/s]
params.rocket.ve  = 8135;           % exit velocity at nozzle [m/s]
params.rocket.target_orbit =  22.922e6; % target orbit [m]

% some parameters for the simulation
params.sim.tspan             = [0 24*3600];   % simulation time [s]
params.sim.plot_during_solve = false;        % set to true to update graphics during 
                                             % numerical integration
params.sim.saveScreenShots   = false;        % set to true to save screenshots
                                             % in subdir "img"

%% initial conditions (position, velocity and mass)
% save rocket status (position, velocity and mass) in new variable q
% q = [x0;x1;v0;v1,m] = [x,v], x = [x0;x1], v=[v0;v1]

% inital flight phase is called 'lift off'
params.rocket.status = 'lift off';

% initial position and velocity
xinit = [params.earth.r; 0];
vinit = [0;0];

% initial mass is dry mass + mass of propellant
minit = params.rocket.m0 + params.rocket.mf;

qinit = [xinit; vinit; minit];

%% solve the ODE

opts = odeset('AbsTol',1e-8,'RelTol',1e-8,...
              'OutputFcn',@rocket_visualization,...
              'Events',@rocket_status);
sol = ode45(@rocket_ode, params.sim.tspan, qinit, opts);

%% plot the animation

if ~params.sim.plot_during_solve
    t = linspace(sol.x(1), sol.x(end), 120);
    q = deval(t,sol);
    
    rocket_visualization(params.sim.tspan(1),qinit, 'init', false, sol.xe);
    for i = 1:length(t)
       ti = t(i);
       qi = q(:,i)';
       rocket_visualization(ti,qi,[],false, sol.xe);
    end
end

run rocket_launch.m

`rocket_ode.m`#

Das Herzstück Ihrer Simulation ist das Lösen eines Systems gewöhnlicher Differentialgleichungen erster Ordnung \(\dot{\mathbf{q}} = \) rocket_ode\((t,\mathbf{q})\) mit dem Löser ode45. Die rechte Seite ist in dieser Datei implementiert. Dabei ist \(\mathbf{q} = [x,y, v_x, v_y, m]^T\) ein Spaltenvektor, der die aktuelle Position \((x,y)\), die aktuelle Geschwindigkeit \((v_x,v_y)\) sowie die aktuelle Masse \(m\) der Rakete beschreibt. Leider fehlt ausgerechnet diese Datei.

`mass_flow_rate.m`#

Teil der Raketen-ODE rocket_ode.m ist die gewöhnliche Differentialgleichung für die Masse \(m\) der Rakete \(\dot{m}=g(m, \tau)\), wobei \(\tau\) die Antriebsdrossel modelliert (\(\tau=0\) bedeutet kein Antrieb, \(\tau=1\) bedeutet maximaler Antrieb). Die Funktion \(g\) ist in mass_flow_rate.m implementiert. Ein Raketenantrieb funktioniert nach dem dritten Newton’schen Prinzip über den schnellen Ausstoß von Treibstoff. Je schneller der Treibstoff ausgestoßen wird, desto höher sind der momentane Impuls und die Gewichtsabnahme der Rakete.

`F_thrust.m`#

In dieser Datei ist die Funktion \(F_{\text{thrust}}(t,\mathbf{q})\) hinterlegt, die abhängig des aktuellen Zustandes den Antrieb der Rakete in \(x\)- und \(y\)-Richtung berechnet. Diese Datei haben Sie bereits von Ihren Kolleg:innen von der Flugplanung erhalten. Der Raketenstart ist aufgeteilt in die Phasen lift off, pitch over, ascend, 2nd stage, und orbit. Je nach Phase wird in F_thrust der Antrieb gedrosselt oder der Angriffswinkel der Rakete geändert. Diese Funktion bestimmt zusammen mit rocket_status.m den Flugplan der Rakete. Optional kann die berechnete Drosselung throttle als zweiter Parameter mit ausgegeben werden. Sie wird unter anderem gebraucht um den Massefluss \(g\) (bzw. mass_flow_rate) zu berechnen.

`rocket_status.m`#

Diese Datei wird von F_thrust.m verwendet. Hier wird je nach Zustand der Rakete nach bestimmten Kriterien die Phase (lift off, pitch over, ascend, 2nd stage, orbit) des Raketenstarts gewechselt. In jeder Phase verhält sich die Funktion F_thrust etwas anders. Die Syntax ist so gewählt, dass die Funktion vom event handler des ODE-Lösers von Matlab aufgerufen wird.

`rocket_visualization.m`#

Diese Datei implementiert die Animation sowie Plots des Raketenstarts. Die Syntax der Funktion ist so gewählt, dass die Funktion während der numerischen Integration in Matlabs ODE Löser in jedem Zeitschritt einen visuellen Output liefern kann.

Show code cell source Hide code cell source

%%file rocket_visualization.m
%% callback function to visualize the rocket launch
%  during numerical integration
function status = rocket_visualization(t,q,flag,simple, event_times)
% This is an OutputFcn for Matlab's ODE solvers. The additional parameter
% simple can be set to true to only print some info to the command window.
% event_times can be passed to determine the transition of flight maneuvers
% from the ode45 output structure


% params contains some simulation as well as physical and rocket parameters 
% and handles is a struct containing handles to the plots
global handles params;

if nargin<4
    if ~params.sim.plot_during_solve
        simple=true;
    else
        simple=false;
    end
end

if params.sim.plot_during_solve && params.sim.saveScreenShots
    warning('You are writing screenshots during an online visualization. This is likely to be slow and create a butt-load of images.')
end

% This happens sometimes. Don't know why
if numel(t)<=0
    status = 0;
    return;
end

if simple
    disp(['t = ' num2str(t(1))])
    status = 0;
    return
end

if ~params.sim.plot_during_solve
    if nargin<5
        rocketstatus = 'unknown';
    else
        % determine rocket status
        status = {'lift off', 'pitch over', 'ascend', '2nd stage', 'orbit'};
        i=1;
        while i <= length(event_times) && t(1) > event_times(i)
            i=i+1;
        end
        rocketstatus = status{i};
    end
else
    rocketstatus = params.rocket.status;
end

if strcmp(flag, 'init')
    
    % initialize plots using initial conditions
    figure('units','normalized','outerposition',[0 0.1 1 1], 'color', 'w')
    
    %% animation
    subplot(2,2,1);
    if t > 0
        s = '+';
    else
        s = '-';
    end
    h = fix(t(1)/3600);
    min = fix(t(1)/60 - h*60);
    sec = abs(t(1) - h*3600 - min*60);
    title(['t = ',s, num2str(h,'%02i'), ':',...
                     num2str(min,'%02i'), ':',...
                     num2str(sec,'%02.3f')])
    
    % plot earth
    plot_earth(t(1),flag);
    set(gca, 'Xtick',[], 'YTick', [], 'ZTick', [], 'box', 'off');
    hold on
    
    % plot target orbit
    r = params.rocket.target_orbit + params.earth.r;
    handles.orbit = rectangle('position',[-r, -r, 2*r, 2*r],...
                              'curvature',[1,1],...
                              'LineStyle','-.',...
                              'EdgeColor',[0 0.4470 0.7410],...
                              'LineWidth',1);

    % current rocket orientation (this does not work well with postprocessed animation)
    handles.currentori = plot([q(1), q(1) + 5e5],...
                              [q(2), q(2)],...
                              'LineWidth',4,...
                              'Color','k');

    % trajectory of rocket
    handles.trajectory = plot(q(1), q(2),'LineWidth',2);
    set(gca,'Fontsize',16)
    xlabel(['status: ',rocketstatus]);
    daspect([1 1 1])
    
    %% mass of fueL
    subplot(2,2,2);
    fuel = (q(5) - params.rocket.m0);
    handles.fuel = plot(t(1),fuel/1000,'LineWidth',2);
    title('mass of propellant')
    xlabel('t [sec]');
    ylabel('propellant mass [t]');
    set(gca,'Fontsize',16)
    
    %% height above sea level
    subplot(2,2,3)
    height = norm(q(1:2)) - params.earth.r;
    handles.orbit  = plot(t(1),params.rocket.target_orbit/1000,...
                          '-.',...
                          'DisplayName','target orbit',...
                          'LineWidth',2);
    hold on
    handles.height = plot(t(1),height/1000,...
                          'DisplayName','trajectory',...
                          'LineWidth',2);
    xlabel('t [sec]');
    ylabel('height [km]');
    title('height above sea level')
    set(gca,'Fontsize',16)
    legend('location','southeast')
    
    %% lateral velocity
    subplot(2,2,4)
    R = params.rocket.target_orbit + params.earth.r;
    veltarget = sqrt(params.G*params.earth.m/R)*3.6;
    %vel = q(3:4)'*[-q(2);q(1)]/norm(q(1:2))*3.6;
    vel = norm(q(3:4));
    handles.targetvel  = plot(t(1),veltarget,...
                              '-.',...
                              'DisplayName','target velocity',...
                              'LineWidth',2);
    hold on
    handles.vel = plot(t(1),vel,...
                       'DisplayName','flight velocity',...
                       'LineWidth',2);
    xlabel('t [sec]');
    ylabel('velocity [km/h]');
    title('lateral velocity')
    set(gca,'Fontsize',16)
    legend('location','southeast')
    
    for i=2:4
        subplot(2,2,i);
        ax = gca;
        outerpos = ax.OuterPosition;
        ti = ax.TightInset; 
        left = outerpos(1) + ti(1);
        bottom = outerpos(2) + ti(2);
        ax_width = outerpos(3) - ti(1) - ti(3);
        ax_height = outerpos(4) - ti(2) - ti(4);
        ax.Position = [left bottom ax_width ax_height]; 
    end
    
    
elseif strcmp(flag, 'done')
    % pass
else
    % do this in every time step of tspan
    % first subplot is animation of trajectory. Title is Simulation time.
    
    % rotate earth and update time
    subplot(2,2,1);
    if t > 0
        s = '+';
    else
        s = '-';
    end
    h = fix(t(1)/3600);
    min = fix(t(1)/60 - h*60);
    sec = abs(t(1) - h*3600 - min*60);
    title(['t = ',s, num2str(h,'%02i'), ':',...
                     num2str(min,'%02i'), ':',...
                     num2str(sec,'%02.3f')])
    xlabel(['status: ',rocketstatus]);
    plot_earth(t(1),flag);

    % update current rocket orientation
    set(handles.currentori, 'Xdata', [q(1), q(1) + params.rocket.dir(1)*5e5]);
    set(handles.currentori, 'Ydata', [q(2), q(2) + params.rocket.dir(2)*5e5]);
    
    % update trajectory of rocket
    X= get(handles.trajectory,'Xdata');
    Y= get(handles.trajectory,'Ydata');
    set(handles.trajectory,'Xdata',[X q(1)]);
    set(handles.trajectory,'Ydata',[Y q(2)]);
    
    % update fuel plot
    fuel = (q(5) - params.rocket.m0);
    if ~strcmp(rocketstatus,'orbit')
        T = get(handles.fuel,'Xdata');
        f = get(handles.fuel,'Ydata');
        set(handles.fuel,'Xdata',[T t(1)]);
        set(handles.fuel,'Ydata',[f fuel/1000]);
    end
    
    % update height above sea level plot
    height = norm(q(1:2)) - params.earth.r;
    T = get(handles.height,'Xdata');
    H = get(handles.height,'Ydata');
    set(handles.height,'Xdata',[T t(1)]);
    set(handles.height,'Ydata',[H height/1000]);
    H = get(handles.orbit,'Ydata');
    set(handles.orbit,'Xdata',[T t(1)]);
    set(handles.orbit,'Ydata',[H params.rocket.target_orbit/1000]);
    
    % update velocity plot
    R = params.rocket.target_orbit + params.earth.r;
    veltarget = sqrt(params.G*params.earth.m/R);
    %vel = q(3:4)*[-q(2);q(1)]/norm(q(1:2));
    vel = norm(q(3:4));
    T = get(handles.targetvel,'Xdata');
    H = get(handles.targetvel,'Ydata');
    set(handles.targetvel,'Xdata',[T t(1)]);
    set(handles.targetvel,'Ydata',[H veltarget*3.6]);
    H = get(handles.vel,'Ydata');
    set(handles.vel,'Xdata',[T t(1)]);
    set(handles.vel,'Ydata',[H vel*3.6]);
    
    drawnow;
    
    if params.sim.saveScreenShots
        step=size(dir(['img/*.png']),1);
        filename = ['img/rocket_launch_',num2str(step, '%03i'),'.png'];
        saveas(gcf, filename);
    end
end
    status = 0;
end

function plot_earth(t, flag)
% modified from Ryan Gray
% https://de.mathworks.com/matlabcentral/fileexchange/13823-3d-earth-example

global params handles;

if strcmp(flag, 'init')
    % initialize the earth surface plot
    
    npanels = 180;   % Number of globe panels around the equator deg/panel = 360/npanels
    alpha   = 1; % globe transparency level, 1 = opaque, through 0 = invisible
    [x, y, z] = ellipsoid(0, 0, 0, params.earth.r, params.earth.r, params.earth.r, npanels);
    handles.globe = surf(x, y, -z, 'FaceColor', 'none', 'EdgeColor', 0.5*[1 1 1]);

    %% Texturemap the globe

    % Load Earth image for texture map
    % https://commons.wikimedia.org/wiki/File:Land_ocean_ice_2048.jpg
    % NASA Goddard Space Flight Center Image by Reto Stöckli (land surface, 
    % shallow water, clouds). Enhancements by Robert Simmon (ocean color, 
    % compositing, 3D globes, animation). Data and technical support: MODIS 
    % Land Group; MODIS Science Data Support Team; MODIS Atmosphere Group; 
    % MODIS Ocean Group Additional data: USGS EROS Data Center (topography); 
    % USGS Terrestrial Remote Sensing Flagstaff Field Center (Antarctica); 
    % Defense Meteorological Satellite Program (city lights). 
    % [CC BY-SA 3.0 (http://creativecommons.org/licenses/by-sa/3.0/)]
    
    worldmap = 'Land_ocean_ice_2048.jpg';
    if ~exist(worldmap, 'file')==1
        websave(worldmap,...
            'http://upload.wikimedia.org/wikipedia/commons/thumb/c/cd/Land_ocean_ice_2048.jpg/1024px-Land_ocean_ice_2048.jpg');
    end
    cdata = imread(worldmap);

    % Set image as color data (cdata) property, and set face color to indicate
    % a texturemap, which Matlab expects to be in cdata. Turn off the mesh edges.

    set(handles.globe, 'FaceColor', 'texturemap', 'CData', cdata, 'FaceAlpha', alpha, 'EdgeColor', 'none');
    
    view([0.5 0.25 0.5])
end

     %rotate earth once around polar axis every 24 hours
    RE = makehgtform('zrotate',(2*pi/(24*3600))*t);

    % rotate earth so that launch position matches H-BRS coordinates
    RZ = makehgtform('zrotate', -7.18*pi/180);
    RY = makehgtform('yrotate', 50.78*pi/180);

    hgx = hgtransform;
    set(hgx,'Matrix', RY*RZ*RE);
    set(handles.globe,'Parent',hgx);

end

Aufgabe 1#

Damit Ihre Simulationssoftware lauffähig ist, füllen Sie die Datei rocket_ode.m mit neuem Inhalt. Die Differentialgleichungen, die gelöst werden müssen, lauten

\[\begin{split} \begin{align} \ddot{\mathbf{z}} &= \frac{\mathbf{F}_{\text{G}} + \mathbf{F}_{\text{thrust}}}{m}, \notag \\ \dot{m} &= g(m, \tau). \notag \end{align} \end{split}\]

Die erste der beiden DGLs beschreibt die Beschleunigung der Rakete, wobei \(\mathbf{z} = [x, y]^T\) die aktuelle Position in karthesischen Koordinaten ist. Die Beschleunigung setzt sich zusammen aus der Gravitationsbeschleunigung

\[ \frac{\mathbf{F}_{\text{G}}}{m} = -\frac{G\cdot M}{\|\mathbf{z}\|^3} \mathbf{z} \]

und der Beschleunigung durch den Raketenantrieb \(\mathbf{F}_{\text{thrust}}/m\). Hierbei ist \(G=6.67408 \cdot 10^{-11} \frac{\textrm{m}^3}{\textrm{kg}\cdot\textrm{s}^2}\) die Gravitationskonstante und \(M=5.972 \cdot 10^{24} \textrm{kg}\) die Erdmasse.

Die zweite Differentialgleichung beschreibt die Masseänderung der Rakete. Der Massefluss hängt von der aktuellen Masse \(m\) der Rakete und der Drosselung \(\tau\) ab. Die Funktionen \(\mathbf{F}_{\text{thrust}}(t,\mathbf{q})\) sowie \(g(m,\tau)\) sind bereits in den Matlab-Dateien F_thrust.m und mass_flow_rate.m implementiert. Die Drosselung (throttle) \(\tau\) wird von F_thrust.m berechnet und als zweiten Ausgabeparameter bereitgestellt:

[Ft, throttle] = F_thrust(t,q)

Tipp

Überführen Sie zunächst dieses System von Differentialgleichungen in ein System erster Ordnung.

\[\begin{split} \begin{align} \dot{\mathbf{q}} &= \textrm{rocket_ode}(t, \mathbf{q})\\ \text{mit }\mathbf{q} &= \begin{pmatrix} x\\ y\\ v_x\\ v_y\\ m \end{pmatrix}. \end{align}\end{split}\]

Aufgabe 2#

Wie schätzen Sie die Modellierungstiefe des Simulationsprogrammes ein? Nennen Sie die wichtigsten Modellvereinfachungen und diskutieren Sie, welche am kritischsten sind. Welche Modellvereinfachungen würden Sie als erstes angehen, wenn Sie das Modell verfeinern würden? Wie schätzen Sie den Aufwand dafür ein?

Zusatzaufgabe#

Bestimmen Sie auf Grundlage Ihrer Simulation welchen Beschleunigungen die Satelliten während des Raketenstarts ausgesetzt sind. Laut Satellitenhersteller sind die kritischen Komponenten für Beschleunigungen bis 10 G ausgelegt. Erfüllt der von Ihnen simulierte Raketenstart diese Anforderungen? In welcher Phase des Raketenstarts sind die Beschleunigungen am höchsten? Erstellen Sie einen Plot, der den Betrag der Beschleunigung in G über die Zeit während des Starts darstellt.

Zusatzinformationen

Der simulierte Raketenstart wurde inspiriert vom Ariane flight VA244 am 25. Juli 2018.

http://www.arianespace.com/mission/ariane-flight-va244/
https://en.wikipedia.org/wiki/Ariane_flight_VA244
https://en.wikipedia.org/wiki/Ariane_5#Variants
Einige der Galileo-Satelliten sind auch in der Datenbank des Space Debris Viewers hinterlegt (allerdings noch nicht die Satelliten 23-24-25-26 von flight VA244).

- Hinweis: Dieses Video ist von YouTube aus eingebunden und nicht Teil des frei lizenzierten Materials! -

Simulation eines Raketenstarts

Contents

Simulation eines Raketenstarts#

Ihre Simulationssoftware#

rocket_launch.m#

rocket_ode.m#

mass_flow_rate.m#

F_thrust.m#

rocket_status.m#

rocket_visualization.m#