Fehler und Effizienz von Lösungsmethoden für Anfangs- und Randwertproblemen aus Flugbahnmodellen

Inhaltsverzeichnis

Erklärung zur Selbständigkeit

Abstrakt

1 Einleitung
1.1 Flugbahnberechnungenfür Risikoanalysen
1.2 ZielsetzungenundHerausforderungen
1.3 Gliederungder Arbeit

2 Grundlagen
2.1 Ballistik
2.2 Koordinatensysteme
2.2.1 ErdfestesKoordinatensystem
2.2.2 GeozentrischesKoordinatensystem
2.2.3 Projektilfestes Koordinatensystem
2.3 Kräfte des STANAG 4355 - Modells
2.3.1 Hilfsgröße I: Aerodynamische Koeffizienten
2.3.2 Hilfsgröße II: LieskeVektor
2.3.3 Schwerkraft
2.3.4 Corioliskraft
2.3.5 Magnus-Kraft
2.3.6 Auftriebskraft
2.3.7 Luftwiderstand
2.3.8 Drallstabilisierung
2.4 AllgemeineszuFlugbahnen
2.4.1 Einhüllende einer Flugbahnschar
2.4.2 Flugzeitenbetrachtung
2.5 Notwendige Definitionen über Differentialgleichungen

3 Differentialgleichungssystem des Flugbahnmodells nach STANAG 4355
3.1 Differentialgleichungssystem nachSTANAG 4355
3.2 Anfangswerte
3.3 Zweiseitiges Randwertproblem
3.4 Randwertproblem mit festen Grenzen - Zeittransformation
3.5 Shooting-Methode zur Auffindung einer Lösung für ein zweiseitiges Randwertproblem

4 Verfahren und Kennzahlen zur Bestimmung der Effizienz von ODE-Solvern
4.1 Kennzahlen numerischer ODE-Solver
4.1.1 Rechenzeit(computingtime)
4.1.2 Speicherplatz(arraystorage)
4.2 Genauigkeit (accuracy,precision)
4.3 Verläßlichkeit (reliability)
4.4 Robustheit (robustness)
4.5 Definition des lokalen Diskretisierungsfehlers
4.6 Verschiedene Verfahren zur Bestimmung des globalen Fehlers
4.6.1 SummationdeslokalenFehlers
4.6.2 Klassische Methode nach Isaacson und Keller für Einschrittverfahren
4.6.3 Schätzung des globalen Fehlers durch Variation des lokalen Fehlers
4.6.4 Schätzung des globalen Fehlers durch Variation der Schrittweite
4.6.5 Defektkorrektur
4.7 Zusammenfassung der Schätzverfahren zur Bestimmung des globalen Fehlers

5 Analyse der ODE-Solver
5.1 ODE-Solver nach Numerical Recipes
5.1.1 Klassisches Runge Kutta Verfahren vierterOrdnung
5.1.2 Adaptives Runge Kutta Verfahren fünfter Ordnung mit Cash- Karp Parametern
5.1.3 Adaptives Runge Kutta Verfahren fünfter Ordnung mit Dormand Prince Parametern
5.1.4 Adaptives Runge Kutta Verfahren achter Ordnung mit Dormand Prince Parametern
5.1.5 Burlisch StoerVerfahren mit modifizierter Mittelpunkt Methode vierterOrdnung
5.1.6 Semi-implizites Extrapolationsverfahren nach Burlisch Stoer
5.1.7 Adaptives Prädiktor-Korrektor Verfahren mit Schrittweitenkontrolle
5.2 Numerische Experimente füreinigeODE-Solver
5.2.1 Rechenzeit vs. Lokaler Fehler
5.2.2 Globaler Fehler vs. Lokaler Fehler
5.2.3 Rechenzeit vs. Globaler Fehler
5.3 Entwicklung eines Verfahrens zum Vergleich der konkurrierenden Ziele Rechenzeit und Genauigkeit
5.4 Darstellungderglobalen Fehlerkurven
5.4.1 Globaler Fehler vs. Zeit
5.4.2 DarstellungderSchrittweitensteuerung
5.4.3 Polarkoordinatendarstellung des globalen Fehlers
5.4.4 Dreidimensionaler Plot des globalen Ortsfehlers
5.4.5 Fehleranteile
5.4.6 Weitere Features der Plots
5.5 Erkenntnisse über den globalen Fehler der Trajektorie

6 Optimierungsmöglichkeiten der Minimierungsalgorithmen
6.1 Guessed Values
6.1.1 Analytische Berechnung des Abschusswinkels für ein Flugbahnmodell im Vakuum bei flacher Erde
6.1.2 Performancesteigerung durch vorgeschaltete Flugbahnmodelle: Quadratisches Luftwiderstandsmodell für gestreckte Flugbahnen
6.2 Performancetestfür verschiedene Minimierungsverfahren
6.2.1 Newtonverfahren
6.2.2 Simulated Annealing
6.2.3 Broyden
6.2.4 Brent
6.2.5 BrentundNewton
6.2.6 VergleichderVerfahren
6.2.7 Zusammenfassung der Ergebnisse

7 Zusammenfassung und Ausblick
7.1 Zusammenfassung
7.2 Ausblick

8 Abbildungsverzeichnis

9 Tabellenverzeichnis

10 Literatur

11 Appendix A: Numerische Experimente und numerische Verfahren
11.1 Performancetestfür vorgeschaltete Flugbahnmodelle
11.2 Herleitung der klassischen Methode zu Abschätzung des globalen Fehlers für ODE-Solver
11.3 Anwendung des Schätzers des globalen Fehlers durch Variation der Schrittweite auf Testprobleme
11.4 Anwendung der Defektkorrektur auf Testprobleme
11.5 Das Programm ODEINT mit Runge-Kutta Quality Control und einem Runge-Kutta 4. Ordnung
11.5.1 Allgemeines zuODEINT
11.5.2 Einbettungsstrategie für RungeKuttaFamilien
11.5.3 Lokale ExtrapolationzurSchätzungdeslokalenFehlers
11.5.4 Runge-KuttaQualityControl
11.5.5 Die Routinen RK4undDERIVS

12 Appendix В: Physikalische Modelle
12.1 Herleltungder Approximation fürdleSchwerkraft
12.2 ErklärungderFormdes Luftwiderstandes
12.3 Analyse der Flugbahndauer für ein lineares Luftwiderstandsmodell
12.4 Erklärung der Auftriebskraft für Projektile

13 Appendix C: Analytische Flugbahnmodelle
13.1 Flugbahn(en) im Vakuum unter allgemeinen Bedingungen
13.2 Berechnung eines gestreckten Flugbahnmodelles mit quadratischem Luftwiderstand

Abstrakt

In dieser Arbeit werden Fehler und Effizienz von Lösungsverfahren für das Flugbahnmodell nach STANAG 4355 analysiert. Für das Anfangswertproblem werden verschiedene ODE-Solver (Klassisches Runge-Kutta, Runge-Kutta Cash-Karp, Runge-Kutta Dormand-Prince, Burlisch Stoer, Semi-implizites Burlisch Stoer, Prädiktor-Korrektor Verfahren) für typische Flugbahnen untersucht. Mit Hilfe einer globalen Fehlerbestimmung mit der Methode der Defektkorrektur und verallgemeinerten Fehlermaßen zeigt sich, dass das Runge-Kutta Dormand-Prince Verfahren fünfter Ordnung am effizientesten arbeitet. Das bisher implementierte Runge-Kutta Cash-Karp Verfahren belegt den zweiten Platz.

Es werden Möglichkeiten zur analytischen Bestimmung der für das zweiseitige Randwertproblem benötigten Startwerte aufgezeigt.

Nach Untersuchung der Minimierungsalgorithmen (Newton, Broyden, Simulated Annealing, Brent) zur Lösung des Randwertproblems konnte festgestellt werden, dass das Broyden-Verfahren dem implementierten Newton-Verfahren überlegen ist.

In this thesis the error and efficiency for methods that solve the trajectory model of STANAG 4355 are analysed. Different ODE-Solver (classical Runge-Kutta, Runge-Kutta Cash-Karp, Runge-Kutta Dormand-Prince, Burlisch Stoer, semi-implizit Burlisch Stoer, Predictor-Corrector method) for the initial value problem are investigated for representative trajectories. Using the defect correction method to determine the global error and mixed measures it is shown that the Runge-Kutta Dormand-Prince method of fifth order works most efficient. On second place Is Implemented Runge-Kutta Cash-Karp.

Analytical start values for the two-side boundary value problem are analytically determined.

When analyzing the minimization algorithms (Newton, Broyden, Simulated Annealing, Brent) that can solve the boundary value problem it was shown that the Broyden-Method is superior to the implemented Newton-Method.

1 Einleitung

In diesem Kapitel wird das Projekt „Fuze Safety Quantitative Risk Analysis" und die Hintergründe für diese Diplomarbeit kurz vorgestellt (Abschnitt 1.1). Die Zielsetzung für die Arbeit werden in Abschnitt 1.2 gegeben und Herausforderungen angedeutet. Eine Übersicht über die restlichen Kapitel der Diplomarbeit vor dem Hintergrund der Fragestellung gibt Abschnitt 1.3.

1.1 Flugbahnberechnungen für Risikoanalysen

Das Fraunhofer-Institut für Kurzzeitdynamik (Ernst-Mach-Institut, kurz: EMl) untersucht unter dem Projektnamen „Fuze Safety Quantitative Risk Analysis" (FSQRA) die Gefährdung von Personen in Überflugsszenarien mit Artillerie- und Mörsergeschossen [2].

Das Modell der Risikoanalyse basiert auf einer repräsentativen Flugbahn des Geschosses und viele tausende Flugbahnen von repräsentativen Fragmenten, die bei einem Schadensereignis durch Detonation des Geschosses entstehen.

Der mathematische Aufwand zur Berechnung der Flugbahnen, wird dabei durch das zugrundeliegende physikalische Modell bestimmt und reicht von trivial integriebarer Gleichungen bis zu gewöhnlichen nichtlinearen Differentialgleichungen zweiter Ordnung.

In dieser Arbeit wurde ein von der NATO standardisiertes modifiziertes Punkt­Masse-Modell verwendet (engl. modified point mass model, MPMM). Das »NATO STANDARDlZATlON AGREEMENT (STANAG) 4355« schreibt 5 Freiheitsgrade (Degrees of Freedom, DoF) für dieses MPMM vor [3].

1.2 Zielsetzungen und Herausforderungen

Das Ziel dieser Diplomarbeit ist es eine bedienzeitsparende numerische Lösung der Differentialgleichungen auf Basis des Flugbahnmodells nach STANAG 4355 bereitzustellen, wobei die dabei auftretende Fehler bekannt und unter Kontrolle sein sollen.

Dafür ist es notwendig alle auftretenden physikalischen Kräfte und der daraus abgeleiteten Anfangs- und zweiseitige Randwertprobleme hinreichend zu untersuchen, um so einen Eindruck für die Komplexität des Problems zu bekommen und um nachfolgende numerische Zusammenhänge zu verstehen.

Um die Genauigkeit der Flugbahn des Projektils zu beschreiben, müssen Kennzahlen eingeführt werden, um die numerischen Abweichungen zur exakten Lösung messen zu können. Hierbei ist es notwendig die verschiedenen Fehlerschätzungen für Anfangswertprobleme zu studieren, um damit die verschiedenen ODE-Solver aus der Numerical Recipes Bibliothek [4] bezüglich ihrer Genauigkeit zu vergleichen. In ähnlicher Weise muss eine Minimierungsaufgabe numerisch gelöst werden.

Zusammenfassend stellt sich die Frage, welche Fehlermaße auf welche Weise für die einzeInen numerischen Verfahren zur Lösung von Anfangs- und Randwertprobleme berechnet werden können, wie aussagekräftig sie sind und ob es möglich ist eine einzige Kennzahl zu konstruieren, die alle Verfahren vergleichbar macht.

Eine Herausforderung der Diplomarbeit ist es den globalen Fehler für Differentialgleichungsverfahren zu beschreiben, da hierbei erst a posteriori die Schwächen des Schätzverfahrens sichtbar werden. Des weiteren ist die Implementierung der gewünschten Visualisierung in C++ mit dem Graphik-Tool TeeChart herausfordernd.

Um das Potential analytischer Lösungen abschätzen zu können wird versucht werden, die auftretenden physikalischen Kräfte analytisch zu untersuchen und z.B. aus der Differentialgleichung zu entkoppeIn. Aber da selbst die Lösung einfacher analytischer Flugbahnmodelle schwierig ist, wie ein analytisch berechnetes quadratisches Luftwiderstandsmodell zeigt, werden diese Untersuchungen nur unterstützende Funktion haben können.

1.3 Gliederung der Arbeit

Im nächsten Kapitel 2 werden die notwendigen Grundlagen zur Bestimmung der Anfangs- und Randwertprobleme nach STANAG 4355 erklärt. Dabei wird auf die einzeInen physikalischen Kräfte eingegangen, die auf das Projektil wirken. Weiterhin werden die relevanten Begriffe und Definitonen von Differentialgleichungen erläutert.

In Kapitel 3 kann dann mit den vorher eingeführten Größen das Differentialgleichungssystem des Flugbahnmodells nach STANAG 4355 definiert und in eine numerisch effiziente Form für ODE-Solver gebracht werden.

Das darauffolgende Kapitel 4 beschäftigt sich mit Verfahren zur Bestimmung des Fehlers und der Effizienz numerischer Lösungsverfahren von

Differentialgleichungen, Insbesondere Rechenzelt, lokaler (Dlskretlslerungs-) Fehler, globaler Fehler und Methoden zur Bestimmung des globalen Fehlers.

Dies Ist notwendig, damit die Vielzahl von ODE-Solver miteinander verglichen werden können.

Die Anwendung eines dieser vorgestellten Verfahren, der Defektkorrektur, wird im Kapitel 5 diskutiert. In diesem Abschnitt werden auch die verwendeten ODE- Solver vorgestellt. Mithilfe der Kennzahlen in Kapitel 4 kann letztendlich das effizienteste ODE-Verfahren für das dreidimensionale Flugbahnmodell bestimmt werden.

In Kapitel 6 werden die benötigten Startwerte des zweiseitigen Randwertproblems analytisch bestimmt und die effizienteste Minimierungsmethode ermittelt.

Abschliessend werden in der Zusammenfassung die numerischen Ergebnisse zusammengetragen und zukünftige mögliche Verbesserungen der numerischen Rechenmethodenaufgezeigt. In Abbildung 1.1 ist die Strukturder Diplomarbeit zusammengefasst.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 1.1 Struktur der Diplomarbeit.

2 Grundlagen

Nach einer Einführung in die Ballistik (Abschnitt 2.1) und die verwendeten Koordinatensysteme (Abschnitt 2.2) soll dieses Kapitel vor allem die Komplexität der Problemstellung darstellen. Dazu werden alle auftretenden physikalischen Kräfte vorgestellt (Abschnitt 2.3) und relevante Begriffe über Differentialgleichungen definiert (Abschnitt 2.5). Einige Kräfte werden durch heuristische Herleitungen, oder durch weiterführende Analysen (z.B. die Corioliskraft im Abschnitt 2.3.4.1) ergänzt. In Abschnitt 2.4 werden allgemeine Eigenschaften von Flugbahnen diskutiertwie zum Beispiel die Flugbahndauer.

2.1 Ballistik

Ballistics ist the science that deals with the motion of projectiles. Modern writers divide the subject into interior, exterior, and terminal ballistics [...]. The modern science of exterior ballistics has evolved as a specialized branch of the dynamics of rigid bodies, moving under the influence of gravitational and aerodynamic forces [5].

Aristoteles, der erste Wissenschaftler, der sich mit der Ballistik beschäftigte, formulierte eine Theorie über die Bewegung von Körpern im lufterfüllten Raum [6]. Diese ersten Überlegungen waren aber mehr philosophischer als physikalischer Natur. Aristoteles Theorie wurde zwar im Spätmittelalter von J. Buridan um 1325 verworfen, aber der Durchbruch einer „echten" berechneten Flugbahn gelang erst sehr viel später Galilei in seinen Untersuchungen zu den Fallgesetzen (1638: „Untersuchungen über zwei Wissenzweige, die Mechanik und die Fallgesetze betreffend") [6]. Diese Überlegungen führten zu dem bekannten „schiefen Wurf". Dieser war Ausgangspunkt für viele mathematische und physikalische Formulierungen.

2.2 Koordinatensysteme

Die Bewegung des Projektils wird hier durch zwei Koordinatensysteme beschrieben, einerseits durch das erdfeste Koordinatensystem, das die Bewegung des Projektils von der Erde aus betrachtet und andererseits durch ein projektilfestes Koordinatensystem.

2.2.1 Erdfestes Koordinatensystem

Betrachten wir ein nach STANAG 4355 definiertes erdfestes mitrotierendes Koordinatensystem das seinen Ursprung auf der Erdoberfläche hat. Die x-Achse bezeichnet in Abbildung 2.1 hierbei die horizontale Schussrichtung. Dagegen steht die z-Achse senkrecht auf der x-Achse nach oben und die y-Achse vervollständigt das rechtsseitige Koordinatensystem. MitXend ist in Abbildung der Zielpunkt der ballistischen Kurve gemeint. Oft ist die Achsenbeschriftung auch in X¡, i = 1,2,3 dargestellt. Bei Ballistischen Problemen wird manchmal die z mit der y-Achse vertauscht. Diese Notation ist aufgrund der ersten mathematischen Flugbahnen entstanden, da diese nur zweidimensional in einem x-y-Koordinatensystem berechnet wurden. Diese Darstellung wird auch von vielen modernen Ballistikern beibehalten (zum Beispiel McCoy: Modern Exterior Ballistics (1998) [5]). In dieser Arbeit wird zu Beginn jeder Rechnung deutlich darauf hingewiesen wie die Achsen bezeichnet werden.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.1 Erdfestes Koordinatensystem.

Geozentrisches Koordinatensystem

Bei dem zweiten wichtigen, nicht rotierenden, Koordinatensystem (siehe Abbildung 2.2) fällt die z-Achse mit der Drehachse der Erde zusammen. Die x- Achse zeigt in Abschussrichtung des Geschosses und die y-Achse beschreibt demnach die rechts-links-Abweichung des Projektils. Mit X0 wird in Abbildung die Lage der Plattform bezeichnet.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.2 Koordinatensystem Im Erdmittelpunkt.

2.2.3 Projektilfestes Koordinatensystem

Als letzter dreidimensionaler euklidischer Raum wird nach STANAG 4355 das projektilfeste Koordinatensystem definiert. Dieses besitzt seinen Ursprung im Masseschwerpunkt des Flugkörpers. Die X-Achse beschreibt die Symmetrieachse des Flugkörpers, siehe Abbildung 2.3. Die Z-Achse steht dabei senkrecht auf dem Flugkörper und die Y-Achse vervollständigt das rechtshändige Koordinatensystem.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.3 Projekt il fes tes Koordinatensystem.

2.3 Kräfte des STANAG 4355 - Modells

Wie das erste Newton'sche Axiom zeigt, wirken auf das Geschoss, nach Abschuss, Kräfte:

Every body perseveres in its state of rest, or of uniform motion in a right line, unless it is compelled to change that state by forces impressed thereon [7].

(Dieses Axiom war eines der größten physikalischen Entdeckungen. Es wurde 1686 von Sir Isaac Newton in seinem Werk: „Philosophiae Naturalis Principia Mathematica" veröffentlich, aber schon 1640 von Galileo Galilei aus Versuchen an schiefen Ebenen abgeleitet [8])

Da die Geschwindigkeit des Geschosses nicht konstant bleibt, muss also (mindestens) eine Kraft auf das Projektil wirken. STANAG 4355 schreibt die in Tabelle 2.1 angegebenen Kräfte vor:

Tabelle 2.1 Kräfte des STANAG 4355 Modells [9, 10].

Abbildung in dieser Leseprobe nicht enthalten

Diese Kräfte werden im nächsten Abschnitt ausführlich betrachtet.

2.3.1 Hilfsgröße I: Aerodynamische Koeffizienten

Ein Problem zur exakten Berechnung der Flugbahn ist die Ermittlung notwendiger Koeffizienten (z.B. der cw -Wert) für die verschiedenen Geschosse:

The main difficulties are due to the need for precise calculations of the aerodynamics and flight mechanics parameters of the projectiles. Some of these parameters, such as drag and moment coefficients, have to be measured in wind tunnels or gun tunnel tests [11].

Diese experimentell bestimmten aerodynamischen Koeffizienten, die Fitting- Factoren und Korrekturpolynome, werden in Fire-Control-Input Daten (FCI- Daten) abgelegt [12]. Diese Koeffizienten sind nach STANAG 4355 nur von der Machzahl abhängig.

Die Machzahl ergibt sich aus dem Verhältnis von Geschwindigkeit und der höhenabhängigen Schallgeschwindigkeitdes Projektils.

Abbildung in dieser Leseprobe nicht enthalten

Aus den FCI-Daten können die Glieder eines abschnittweisen definierten Polynoms vom Grad < 4 abgelesen und so der aeroballistische Koeffizient Ck an der Stelle M berechnet werden:

Abbildung in dieser Leseprobe nicht enthalten

Alle in Tabelle 2.2 angegebenen aerodynamischen Koeffizienten werden auf diese Weise berechnet.

Tabelle 2.2 Aerodynamische Koeffizienten.

Abbildung in dieser Leseprobe nicht enthalten

2.3.2. Hilfsgröße II: Lieske Vektor

Der Winkel zwischen der vertikalen Achse im projektilfesten Koordinatensytem und der Flugbahn ist unter dem Namen ,,yaw of repose" bekannt [1]. Sein Formelzeichen ist α

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.4 Yaw of repose im projektilfesten Koordinatensystem (1]

Der Lieske-Vektor ae, dessen 12 - Norm eine Approximation für den yaw of repose darstellt, berechnet sich nach STANAG 4355:

Abbildung in dieser Leseprobe nicht enthalten

Mit:

- repräsentativer Durchmesser d [m] des Projektils,
- von der Höhe über Normal Null abhängige Luftdichte p(h) [kg/m3]
- Trägheitsmoment ix, welches sich auf die Symmetrie-/Längsachse des projektilfesten Koordinatensystems bezieht,
- Kreisfrequenz p(t ) [1/min],
- Beschleunigungsvektor U(t) [m/s2],
- und Geschwindigkeitsvektor V(t) [m/s].

Hierbei istzu beachten, dass V(t) = u(t)-w(t)

Abbildung in dieser Leseprobe nicht enthalten

ist. Das aktuelle Modell berücksichtigt den Wind nicht, also kann w(t) = 0 gesetzt werden.

Abbildung in dieser Leseprobe nicht enthalten

Folgen eines yaw of repose > 0o: [1,13]:

- Erhöht den Luftwiderstand und reduziert damit die Reichweite
- Es entsteht Seitendrift
- Führt zu (verstärktem) Magnus-Effekt

Die Gleichung (2.2) beschreibt einen impliziten Ausdruck, der nach der Diplomarbeit von Herrn Richter [12], mit folgender Iterationsvorschrift berechnet werden kann: erklären, dass eine niedrige Geschwindigkeit die Stabilität eines Flugkörpers senkt und dieser so mehr ins „TaumeIn" geraten kann. Während der Bewegung zurück zur Erdoberfläche nimmt die Geschwindigkeit wieder zu und das Projektil stabilisiert sich wieder.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.5 Yaw ๙ repose M80 bullet fired at 32° [1].

Schwerkraft

Eine allgemein bekannte Kraft, die auf das Geschoss wirkt ist die Schwerkraft, die für das charakteristische Aussehen der Flugbahn, in Form einer Parabel, verantwortlich ist. Es ist zu beachten, dass die Erde nicht flach und somit der konstante Vektor, bei einem erdfesten Koordinatensystem nach STANAG 4355,

Abbildung in dieser Leseprobe nicht enthalten

nicht ausreichend präzise für eine Approximation der Gravitation und die Zentrifugalkraft ist.

Vielmehr eignet sich folgende Formel nach STANAG [9], für „erdnahe" Flugkörper über der Erdoberfläche bei Idealisierung der Erde als Kugel, deren Herleitung im Appendix 12.1 zu finden ist.

Abbildung in dieser Leseprobe nicht enthalten

Hierbei wirddieauftretendeZentrifugalkraftin g0 beachtet:

Abbildung in dieser Leseprobe nicht enthalten

2.3.4 Corioliskraft

Durch die Erdrotation entsteht nicht nur die Zentrifugalkraft, sondern auch die Corioliskraft. Diese Scheinkraft beeinflusst Flugbahnen mit einer hohen Schussweite.

Für die Corioliskraft ergibt sich für das erdfeste Bezugssystem nach STANAG 4355 [9]:

Abbildung in dieser Leseprobe nicht enthalten

Dabei bezeichnet owie üblich die Kreisfrequenz der Erddrehung

Abbildung in dieser Leseprobe nicht enthalten

2.3.4.1 Koordinatentransformation zur analytischen Berechnung der Corioliskraft

Es ist möglich, die Corioliskraft in ein anderes dreidimensionales, kanonisches Koordinatensystem zu transformieren.

Abbildung in dieser Leseprobe nicht enthalten

Jeder Punkt im erdfesten System nach STANAG kann mit folgender Transformation in das geozentrische Koordinatensystem (siehe Abschnitt 2.2.2) abgebildet werden.

Abbildung in dieser Leseprobe nicht enthalten

Das erdfeste Koordinatensystem (rot) nach STANAG 4355 ist in folgender Abildung einbeschrieben:

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 2.6 Beispiel einer Koordinatensystemtransformation.

MitTransformationssformel (2.7)werden folgende FormeIn errechnet:

Abbildung in dieser Leseprobe nicht enthalten

Offensichtlich kann mit Hilfe der obigen FormeIn jeder Punkt im rotfarbige Koordinatensystem in das blaue transformiert werden.

Bemerkung:

Die Inverse der gesamten Transformation berechnet sich durch Anwendung der inversen Transformationsmatrizen in umgekehrter Reihenfolge [14], wobei: R-1(a) = R(-α)

T-1 entsteht durch das Einsetzen des negativen Translationsvektors.

Abbildung in dieser Leseprobe nicht enthalten

Im Koordinatensystem, das seinen Ursprung im Erdmittelpunkt hat, gilt für die Kraftkomponente der Corioliskraft nach Athen [15] (-2ωΥ,2ωΧ, 0)T . Hierbei sei ω dieWinkelgeschwindigkeit derErde.

Wären zu jedem Zeitpunkt die parametrischen Kurven des Geschosses ohne Corioliskraft bekannt mit (ξ(ί ), η(ί ), r(t )), dann beschreibt das Differentialgleichungssystem (2.8) die Bewegung des Projektils (nach Athen [15]).

Abbildung in dieser Leseprobe nicht enthalten

Trivialerweise wird in Richtung der Erdachse durch die Corioliskraft keine Veränderung herbeigeführt.

Wird die komplexe Variable

Abbildung in dieser Leseprobe nicht enthalten

definiert, entsteht aus den ersten beiden Zeilen von (2.8): [15]

Abbildung in dieser Leseprobe nicht enthalten

und der Lösung:

Abbildung in dieser Leseprobe nicht enthalten

Dieses Integral lässt sich leicht in Real- und Imaginärteil zerlegen.

Abbildung in dieser Leseprobe nicht enthalten

Somit ist es möglich, die Corioliskraft getrennt von anderen Kräften zu betrachten. Leider ist diese Vorgehensweise kein Durchbruch für eine analytische Lösung, da unser Flugbahnmodell auch ohne die Corioliskraft nicht geschlossen lösbar ist. Auch für die numerische Behandlung ist diese Methode nicht geeignet, da zwei Koordinatentransformationen (in jedem Zeitschritt) durchgeführt werden müssen.

2.3.5 Magnus-Kraft

Der Magnus-Effekt wurde 1852 von dem deutschen Physiker und Chemiker Heinrich Gustav Magnus entdeckt. Dieser Effekt tritt in Erscheinung, wenn sich ein rotierendes Objekt in einem strömenden Medium befindet. Dadurch entsteht eine Kraft senkrecht zur Strömung.

Bei Geschossen entsteht der Magnus-Effekt dadurch, dass das Geschoss die vorbeistreichenden Luftteilchen im Sinne seiner Rechtsrotation mitreißt. Die Geschossachse liegt oberhalb der Flugbahntangente (Strömungsrichtung der Luft). [...] Der so entstehende Druckunterschied auf beiden Seiten muß im Allgemeinen in einer Linksabweichung des Geschosses in Erscheinung treten. Es kann jedoch auch der Fall eintreten, dass die geschilderten Abweichungen nach der entgegengesetzten Seite erfolgen, nämlich dann, wenn die Geschossspitze unterhalb der Bahntangete liegt [15].

Dadurch kann das Projektil sich aber auch instabil verhalten:

The spin velocity required to stabilize the projectiles tends to induce Magnus effects which can lead to dynamic instabilities that influence the flight path [16].

Folgende Grafik macht den Magnus-Effekt in Abhängigkeit seiner Kreisfrequenz Ω und die daraus resultierende Kraft deutlich. Dargestellt ist ein Zylinder. Der Schnitt verläuft senkrecht zur Symmetrieachse. Die Rotation ist im Uhrzeigersinn:

Abbildung in dieser Leseprobe nicht enthalten

Nach STANAG 4355 geht die Magnus-Kraft mit folgender Formel In das MPMM ein (siehe auch McCoy [5]):

Abbildung in dieser Leseprobe nicht enthalten

Mit der mittleren Querschnittsfläche[Abbildung in dieser Leseprobe nicht enthalten]die restlichen Größen wurden in der Tabelle 2.2 und Gleichung (2.2) erklärt. Gleichung (2.11) zeigt, dass die Magnus-Kraft zunemmt, je mehr die Geschossachse von der Flugbahntangente abweicht. Mit Hilfe von (2.2) erkennt man noch eine direkte Proportionalität zur Kreisfrequenz p und zum Trägheitsmoment.

2.3.6 Auftriebskraft

Durch das Luftpolster, dass auf der Vorderseite eines symmetrischen Geschosses entsteht ein Druck auf das Projektil. Diese Kraft wird als die Auftriebskraft bezeichnet. Die Auftriebskraft entsteht aufgrund der Abweichung zwischen der Nase des Projektils und seiner Flugbahn (siehe Abschnitt 2.3.2) [15, 18].

Abbildung in dieser Leseprobe nicht enthalten

Für idealisierte symmetrische Geschosse deren Nase immer in Richtung der Flugbahn zeigen fällt diese Kraft weg.

Durch den Auftrieb wird die Reichweite von Geschossen erhöht, dabei sei vermerkt, dass eine Auftriebserhöhung auch immer eine Widerstandserhöhung bewirkt [19]. Nach STANAG 4355 ist diese Kraft wie folgt definiert:

Abbildung in dieser Leseprobe nicht enthalten

Mit dem Fitting Faktor fL, der von der Abgangsgeschwindigkeit und damit von derTreibladung abhängt.

Im Anhang ist unter die Auftriebskraft für Flugzeuge und asymetrische gewölbte Flugkörper erklärt, bei denen diese einfache Anschauung nicht ausreicht.

2.3.7 Luftwiderstand

Eine Erklärung für die Form des Luftwiderstandes ist in Appendix 12.2 zu finden. Nach STANAG 4355 wird folgende Gleichung verwendet:

Abbildung in dieser Leseprobe nicht enthalten

mitdem dimensionslosen Fitting Factor fD .

Hier fallen sofort die Unterschiede zu der Herleitung im Anhang auf:

Der cw-Wert wurde offensichtlich mit dem Ausdruck[Abbildung in dieser Leseprobe nicht enthalten]gleichgesetzt. Dieser hängt nun von der Machzahl, dem Lieske-Vektor und dem konstanten yaw drag factor Qd , der einen Fitting Factor für die Machzahl darstellt (siehe STANAG 4355 [20]), ab.

2.3.8 Drallstabilisierung

Bei Glattrohrkanonen fällt diese Krafteinwirkung weg, da diese Projektile flügelstabilisiert sind.

Der Gegensatz zur Glattrohrkanone sind bei der Zugrohrkanone die Geschosse drallstabilisiert. Bei diesem Geschützt sind Nuten spiralförmig in den Lauf eingefräst. Die geschossführenden Flächen nennt man "Felder", die Nuten "Züge". Diese spiralförmig eingefrästen Züge geben dem Geschoss einen Drall der zur Flugbahnstabilisierung dient.

Wie sich der Drall auf die Flugbahn eines Baseballs auswirkt wird in Abbildung 2.9 deutlich:

Abbildung in dieser Leseprobe nicht enthalten

Die durchgehende Linie, genzeichnet eine Flugbahn ohne Drall und einer Abschlagsgeschwindigkeit von 100 mph. Die kurz gestrichelte Linie zeigt die Reichweite des Baseballs bei 1000 Umdrehungen pro Minute. Die Kurve in der Mitte der zwei anderen zeigt dagegen die Weite bei 2000 1/min.

Dadurch wird deutlich, dass die Kreisfrequenz einen großen Einfluss auf die Reichweite des Flugkörpers besitzt.

In STANAG 4355 wird für die Berechnung der Frequenz folgende Formel als Rechnungsgrundlage genommen:

Abbildung in dieser Leseprobe nicht enthalten

2.4 Allgemeines zu Flugbahnen

Es gibt einige Eigenschaften, die sich über Flugbahnmodelle aussagen lassen, ohne die Differentialgleichungen zu kennen.

2.4.1 Einhüllende einer Flugbahnschar

Eine Flugbahnschar besitzt unter gleichbleibenden Bedingungen und einer konstanten Abschussgeschwindigkeit v0 > 0 für variable Schußwinkel φ, θ eine einhüllende Oberfläche. Die Einhüllende selbst ist eine eindeutige Flugbahn. Innerhalb dieser Einhüllenden kann jeder Punkt von zwei Flugbahnen und außerhalb von keiner getroffen werden (siehe Molitz und Strobel [22]).

Anmerkung:

Diese Eigenschaft wird verwendet, falls zwei Geschosse simultan auf ein Ziel einschlagen sollen. Dafür wird zuerst ein Hochschuss (die Flugbahn mit dem größeren Abschusswinkel) und danach der Direktschuss abgefeuert.

2.4.2 Flugzeitenbetrachtung

Eine weitere interessante Eigenschaft der klassischen Flugbahnmodelle ist, dass die Flugbahndauer bis zur maximalen Höhe größer ist als die Zeit, die das Geschoss von seinem höchsten Punkt zurück zur Erde benötigt. Es spielt nach F. Bauer [23] keine Rolle welches Modell betrachtet wird, solange für die Force Function des Geschosses f(v) gilt:

Abbildung in dieser Leseprobe nicht enthalten

Im Anhang wird unter 12.3 ein Beweis für ein lineares Luftwiderstandsmodell geführt.

2.5 Notwendige Definitionen über Differentialgleichungen

In diesem Abschnitt wird auf die wichtigsten und nur relevanten Begriffe über Differentialgleichungen eingegangen.

1.) Definition einer Differentialgleichung

Jede Gleichung, die mindestens einen Differentialquotienten enthält heißt Differentialgleichung (kurz: DGL). Sie heißt von n-ter Ordnung, falls die höchste vorkommende Ableitung von Ordnung n ist. [24]

2.) Differentialgleichungssystem

Allgemein spricht man über ein „System von Differentialgleichungen erster Ordnung" wenn folgende Form vorliegt:

Abbildung in dieser Leseprobe nicht enthalten

Dabei bezeichnet x, die Ableitung der Funktion x, (t) nach t.

Abbildung in dieser Leseprobe nicht enthalten

gesetzt, so kann man ein solches System auch in Vektordarstellung wie folgt aufschreiben (siehe Hermann [25]):

Abbildung in dieser Leseprobe nicht enthalten

Diese platzsparende Notation wird im folgendem verwendet, da nach Hermann jede Differentialgleichung n-ter Ordnung in ein solches Differentialgleichungssystem (2.14) mit n Unbekannten durch Substitution transformiert werden kann.

3.) Anfangswertproblem und Randwertproblem

Im Allgemeinen wird zwischen Anfangswertproblemen (AWP), bei denen die Lösung nur am linken Rand t = t0 eines Intervalls vorgegeben ist (Anfangswert) und zweiseitge Randwertproblemen (RWP), bei denen für die zwei Ränder t = t0, t = tend Bedingungen vorgegeben sind, unterschieden.

4. ) Lineare und nichtlineare Differentialgleichungen

Eine lineare Differentialgleichung liegt vor falls die Differentialgleichung in die normierte Form

Abbildung in dieser Leseprobe nicht enthalten

überführbar ist. Ein lineares Randwertproblem liegt vor, wenn sowohl die DGLs als auch die Randbedingungen linear sind [25].

Ansonsten istdie Differentialgleichung nichtlinear und verliert viele „schöne" Eigenschaften. Nichtlineare Differentialgleichungen besitzen im allgemeinem meist keine eindeutige Lösung und oft keine analytisch (geschlossene) Lösung. Auch Beweise die die Existenz von mindestens einer Lösung garantieren sind oft sehr schwierig zu führen. Für unsere nichtlineare Flugkurven ist aber bei einer hinreichenden Startgeschwindigkeit mindestens eine Lösung des zweiseitigem RWPs vorhanden, so dass wir die Annahme der Existenz einer Lösung für die nachfolgenden Kapitel vorraussetzen.

3 Differentialgleichungssystem des Flugbahnmodells nach STANAG 4355

In diesem Kapitel wird das Differentialgleichungssystem für das Flugbahnmodell nach STANAG 4355 aufgestellt (Abschnitt 3.1 bis Abschnitt 3.4). Durch geschickte Transformation in Abschnitt 3.4 wird der numerische Aufwand der DGLs reduziert. Um eine spätere Analyse durch die ODE-Solver zu ermöglichen ist eine genaue Kenntnis der zugrundeliegenden Differentialgleichung unabdingbar. Außerdem wird kurz der Algorithmus zur Lösung des Randwertproblems vorgestellt (Abschnitt 3.5).

3.1 Differentialgleichungssystem nach STANAG 4355

Betrachten wir die Bewegung des Projektils in drei Dimensionen, dann sei[Abbildung in dieser Leseprobe nicht enthalten]der Geschwindigkeitsvektorund[Abbildung in dieser Leseprobe nicht enthalten]der

Ortsvektor zum Zeitpunkt[Abbildung in dieser Leseprobe nicht enthalten]bzgl des Koordinatenssystems. Für die Beschleunigung des Geschosses X(t) gilt nun:

Abbildung in dieser Leseprobe nicht enthalten

Wobei die Kräftef(i), i = i,...,5 in den Abschnitten 2.3.3 bis 2.3.7 eingeführt wurden.

Wie bekannt, können diese drei DGLs zweiter Ordnung in ein DGL-System mit sechs Gleichungen transformiert werden.

Wird nun auch der Drall aus Abschnittp(t) aus Abschnitt 2.3.8 in das DGL- System erster Ordnung eingeführt:

Abbildung in dieser Leseprobe nicht enthalten

3.2 Anfangswerte

Für eine Lösung des DGL-Systems müssen noch die Anfangsbedingungen der Form

Abbildung in dieser Leseprobe nicht enthalten

bestimmt werden. Sind diese Anfangsbedingungen geeignet gewählt, sowird ein Zielpunkt

Abbildung in dieser Leseprobe nicht enthalten

zum Zeitpunkt tend getroffen. Im erdfesten Koordinatensystem nach STANAG 4355 entspricht die Koordinate der Plattform.

Abbildung in dieser Leseprobe nicht enthalten

mit[Abbildung in dieser Leseprobe nicht enthalten]Höhe der Plattform über normal Null.

Abbildung in dieser Leseprobe nicht enthalten

Das Projektil beginnt seinen Flug jedoch nicht am Plattformstandpunkt selbst, sondern dort, wo der Flugkörper die Plattform, z.B. den Lauf verlässt. Mit einer plattformpezifischen Länge lW beschreibtfolgende Gleichung die Projektilkoordinaten zum Zeitpunkt t0

Abbildung in dieser Leseprobe nicht enthalten

Mit

QE = vertikalem Abschusswinkel (oft als «bezeichnet) zur X-Achse.

ΔΑΖ - horizontalem Abschusswinkel (oft auch als ^bezeichnet) zur X-Achse. Somit steht lediglich noch die Anfangsbedingung für die Kreisfrequenz des

Dralls aus. Diese ¡stzum Zeitpunkt t0 gemäß STANAG 4355 wie nachstehend definiert:

Abbildung in dieser Leseprobe nicht enthalten

wobei u0 für die Anfangsgeschwindigkeitdes Flugkörper(ohne den Einfluss des Windes) steht. tcd ist die Strecke, die für eine Umdrehung im Rohr durchlaufen werden muss.[Abbildung in dieser Leseprobe nicht enthalten]ist die dafür benötigte Zeit.[Abbildung in dieser Leseprobe nicht enthalten]wobei ls als

Dralllänge bezeichnet wird (siehe Richter [12]).

Somit gilt für die Anfangswerte:

Abbildung in dieser Leseprobe nicht enthalten

3.3 Zweiseitiges Randwertproblem

Da wir ein gewähltes Ziel zum Zeitpunkt tend treffen wollen, erhalten wir ein zweiseitiges Randwertproblem. Nach STANAG 4355 werden nur die Abschusswinkel QE (vertikaler Abschusswinkel) und ΔΑΖ (horizontaler Abschusswinkel) variiert um das Ziel zu treffen.

Das gesamte Gleichungssystem mit der linken und rechten Randbedingung bei [Abbildung in dieser Leseprobe nicht enthalten]tend sei somit:

Abbildung in dieser Leseprobe nicht enthalten

Die Flugzeit beträgt:

Abbildung in dieser Leseprobe nicht enthalten

Meistsetzt man [Abbildung in dieser Leseprobe nicht enthalten] und [Abbildung in dieser Leseprobe nicht enthalten]

3.4 Randwertproblem mit festen Grenzen - Zeittransformation

Hier wird gezeigt, wie mit Hilfe einer linearen Zeittransformation das zweiseitige Randwertproblem mit variablen Grenzen in ein zweiseitiges Randwertproblem mit festen Grenzen umgewandeltwerden kann. Dann werden für die Parametrisierung der Anfangsbedingung bereits alle zu bestimmenden Größen an der linken Grenze verwendet. Ein RWP mit festen Grenzen erfordert weniger Aufwand als die numerische Umsetzung eines zweiseitigen Randwertproblems mit variablen Grenzen.

Wir definieren:

Abbildung in dieser Leseprobe nicht enthalten

Wie man sieht, ergibt sich für [Abbildung in dieser Leseprobe nicht enthalten] die Abschusszeit und für[Abbildung in dieser Leseprobe nicht enthalten]die Auftreffzeit.

Um nun die Zeitskalierung auf das DGL-System anwenden zu können, muss zunächstzu ddt ein äquivalentes Differential gefunden werden:

Sei im folgenden nun f eV c Mn ein Platzhalter für eine beliebige zeitabhängige Vektorfunktion.

Dann gilt (siehe Richter [12]):

Abbildung in dieser Leseprobe nicht enthalten

Und mit:

Abbildung in dieser Leseprobe nicht enthalten

folgt:

Abbildung in dieser Leseprobe nicht enthalten

Aus Übersichtsgründen werden Größen mit transformierter Zeit wie folgt geschrieben.

Abbildung in dieser Leseprobe nicht enthalten

Werden diese Überlegungen nun verwendet, erhält man folgendes DGL- System:

Abbildung in dieser Leseprobe nicht enthalten

(3.2)

3.5 Shooting-Methode zur Auffindung einer Lösung für ein zweiseitiges Randwertproblem

Nach STANAG 4355 werden die Abschusswinkel geeignet variiert um eine Lösung für das zweiseitige Randwertproblem zu ermitteIn. Dazu wird ein Startwinkelpaar geschätzt oder durch den Benutzer vorgegeben. (Wie dies im einzeInen geschieht, wird im Kapitel 6.1 erklärt).

Es müssen zudem die Randbedingungen (boundary conditions), die sich als

Abbildung in dieser Leseprobe nicht enthalten

schreiben lassen, erfüllt werden [4].

Das folgende Vorgehen ist unter dem Namen „Shooting Verfahren" oder „Einfach-Schießverfahren" bekannt. Auf dem linken Rand werden die geschätzten Startwerte eingesetzt. Dann wird das Differentialgleichungssystem mit einem ODE-Solver gelöst und ein Fehlermaß (Score) für den rechten Rand an fendermittel. Dieser Score beschreibt im wesentlichen eine Distanzfunktion zwischen dem gewünschten Ziel und dem tatsächlichen Auftreffpunkt. Nun werden die Abschusswinkel durch einen Minimierungsalgorithmus (z.B. das Newton-Verfahren) geschickt variiert um so einen besseren Score für den nächsten Schuss zu erhalten.

Dies wird so lange iteriert, bis ein Abbruchkriterium (Fehlerschranke oder maximale Anzahl an Iterationen) erfüllt ist. In Abbildung 3.1 ist der Programmablauf dargestellt und in Abbildung 3.2 die Iteration selbst.

Abbildung in dieser Leseprobe nicht enthalten

Abbildung 3.1 Allgemeiner Programmablauf für einen Minimierungsalgorithmus zur Ermittlung einer Lösung des zweiseitigem RWPs.

Abbildung in dieser Leseprobe nicht enthalten

4 Verfahren und Kennzahlen zur Bestimmung der Effizienz von ODE-Solvern

In diesem Kapitel werden die erforderlichen Kennzahlen zur Bestimmung der Effizienz von ODE-Solver vorgestellt (Abschnitt 4.1). Nach einer Definition des lokalen Fehlers in Abschnitt 4.5 werden einige Verfahren zur Fehlerschätzung des globalen Fehlers eingeführt (Abschnitt 4.6). Am Ende dieses Kapitels werden alle Verfahren nocheinmal in einer Tabelle im Abschnitt 4.7 zusammengefasst.

4.1 Kennzahlen numerischer ODE-Solver

Da das Differentialgleichungssystem nun bestimmt ist kann es numerisch gelöst werden. (Der Versuch der analytischen Lösung scheitert an der „echten" Nichtlinearität des DGLs-Systems. Eine Vereinfachung der Differentialgleichung ist ebenfalls sehr aufwändig wie z.B. die Behandlung der Corioliskraft, Abschnitt 2.3.4, zeigt).

Differentialgleichungsverfahren zur Lösung von expliziten ODEs sind nach P. Deuflhard [26] vergleichbar bezüglich:

4.1.1 Rechenzeit (computing time)

Ein wichtiges Kriterium, das allerdings erst durch eine Gegenüberstellung von Genauigkeit vs. Rechenzeit eine Aussage über ein optimales Verfahren ermöglicht.

Abbildung in dieser Leseprobe nicht enthalten

Für nichtsteife Verfahren gilt:

Abbildung in dieser Leseprobe nicht enthalten

[...]