# Kapitel 3_c: Raketengrundgleichung
# © 2015 Friedrich U. Mathiak,
# mathiak@mechanik-info.de
# 
> restart:
# Wir berechnen zunächst aus der Raketengrundgleichung die Geschwindigkeit v(t) . Die Differenzialgleichung 1. Ordnung für v(t) lautet (Buch TM 3)
> ode1:=diff(v(t),t) - mu0/(m0-mu0*t)*vr + g;
# Die Geschwindigkeit zum Startzeitpunkt t = t0 ist
> ABV:=v(t0) = v0; #Anfangsbedingung
;
> dsolve({ode1,ABV},v(t)): assign(%):
> v:=unapply(v(t),t);
# Die Beschleunigung ergibt sich aus der Ableitung der Geschwindigkeit nach der Zeit t
> a:= unapply(diff(v(t),t),t);
# Die Steighöhe x(t) ermitteln wir durch Integration der  Differenzialgleichung 
> ode2:=diff(x(t),t)-v(t);
# Die Anfangsbedingung für die Steighöhe ist:
> ABX:=x(t0)=x0;
> dsolve({ode2,ABX},x(t)): assign(%):
> x:=simplify(unapply(x(t),t));
# Wir prüfen die Anfangsbedingungen 
> v(t0); x(t0);
# Ist mL die Leermasse der Rakete zur Brennschlusszeit tB, dann ist
> tB:=(m0-mL)/mu0;
# Die größte Steighöhe xmax und die größte Geschwindigkeit vmax errechnen sich bei Brennschluss zu
> xmax:=x(tB);
> vmax:=v(tB);
# Wir schreiben eine Prozedur, die die obigen Beziehungen automatisch auswertet:
> Proc_Calc_17:=proc(vr,mu0,m0,mL,a0,v0,x0,t0)
> #----------------------------------------------------------------------
> # Eingabe: 
> #          vr: Relativgeschwindigkeit der heißen Brenngase in [m/s]
> #         mu0: Konstanter Massenstrom dm/dt                in [kg/s]
> #          m0: Startmasse der Rakete                       in [kg]
> #          mL: Leermasse                                   in [kg]
> #          v0: Anfangsgeschwindigkeit                      in [m/s]
> #          x0: Anfangshöhe                                 in [m]
> #          a0: Anfangsbeschleunigung                       in [m/s^^2]
> #          t0: Startzeit                                   in [s]
> # Ausgabe:  
> #          a: Beschleunigung  (global vereinbart)          in [m/s^^2]
> #          v: Geschwindigkeit (global vereinbart)          in [m/s]
> #          x: Flughöhe        (global vereinbart)          in [m]
> #         tB: Brennschlusszeit                             in [s]
> #          S: Schubkraft                                   in [N]
> #         a0: Beschleunigung  bei Start                    in [m/s^2]
> #         aB: Beschleunigung  bei Brennschluss             in [m/s^2]
> #         vB: Geschwindigkeit bei Brennschluss             in [m/s]
> #         xB: Steighöhe       bei Brennschluss             in [m]
> #---------------------------------------------------------------------
> local g,t;
> global a,v,x,tB,S,G0,aB,vB,xB;
> description "Bewegungsvorgang einer senkrecht startenden Rakete in Erdnähe ohne Berücksichtigung des Luftwiderstands";
> g:=9.81;
> a:= unapply(mu0*vr/(-mu0*t+m0)-g + a0,t);
> v:= unapply(-vr*ln(-mu0*t+m0)-g*t+v0+vr*ln(-mu0*t0+m0)+g*t0,t);
> x:= unapply((-vr*t+vr*m0/mu0)*ln(-mu0*t+m0)+(vr*t-vr*m0/mu0)*ln(-mu0*t0+m0)+vr*t-(1/2)*g*t^2+(g*t0+v0)*t+x0-vr*t0+(1/2)*g*t0^2-(g*t0+v0)*t0,t);
> tB:=(m0-mL)/mu0;   #Brennschlusszeit
> S := mu0*vr;       #Schubkraft
> G0:= g*m0;         #Startgewicht
> aB:= a(tB);        #Beschleunigung  bei Brennschluss
> vB:= v(tB);        #Geschwindigkeit bei Brennschluss
> xB:= x(tB);        #Höhe bei Brennschluss
> #Ausgabe der Ergebnisse
> print(`Startgewicht G0 in [N] = `,G0);
> print(`Brennschlusszeit tB in [s] = `,tB);
> print(`Schubkraft S in [N] = `,S);
> print(`Flughöhe bei Brennschluss xB in [m] = `,xB);
> print(`Geschwindigkeit bei Brennschluss vB in [m/s] =`,vB);
> print(`Beschleunigung beim Start a0 in [m/s**2] = `,a(t0));
> print(`Beschleunigung bei Brennschluss aB in [m/s**2] = `,aB);
> end proc:
# Beispiel 3-4: Start einer einstufigen Rakete
# Mit den Werten des Beispiels 
> vr:=2.22*10^3; mu0:=1.5*10^4; m0:=2.95*10^6; mL:=1.0*10^6;a0:=0.; v0:=0.; x0:=0.; t0:=0.;
# rufen wir die Prozedur auf und erhalten:
> Proc_Calc_17(vr,mu0,m0,mL,a0,v0,x0,t0);
# Wir geben die Zustandsgrößen grafisch aus:
> para:=t=t0..t0+tB,titlefont=["ARIAL", 15],labeldirections=["horizontal","vertical"],axesfont = ["HELVETICA","ROMAN",8],labelfont = ["HELVETICA", 12],axes=boxed,gridlines=true: 
> p1:=plot(x(t-t0),para,labels=["Zeit t [s]", "Flughöhe x [m]"]):
> p2:=plot(v(t-t0),para,labels=["Zeit t [s]", "Geschwindigkeit v [m/s]"]):
> p3:=plot(a(t-t0),para,labels=["Zeit t [s]", "Beschleunigung a [m/s**2]"]):
> plots[display](Matrix(1,3,[p1,p2,p3]));
# Aufgrund der geringer gewordenen Masse erreicht die Rakete  in den letzten
> dt:=10.:
# Sekunden vor Brennschluss  einen Zuwachs der Steighöhe in [m] von
> deltax:=x(tB-t0)- x(tB-dt-t0);
# sowie einen Zuwachs in der Geschwindigkeit in [m/s] von 
> deltav:=v(tB-t0)- v(tB-dt-t0);
# Beispiel 3-6: Berechnung des Bewegungszustandes einer zweistufigen Rakete. Die Nutzmasse mN  beträgt  120 kg.
> mN:=120.0;
# Technische Daten der  1. Stufe:
> vr1:=1.1*10^3; mu01:=15.;m01:=1100.0; mL1:=50.; mT1:=m01-mL1; 
# Technische Daten der  2. Stufe:
> vr2:=1.4*10^3; mu02:=3. ;m02:= 290.0; mL2:=20.; mT2:=m02-mL2; 
# Wir simulieren im ertsen Schritt den Flug der Rakete vom  Start bis zum Brennschluss der 1. Stufe. Die Rakete startet zum Zeitpunkt t0 = 0 aus der Ruhelage mit v0 = 0 und x0 = 0. Dann gilt:
> a0:=0.; v0:=0.; x0:=0.; t0:=0.;
> m0:=m01 + m02 + mN; mL:=mL1 + m02 + mN;
# Wir rufen die Prozedur  Proc_Calc_17  mit diesen Parametern auf:
> Proc_Calc_17(vr1,mu01,m0,mL,a0,v0,x0,t0);
# Wir geben die Zustandsgrößen vom  Start bis zum Brennschluss der 1. Stufe grafisch aus:
> para:=t=0..tB,titlefont=["ARIAL", 15],labeldirections=["horizontal","vertical"],axesfont = ["HELVETICA","ROMAN",8],labelfont = ["HELVETICA", 12],axes=boxed,gridlines=true: 
> p11:=plot(x(t),para,labels=["Zeit t [s]", "Flughöhe x [m]"]):
> p21:=plot(v(t),para,labels=["Zeit t [s]", "Geschwindigkeit v [m/s]"]):
> p31:=plot(a(t),para,labels=["Zeit t [s]", "Beschleunigung a [m/s**2]"]):
> plots[display](Matrix(1,3,[p11,p21,p31]));
# Für die grafische Darstellung der Zustandsgrößen beider Brennstufen werden die Zustandsgrößen der ersten Stufe  abgespeichert.
# Wir starten die 2. Stufe. Hier gilt nach dem Abwurf der 1. Stufe
> a0:=aB; v0:=vB; x0:=xB; tB1:=tB;
> m0:=m02 + mN; mL:= mL2 + mN;
> Proc_Calc_17(vr2,mu02,m0,mL,a0,v0,x0,0);
> tB2:=tB:
# Wir geben die Zustandsgrößen nach Abwurf der 1. Stufe bis zum Brennschluss der 2. Stufe grafisch aus:
> para:=t=tB1..tB1+tB2,titlefont=["ARIAL", 15],labeldirections=["horizontal","vertical"],axesfont = ["HELVETICA","ROMAN",8],labelfont = ["HELVETICA", 12],axes=boxed,gridlines=true: 
> p12:=plot(x(t-tB1),para,labels=["Zeit t [s]", "Flughöhe x [m]"]):
> p22:=plot(v(t-tB1),para,labels=["Zeit t [s]", "Geschwindigkeit v [m/s]"]):
> p32:=plot(a(t-tB1),para,labels=["Zeit t [s]", "Beschleunigung a [m/s**2]"]):
> plots[display](Matrix(1,3,[p12,p22,p32]));
# Es folgt eine zusammenfassende Ausgabe der Zustandsgrößen für den Zeitraum beider Brennstufen:
> p1:=plots[display]([p11,p12]):
> p2:=plots[display]([p21,p22]):
> p3:=plots[display]([p31,p32]):
> plots[display](Matrix(1,3,[p1,p2,p3]));
> 
;
