Kaltwasserrakete mit DGL berechnet

Flughöhe, Geschwindigkeit einer Wasserrakete durch Iteration der Differenzialgleichungen mit gnuplot berechnet.

DGL und Iterationen sind sicher nicht die Domäne von gnuplot, besser geeignet sind dafür sicher z.B. gnu octave, dort hat man entsprechende Funktionen auch mit besseren Verfahren.

Trotzdem, mit Gnuplot kann man leicht das einfachste Verfahren, das Polygonverfahren anwenden. Die Abweichung zur Realität ist am Anfang ganz akzeptabel, aber diese theoretische Gipfelhöhe erreicht die Rakete nur bei perfekter Bahnstabilisierung. Das ist nicht so einfach, in der Praxis kommt sie ins Taumeln, oder fliegt eine gekrümmte Bahn. Ab dann ist die Annahme des cw Wertes natürlich komplett daneben.

Skripte

Skript Download

#!/usr/bin/gnuplot
# Berechnung der Geschwindigkeit, der Flughöhe und ~Zeit einer Kaltwasserrakete bei senkrechtem Start
# Zahlenbeispiel mit einer Mehrweg-PET Fasche.
# (c) 2015, Joachim Schwender
set decimalsign locale
set decimal ","
dt = 0.001;       # Die Schrittweiter für die Berechnung [s]
#------- Raketeneigenschaften --------
D_d = 10.0/1000;        # Austrittsdurchmesser der Düse [m]
D_g = 70.0/1000;        # Durchmesser der Rakete [m]
m_geh = 0.049+0.15;     # Gehäusemasse [kg] + Ballast
vol_rak = 0.00155; # Gesamtvolumen der Rakete [m³]
cw = 0.60;         # Strömungswiderstandsbeiwert cw
#-------- Treibladung -----------
fuellgrad = 0.56;   # Füllgrad der Flasche, also Volumenanteil mit Wasser am Gesamtvolumen [%]
Druck = 16.0*100000.0;    # Druck der Gasfüllung beim Start der Rakete [Pa], Hinweis: 100000 Pa = 1 bar
Titel1 = gprintf("Kaltwasserrakete: Füllgrad %.2f", fuellgrad)
Titel2 = gprintf("Druck %.3s %cPa", Druck)
#
# ------ Konstanten ------
dichte_luft = 1.234;    # Dichte der Luft [kg/m³]
dichte_wasser = 1000.0 ;  # Dichte des Wassers [kg/m³]
g = 9.81;               # Erdbeschleunigung
k_luft = 1.4;           # isentropenexponent der Luft (trocken, 0°C), 
#                 zur Berechnung der adiabatischen Expansion des Treibgases
#------------------------
A_d = D_d**2.0 * pi/4.0;    # Fläche der Düse [m²]
A_g = D_g**2.0 * pi/4.0;    # Fläche des Raketengehäuses [m²]
vol_fl = fuellgrad * vol_rak;             # Füllvolumen der Flüssigkeit
vol_luft = (1.0-fuellgrad)*vol_rak;         # Gasvolumen in der Rakete
Energie = vol_rak * (1.0-fuellgrad)*Druck;  # Energie des Raketentreibsatzes der Rakete [J]
#      Startmasse der Rakete [kg], incl. Gasfüllung bei Anfangsdruck
masse = m_geh + vol_fl*dichte_wasser + vol_luft*dichte_luft*Druck/100000;  
#
#
set print "wasserrakete.dat"
load 'wasserrakete.func'
# Die Gipfelhöhe ermitteln
plot "wasserrakete.dat" using 1:2 with dots
set print
#    Da "set decimal locale" nicht auf Ausgaben wirkt miss dieser 
#    Kunstgriff her um den Dezimalseparator konsequent umzusetzen
system sprintf("mv %s /tmp/tttt; sed \"s/\\./,/g\" /tmp/tttt > %s", "wasserrakete.dat", "wasserrakete.dat")
hmax = GPVAL_DATA_Y_MAX
set print "hmax.dat"
print fuellgrad, hmax
set print
system sprintf("mv %s /tmp/tttt; sed \"s/\\./,/g\" /tmp/tttt > %s", "hmax.dat", "hmax.dat")
set title Titel1.Titel2
# Farbe der Kurve repräsentiert die Geschwindigkeit
set palette defined (0 'blue', 0.5 'red', 1 'yellow')
set style line 40 lc palette

# Papiergrösse leider in Inch! 
#  Der Plot muss etwas verkeinert und verschoben werden damit das nicht bis an den Rand geht
set size 0.95,0.88
set origin 0.03,0.06
set terminal pdf size 29.7cm,21.0cm
set output 'wasserrakete.pdf'
set ylabel gprintf("Gipfelhöhe: %.2f m",hmax)
set xlabel "Zeit [s]"
set yrange [0:100]
set xtics 0.5
set ytics 2
set grid xtics ytics
unset key
plot "wasserrakete.dat" using 1:2:3 title "Flughöhe [m]" with lines ls 40 lw 2, \
     "wasserrakete.dat" using 1:3 title "Geschwindigkeit[m/s]" with lines lw 2

set terminal png size 1800,1000
set output 'wasserrakete.png'
replot

Skript Download

# Wir beginnen die Differentialgleichungen in kleinen Schritten 
# bei Flughöhe=0 Geschwindigkeit=0 und Zeit=0
t = 0.0; h = 0.0; v = 0.0
# Wir rechnen bis die Rakete wieder unten ankommt also die Flughöhe h=0 wird
while h>=0 {
    v_aus = sqrt(2.0*Druck/dichte_wasser);      #  Die Strahlaustrittsgeschwindigkeit
    F_luft = A_g * abs(v)**2.0 * dichte_luft/2.0*cw;   #  Die Luftwiederstandskraft der Rakete
    dStrahlVol = A_d * v_aus*dt;                #  Volumeneinheit die in der Zeit dt 
    if (vol_fl > 0) {                           # ausgestossen wird, und damit das Füllvolumen mindert
	dmasse = dStrahlVol*dichte_wasser;        #  Masseneinheit die in der Zeit dt ausgestossen wird
	F_schub = v_aus*dmasse/dt - g*masse - F_luft;   # Die effektive Schubkraft ist die Rückstosskraft,
	                                          # vermindert um die Schwerkraft und die Luftwiederstandskraft
	masse = masse - dmasse;                   # Die Masse verringert sich um die ausgestossene Wassermasse
	Druck = Druck * (1.0/(1.0+dStrahlVol/vol_luft))**k_luft;      # Der Druck ändert sich
	                                               # mit der Volumenzunahme des Treibgases adiabatisch
	vol_luft = vol_luft + dStrahlVol;      # Das Druckluftvolumen nimmt um die ausgestossenen Wassermenge zu
	vol_fl = vol_fl - dStrahlVol
    } else {
        masse = m_geh                   # Wenn das Wasser komplett ausgestossen ist 
                                        #  bleibt die Masse konstant und es gibt keinen Schub mehr
	F_schub = - g*masse - F_luft;   # Es bleibt die Schwerkraft und die Strömungswiderstandsraft
    }
    v = v + F_schub/masse*dt;
    h = h + v*dt;              # Die zurückgelegte Flugstrecke oder bei senkrechtem Flug die Steighöhe
    t = t + dt
    print t, h, v, masse, Druck
    # , v_aus, F_luft, F_schub, vol_fl
}