Wenn Ihr Browser die Formeln nicht richtig anzeigt, verwenden Sie bitte das PDF.

Bestimmung von Geschossflugbahnen mit dem Analogrechner

Rainer Glaschick, Paderborn
rainer@glaschick.de
2015-09-24 korr. 2021-10-06

Hintergrund

Im 2. Weltkrieg waren die Berechnung von Geschossflugbahnen eine militärisch wichtige Aufgabe; der erste elektronische Rechner (ENIAC) war zu diesem Zweck gebaut. Es gab damals noch keine elektronischen Analogrechner, und der mechanische differential analyzer war wohl zu ungenau.

Die nachfolgend beschriebene Lösung ist übrigens auf einem einfachen Analogrechner auch nicht befriedigend.

Da der Luftwiderstand eines Geschosses keiner einfachen Funktion folgt, ist, wie so häufig, eine geschlossene Integration nicht möglich, so dass es naheliegend ist, dies Beispiel auf einem der üblichen Analogrechner zu programmieren, insbesondere um den Einsatz des Funktionsgenerators für den Luftwiderstand zu zeigen.

Die Ergebnisse sind nicht spektakulär und eher für ein Fachpublikum interessant.

Ausgangssituation

Reed und GornJuncosa beziehen ihre Ausführungen auf den ENIAC und geben folgendes Formelsystem an:

Hierin sind:

Gegenüber dem Original ist hier `E` noch mit `x`, `y` und `z` indiziert, denn es muss selbstverständlich die Komponente in der jeweiligen Koordinatenrichtung verwendet werden, also (mit dem Winkel `phi` der Flugrichtung zur Horizontalen):

Da die Verzögerung immer der Geschwindigkeit entgegen gerichtet ist, ist der Vorzeichenwechsel des Sinus am Scheitel richtig, sofern E(v) immer positiv ist.

Nicht berücksichtigt ist in der Formel die von der Temperatur und der Höhe abhängige Dichte der Luft, sowie ein durch gezogenen Läufe erreichter Drall, der zwar ein Zylindergeschoss stablilisiert, aber auch ein Drehmoment benötigt, um das Geschoss von der Abschussrichtung in die Landerichtung zu drehen.

Der Luftwiderstand `E` hat die Dimension einer Beschleunigung, d.h. ein Quotient aus Kraft und Masse; die Funktion ist damit spezifisch für ein Geschoss und berücksichtigt so seine Form und Masse.

Empirisch wurde herausgefunden, dass der Luftwiderstand unterhalb Schallgeschwidigkeit `v_s` gut durch eine quadratische Abhängigkeit dargestellt werden kann:

wobei die — geschossspezifische — Konstante `K` die Dimension `m^-1` hat. Da die Mündungsgeschwindigkeiten üblicherweise über der Schallgeschwindigkeit liegen, wird `K` als Funktion von `v` dargestellt (die für kleine Geschwindigkeiten nahezu konstant ist):

Hierzu zeigt Athen auf S.19 folgenden Graphen für `K(v)`, wobei für die Y-Werte — wie sich aus dem Vergleich mit den Tabellen ergibt — auf den Bereich 0..10 skaliert wurden, also relativ zueinander sind:

Lediglich die Kurve für Siacci zeigt die horizontale Gerade für eine quadratische Abhängigkeit. Offenbar setzt der Anstieg des Widerstandes gegenüber der quadratischen Abhängigkeit bereits deutlich vor der Schallgeschwindigkeit ein.

Im Anhang von Athen findet sich auch eine Tabelle, deren Werte `K(v)` von Siacci folgenden Graph ergeben, in dem auch die Verzögerung `E(v)=K(v)*v^2` des Geschosses dargestellt ist:

Ersichtlich geht der Luftwiderstand eines Geschosses bei Schallgeschwindigkeit von einer Parabel in eine Grade mit ins Negative verschobenem Ursprung über.

Athen zitiert folgende von Siacci entwickelte Formel:

Deren Graph ist gleichfalls in der obige Graphik gezeigt und liegt so genau auf den Tabellenwerten, dass die Tabelle vermutlich daraus berechnet wurde. Eine Analyse der Formel (Details im Anhang) zeigt, dass diese in der Tat eine in eine Gerade übergehende Parabel ist; es könnte also sein, dass die mit Siacci bezeichnete Kurve in obiger Darstellung eine theoretische ist, während — wie auch an anderer Stelle erwähnt — die von Krupp empirisch ist und nicht in Form eines analytischen Ausdrucks vorliegt.

Bei Lorenz auf S.92 wird eine Tabelle für das S-Geschoss (Infanteriegewehr) angegeben (zwei Spalten von mir hinzugefügt)1:

Hierbei sind:

Hier werden die Werte der 2. und 3. Spalte auf den (unbekannten) Querschnitt F des Geschosses bezogen.

Interessant ist die vierte Spalte, die dimensionierte Werte für die Verzögerung enthält und gleichfalls einen in etwa linearen Anstieg oberhalb der Schallgeschwindigkeit zeigt:

Offenbar handelt es sich bei dem S-Geschoss um ein kleines Geschoss mit einer geringen Masse pro Querschnitt, so dass die Verzögerungen recht groß sind; bei einer Mündungsgeschwindigkeit von 500m/s² wird das Geschoss überschlägig während der ersten Sekunde auf 250m/s² abgebremst, und hat dabei weniger als 500m zurückgelegt, während Artilleriegeschosse Flugzeiten um die 20s bei und Reichweiten von 6km und einer Endgeschwindigkeit von 250m/s² erreichen (Athen, S.239).

Da die Kurve der Verzögerung zwar bei kleinen Geschwindigkeiten in eine Parabel übergeht, dieser Bereich aber gar nicht ausgenutzt wird, weil nur eine signifikante Aufprallgeschwindigkeit zerstörend wirkt, könnte eine erste Annäherung mit einer der Geschwindigkeit proportionalen Verzögerung versucht werden; jedoch würde hier eine positive Beschleunigung unterhalb der Schallgeschwindigkeit resultieren, die tunlichst vermieden werden sollte.

Man kann daher, wie in der Formel von Siacci, ausnutzen, dass sich der Ausdruck `sqrt(k+x^2)` für große x (`x>sqrt(k)`) einer Geraden nähert, und hat dann eine einfache Formel, die jedoch für Endgeschwindigkeiten unterhalb der Schallgeschwindigkeit zu große Verzögerungen bewirkt; ihr Graph ist oben bereits eingezeichnet.

Auf S.239 gibt Athen am Beispiel einer Schusstafel an, dass ein Artilleriegeschoss bei einem Abschusswinkel von 12° nach einer Flugzeit von 17s mit einer Endgeschwindigkeit von 257 m/s ein Ziel in 5,4 km Entfernung erreicht. Folgende Überlegungen dazu:

Angenommen, die Mündungsgeschwindigkeit sei 700m/s. Die y-Komponente (in der Senkrechten) der Abfluggeschwindigkeit ergibt sich durch Multiplikation mit sin 12° zu 145m/s². Sie wird durch die Gravitation und den Luftwiderstand im Scheitelpunkt auf Null reduziert, wobei die Bremsung auch in y-Richtung sich nach der Gesamtgeschwindigkeit und nicht nach deren y-Komponente bestimmt. Wenn der Scheitelpunkt (angenommen) nach 7s erreicht ist, dann ist der Beitrag der Graviation eine Geschwindigkeitsreduktion von 70m/s, verbleiben 75m/s für den Luftwiderstand, dessen y-Komponente allein schon deshalb abnimmt, weil im Scheitelpunkt das Geschoss allein in x-Richtung fliegt. Für eine grobe Abschätzung sei daher die y-Komponente des Luftwiderstands zu Anfang doppelt so groß wie die Erdanziehung, also 20m/s², was für den Gesamtwiderstand in Flugrichtung 96m/s² ergibt. Verbleiben 10s bis zu Aufprall, der mit 98m/s erfolgen würde, wenn weiter kein Luftwiderstand vorhanden wäre, der die Wirkung der Gravitation schwächt. Eine Endgeschwindigkeit von 257m/s bedeutet daher einen Auftreffwinkel von 22°.

Für die x-Richtung ist die Komponente der Abfluggeschwindigkeit 684m/s und die Auftreffgeschwindigkeit 240m/s (bei den angenommen 22°). Die Geschwindigkeit nimmt zunächst exponentiell ab, solange der Widerstand proportional der Geschwindigkeit ist; im Scheitelpunkt ist der Widerstand nur in x-Richtung wirksam. Wäre der Widerstand konstant, dann ist bei 17s für eine Geschwindigkeitsverminderung von 444m/s eine mittlere Verzögerung von 26m/s² gegeben, also anfänglich deutlich größer, wenn auch eher weniger als die aus der y-Richtung geschätzten 98m/s². Vorbehaltlich der Simulation ist offenbar bei Artilleriegeschossen auf Grund der relativ zum Querschnitt größeren Masse von Verzögerungen im Bereich bis 100m/s², also der 10-fachen Erdbeschleunigung auszugehen.

Aufstellung der Differentialgleichungen

Für die exemplarische Darstellung auf dem Analogrechner sollen Windeffekte, die Abhängigkeit vom geographischen Ort usw. nicht berücksichtigt werden und auch nur die x- und y-Richtung (Schussrichtung und Höhe) betrachtet werden; solange die Erdrotation nicht berücksichtigt wird, bleibt, da die Anfangswerte Null sind, der Endwert von `z` (Rechts-Links-Abweichung) auch Null.

Ohne Luftwiderstand ist die Rechenschaltung trivial: Die x-Richtung ist proportional der Zeit, benötigt also allenfalls einen Integrator, und die y-Richtung ist das Integral über das Integral der Konstanten, der Gravitation, ist also eine Parabel, deren Darstellung mit zwei Integratoren erfolgen kann. Dies ist die nachfolgende Rechenschaltung, wobei die Verbindungen von den Multiplizierern zu den Integratoren nicht verbunden werden und die Skalierung geändert werden muss.

Die Ableitungen werden als eigene Variablen dargestellt und die Abhängigkeit von der Zeit nicht geschieben:

Damit lauten die zu lösenden Gleichungen:

Dabei werden als weitere Variable die Geschwindigkeit `v` in Flugrichtung, nach der sich der Luftwiderstand bestimmt, sowie der Winkel `alpha` verwendet, der den Winkel der Fluggeschwindigkeit zur Horizontalen angibt. Da `sin alpha` im absteigenden Ast negativ wird, hat die Komponente des Luftwiderstands weiterhin dasselbe Vorzeichen wie die Geschwindigkeit und bremst das Geschoss.

Erfreulicherweise kann damit (für `v != 0`) der Winkel `alpha` einfach eliminiert werden:

Setzt man `R(v^2) = (E(v)) / v` (mit der Dimension 1/s), d.h. arbeitet man die Wurzel und die Division durch `v` in die Luftwiderstandsfunkion ein, so ergibt sich:

In Integralschreibweise ergibt das mit der Integration der Geschwindigkeiten zu den Wegen:

Damit werden — ausser Integrationen und Additionen — zwei Quadrierungen und zwei Multiplikationen benötigt, und nur eine einzige Funktionstabelle für `R(v^2)`.

Variante 1

Auf dem EAI Mini-AC mit drei Multiplizieren und zwei Funktionsgebern können somit die beiden Quadrierungen durch jeweils einen Multiplizierer und eine Funktionstabelle bewirkt werden.

Das ergibt folgende Rechenschaltung:

Da die Anfangwerte ja Geschwindigkeitskomponenten sind, ist die Summe ihrer Quadrate immer ≤1, so dass für den Addierer kein Überlauf zu befürchten ist, weil zudem die maximalen Geschwindigkeiten am Anfang auftreten2. Ungünstig ist, dass die Quadrate der Geschwindigkeiten gebildet werden, die bei kleinen Geschwindigkeiten sehr klein werden und damit kleine Eingangswerte für die Widerstandsfunktion bedeuten; dies wird am Schluss noch genauer dargestellt. Da die kleinen Geschwindigkeiten erst am Schluss der Flugbahn auftreten, wirken sich die Ungenauigkeiten nur dort aus und werden nicht propagiert.

Für kleine Geschwindigkeiten mit quadratischer Widerstandsfunktion wird `R(v^2) = K*v`, also `R(x) = K * sqrt(x)`, so dass, wenn dies vorteilhaft wäre, anstellle eines Funktionsgenerators eine Quadratwurzel verwendet werden könnte.

Variante 2

Zwar ist die vorherige Version übersichtlich, hat aber den Nachteil, dass die Funktion `R(v^2)` im unteren Bereich ausgeprägt ist und daher schwer stabil einstellbar ist.

Es ist aber möglich, die Länge `v=sqrt(x^2 + y^2)` eines Vektors wegen `v=x + y^2/(v+x)` mit einem Quadrat und einer Division zu ermitteln, und dann `R(v)=(E(v))/v` zu verwenden:

Anmerkungen

Während die Anfangswerte für die Position einfach als (0,0) gesetzt werden können, sind die Anfangswerte für `v_x` und `v_y` die Komponenten der Anfangsgeschwindigkeit `v_0` mit dem Abschusswinkel `alpha`. Da wenige Analogrechner feste Funktionsgeneratoren für beispielsweise den Cosinus (oder eine andere Einrichtung zur Umwandlung von Polarkoodinaten in cartesische, z.B. gekoppelte sin/cos Potentiometer) haben, müssen diese Anfangswerte extern berechnet und dann eingestellt werden.

In der Literatur werden Lösungen für kleine Geschwindigkeiten mit quadratischer Abhängigkeit des Widerstandes beschrieben:

WassGarner benötigt den Funktionsgenerator nur optional (und dann für `K(v)`), benutzt einen Arcus-Tangens-Resolver, um den Winkel zu bestimmen, und Sinus und Cosinus-Funktionen für die Komponenten in x- und y-Richtung.

KornKorn verwendet anstelle der Komponenten in x- und y-Richtung den Winkel und die Geschwindigkeit in Flugrichtung und benötigt daher den Sinus und den Cosinus für die Komponenten der Gravitation in Flugrichtung und als Drehmoment.

Dimensionierung und Skalierung

Da die elektronischen Analogrechner so beschriftet sind, dass die Rechengrößen als Zahlen zwischen -1 und +1 auftreten, wird der Begriff der Maschineneinheit (z.B. 1ME=10V) nicht benötigt; zur Skalierung wird ein Faktor gleicher Dimension wie die zu skalierende Größe verwendet, mit dem die Rechengrößen multipliziert werden, um die Realgrößen zu erhalten, z.B. 100km für eine Länge. Der Skalierungsfaktor ist meist das Maximum der Absolutwerte der physikalischen Größe, durch ein Dach gekennzeichnet, während die dimensionslose Rechenvariable im Bereich -1..+1 mit einem Unterstrich bezeichnet wird:

Beispielsweise wird in einer Gleichung nach der Substitution der Variablen durch — in diesem Fall — `hat a_x` dividiert:

Ersichtlich hat `ul k` dieselbe Dimension wie `k`, in diesem Fall sind beide dimensionslos.

Rechnung in Echzeit heisst dabei lediglich, dass `hat t = 1 s` und ausser den richtigen Dimensionen keine Wirkung hat.

Einsetzen der noch symbolischen Skalierungsfaktoren in obige Formeln ergibt:

wobei

Die maximale Geschwindigkeit beim Abfeuern ist Mach 3, also ist der Maximalwert der Geschwindigkeit und damit der Skalierungsfaktor für Geschwindigkeiten (zweckmäßig in beiden Richtungen gleich):

Für die Flugzeiten könnten 20 sec ausreichen, die in 2 sec ablaufen sollen:

Die maximale Reichweite könnte bei 20km liegen:

Wegen der Gravitation wird zwar die Höhe geringer sein, aber für die Darstellung auf dem Plotter sind unterschiedliche Skalierungen nicht sinnvoll, daher ist auch in y-Richtung:

Als maximale Beschleunigung wird die ca. zehnfache Erdbeschleunigung genommen:

Das ergibt:

Die Anfangswerte für die Geschwindigkeiten in x- und y-Richtung sind dabei für einen Abschusswinkel von 30° und einer Anfangsgeschwindigkeit von 700m/s (Mach 2):

Soll die Simulation in 5s statt in 2s erfolgen, so wird

Bleibt die Funktion `ul R(ul v^2)`. In der folgenden Tabelle wurden von der Siacci-Funktion zehn Werte ausgewählt, die Bereiche auf 0..1 skaliert und in der zweiten Spalte das Quadrat der ersten und der vierten Spalte das Produkt der ersten mit der dritten eingefügt:

Die entsprechenden Graphen seien hier wiederholt:

Die benötigte Funktion `R(v^2)` ergibt sich dann als Graph der vierten in Abhängigkeit von der zweiten Spalte:

Die Koeffizienten für die beiden zusätzlich dargestellten Annäherungen wurden durch Ausprobieren gefunden. Der Kreuzungspunkt der beiden Funktionen liegt bei `0.082=sqrt(0.286)` und damit — wie bereits oben erwähnt — deutlich unterhalb der Schallgeschwindigkeit.

Für die Funktionen ergibt sich aus der Theorie nach Siacci für `E(v)` (mit der Schallgeschwindigkeit `v_s`):

und damit wegen `R(v^2) = (E(v))/v` bzw. `R(x) = (E(sqrt(x)))/sqrt(x)`:

Dass die Funktion `R(v^2)` zwischen 0 und 0.1 sehr ausgeprägt ist, hängt auch damit zusammen, dass Argument das Quadrat der Geschwindigkeit ist; damit sind die Rechenwerte immer wesentlich kleiner als die Geschwindigkeiten; dies ist — wie oben erwähnt — ungünstig. Allerdings solle ja die Endgeschwindigkeiten noch recht hoch sein; eine Endgeschwindigkeit von 0.25 entsprechend 250m/s² bedeutet den Wert von 0.0625, also gut 6% vom Rechenbereich und sollte somit noch brauchbare Ergebnisse liefern.

Um den Flug im luftleeren Raum zu simulieren, sind andere Skalierungen notwendig. Nimmt man hier — wegen der größeren Höhe — einen Abschusswinkel von 45° und eine Anfangsgeschwindigkeit von 710m/s, dann sind die Anfangsgeschwindigkeiten in x- und y-Richtung 500m/s. Die Flugzeit ist das Doppelte der Zeit, die die Graviation benötigt, um die Anfangsgeschwindigkeit in y-Richtung auf Geschwindigkeit Null im Scheitelpunkt zu reduzieren, also etwa 2*50s=100s. Da der Flug in x-Richtung ungebremst ist, sind das 50km, so dass die Skalierungsfaktoren sind:

Soll die Simulation 5s dauern, so wird:

Auf einem EAI MiniAC mit der zweiten Variante (Quadrierer und Dividier) entstanden die folgenden Ergebnisse:

Die Knicke bei der Berechnung der Geschwindigkeit sind nicht plausibel; die Ursache wurde noch nicht gefunden, könnte aber daran liegen, dass das Quadrat durch den Funktionsgenerator gebildet wird und dieser nicht sehr stabil ist. Vertauschen von x und y (d.h. via y² und nicht via x²) ergab schlechtere Ergebnisse.

Anhang

Siaccis Formel

Die Formel gewinnt an Klarheit, wenn sie ein wenig umgeformt wird und — ohne den Verlust an Genauigkeit im Detail zu prüfen — die Konstanten rundet:

Die Formel besteht aus drei Teilen:

Ein Blick auf die Graphen ist aufschlussreich:

Athens Musterbahn

Auf S.75 gibt Athen ein mit den zuvor entwickelten Formeln die Zahlen für eine Flugbahn an, die er als Musterbahn bezeichnet. Es handele sich um (Artillerie-) Geschoss Kaliber 10 cm mit einem Abschusswinkel von 12° und einer Mündungsgeschwindigkeit von 381,6 m/s, dass nach einer Flugzeit von 17s mit einer Endgeschwindigkeit von 257 m/s ein Ziel in 5,4km Entfernung erreicht (alle Werte in m, s, m/s und Grad):

Die Daten wurden in der folgendenden Tabelle umsortiert und die Geschwindigkeit nach dem Scheitelpunkt negativ notiert:

Graphisch dargestellt (links Weg in m, rechts Geschwindigkeit in m/s):

Die Flugbahn ist nicht stark durch den Luftwiderstand beeinflusst; auf dem absteigenden Ast ist die Beschleunigung durchschnittlich 8,67 m/s² (191,5 m/s auf 22,08s), also nur um 1,1 m/s² kleiner als die Erdbeschleunigung.

Zudem ist die Flugzeit im absteigenden Ast größer als im aufsteigenden, obwohl eher eine Verkürzung der Parabel im absteigenden Ast erwartet wird. Athen gibt nicht an, ob diese berechnete Bahn experimentiell bestätigt wurde.

Literatur

In älteren Büchern wird häufig für eine Kraft auch die Einheit kg anstelle von kp (=9,81 kgm/s² = 9,81 N) verwendet; ob eine Kraft oder Masse gemeint ist, kann nur aus dem Zusammenhang erschlossen werden.

Das Buch von Siacci ist aufgeführt, obwohl es nicht direkt zitiert wird (es konnte nicht einfach beschafft werden), weil es von Lorenz und Athen zitiert wird.

Reed:
Harry L. Reed, Jr.: Firing Table Computations on the Eniac. Proceedings of the 1952 ACM National Meeting (Pittsburgh), pp. 103-106 (1952)
GornJuncosa:
S. Gorn; M.L. Juncosa: On the Computational Procedures for Firing and Bombind Tables. BRL Report No. 889, Aberdeen Proving Ground, Maryland USA (1954)
Athen:
Hermann Athen: Ballistik. Von Quelle und Meyer Leipzig (1941)
Lorenz:
H. Lorenz: Ballistik. R. Oldenbourg Berlin (1917)
Siacci:
F. Siacci: Corso di Balistica. 3 Bände, Rom 1870 bis 1884, 2. Auflage Turin 1888 (französisch: Balistique extérieure. Paris 1892) (deutsch: Ballistik und Praxis. Berlin 1882)
Jackson:
A. S. Jackson: Analog Computation. McGraw-Hill (1960)
KornKorn:
Granino A. Korn, Therea M. Korn: Electonic Analog Computers. McGraw-Hill (1956) (auch in deutscher Übersetzung)
WassGarner:
C.A.A. Wass, K.C. Garner: Introduction to Electonic Analogue Computers. Pergamon Press (1965)

1Dezimalkommas und Dezimalpunkt werden beide gleichermaßen verwendet, folglich wird weder Punkt noch Komma zur Gruppierung von Dezimalziffern verwendet
2Durch eine Rückführung vom Ausgang auf den Eingang die Eingänge mit dem Faktor 0.5 zu bewerten, ist daher nicht sinnvoll, weil es die Eingangswerte zu dem variablen Funktionsgenerator unnötig klein macht.