# Kapitel 4_b: 
# Fallschirmsprung mit Luftwiderstand (Turbulente Dämpfung)
# © 2015 Friedrich U. Mathiak,
# mathiak@mechanik-info.de
# 
> restart:
# Beispiel 4-5:
# 
# 
# Wir lösen zunächst das Problem in der Freifallzone. In diesem Bereich ist der Koeffizient  k = cwÅAÅrho/(2m)  konstant. 
# Die zugehörige Differenzialgleichung für die Geschwindigkeit lautet 
> ode2:=diff(x(t),t,t)-k*diff(x(t),t)^2 + g;
# Der Anfangswert ist
> AW:=x(0)= h,D(x)(0)=0:
# Als Lösung erhalten wir: 
> dsolve({ode2,AW}); assign(%);
> loe:=allvalues(x(t));
> x(t):=loe[2];
# Die Ableitung von x(t) nach der Zeit t ergibt die Geschwindigkeit 
> v(t):=simplify(diff(x(t),t));
# Grenzgeschwindigkeit
> vgr:=limit(v(t),t=infinity);
# Die Ableitung der Geschwindigkeit v(t) nach der Zeit ergibt die Beschleunigung
> a(t):=simplify(diff(v(t),t));
# Die Öffnung des Fallschirms erfolgt in der Höhe h1. Die Freifallzeit ermitteln wir aus der Beziehung
> loe:=solve(h1-x(t),t);
# Wir erhalten 2 Lösungen. Mit den Werten unseres Beispiels 
> h:=3000.; h1:=1500.; g:=9.81;
> cw:=1.3; A:=0.5;rho:=1.2; m:=85.;
> k:=cw*A*rho/2/m;
# können wir weiter konkretisieren
> loe;
> ts:=loe[1];
# Wir geben die Zustandsgrößen x(t) , v(t) und a(t) für die Freifallzone mit den Werten des Beispiels aus.
> para:=t=0..ts,axes=boxed,gridlines=true,titlefont=["ARIAL", 15],labeldirections=["horizontal","vertical"],axesfont = ["HELVETICA","ROMAN",8]:
> p1:=plot(x(t),para,title = "\nWeg\n",labels=["Zeit t", "Weg x(t) [m]"]):
> p2:=plot(v(t),para,title = "\nGeschwindigkeit\n",labels=["Zeit t", "Geschwindigkeit v(t) [m/s]"]):
> p3:=plot(a(t),para,title = "\nBeschleunigung\n",labels=["Zeit t", "Beschleunigung a(t) [m/s^2"]):
> plots[display](Matrix(1,3,([p1,p2,p3])));
# Die Grenzgeschwindigkeit beträgt
> vgr;
> restart:
# Im Fogenden wird die nummerische Lösung vom Absprung bis zur Landung bereitgestellt. Wir notieren zunächst die maßgebende Differenzialgleichung 2. Ordnung
> ode2:=diff(x(t),t,t) - cw*rho/2/m*A*diff(x(t),t)^2 + g;
# und beachten, dass die Schattenfläche A eine Funktion von x ist.:
> Alin:=A2-(A2-A1)*(x(t)-h1+Dh)/Dh;
> A:=piecewise(x(t) >= 0 and x(t) < h1-Dh ,A2, x(t) >= h1-Dh and x(t) < h1,Alin,A1);
# Es sind folgende Anfangswerte zu berücksichtigen
> AW:=x(0)=h,D(x)(0)=0;
# Mit den Werten unseres Beispiels 
> cw:=1.3; g:=9.81; rho:=1.2; m:=85.;
> A1:=0.5; A2:=15.; h:=3000.; h1:=1500.; Dh:=200.;k:=cw*A*rho/2/m;
# Wir stellen zur Kontrolle die Funktion A(x)   grafisch dar:
> plot(A,x=0..h,y=0..A2,gridlines=true,axes=boxed,gridlines=true,titlefont=["ARIAL", 15],labeldirections=["horizontal","vertical"],axesfont = ["HELVETICA","ROMAN",8],labels=["Weg x [m]", "A(x) [m^2]"]);
# Wir lösen nun die Differenzialgleichung  ode2  nummerisch:
> loe:=dsolve({ode2,AW},numeric,output=listprocedure);
# Zur Kontrolle der Anfangswerte geben wir die Lösung für t = 0 aus:
> loe(0);
# Die Berechnung endet für x(t = te) = 0, denn dann ist der Springer gelandet. Um diesen Zeitpunkt te zu ermitteln, geben wir die Lösung für die Höhe x(t) zunächst grafisch aus.
> plots[odeplot](loe,[t,x(t)],t=0..220,axes=boxed,gridlines=true,titlefont=["ARIAL", 15],labeldirections=["horizontal","vertical"],axesfont = ["HELVETICA","ROMAN",8],title= "\nHöhe x(t)\n",labels=["Zeit t [s]", "Weg  x [m]"]);
# Die Nullstelle von x(t)  wird iterativ ermittelt; sie liegt augenscheinlich zwischen  200 < t < 210. Wir beginnen die Suche bei t = 200 und erhöhen diesen Wert jeweils um 1%, bis wir im Ergebnis einen Vorzeichenwechsel feststellen. 
> t0:=200.:
> for i to 5 do
>   t0:=t0*1.01:
>   print(i,loe[2](t0));
> end do:
# Wir wählen te = 204,36 und prüfen damit die Randbedingung für x(t = te):
> te:=204.36:
> loe(te);
# Mit diesem Wert ist eine ausreichende Genauigkeit erreicht, und die Zustandsgrößen werden in einem Matrixplot ausgegeben.
> para:=t=0..te,axes=boxed,gridlines=true,titlefont=["ARIAL", 15],labeldirections=["horizontal","vertical"],axesfont = ["HELVETICA","ROMAN",8]:
> p1:=plots[odeplot](loe,[t,x(t)],para,title = "\nWeg\n",labels=["Zeit t", "Weg x(t) [m]"]):
> p2:=plots[odeplot](loe,[t,diff(x(t),t)],para,title="\Geschwindigkeit v(t)\n",labels=["Zeit t","Geschwindigkeit v(t) [m/s]"]):
> p3:=plots[odeplot](loe,[t,k*diff(x(t),t)^2 - g],para,title="\Beschleunigung a(t)\n",labels=["Zeit t","Beschleunigung a(t) [m/s^2]"]):
> plots[display](Matrix(1,3,[p1,p2,p3]));
# Zum Vergleich mit dem analytischen Ergebnis der Freifallzone berechnen wir den Zeitpunkt, zu dem der Fallschirm in der Höhe h = 1500 m  geöffnet wird. Für den Faktor k  mit A = A1 = 0.5 m2 gilt:
> k:=cw*rho*A1/2/m;
> tau:=arccosh(exp(k*(h-h1)));  #Dimensionslose Zeit
;
# Der Springer besitzt dort die Geschwindigkeit
> v:=-sqrt(g/k)*tanh(tau);
# Die nummerische Rechnung liefert praktisch dasselbe Ergebnis:
> loe(tau/sqrt(g*k));
> 
;
