In dieser Arbeit werden Fehler und Effizienz von Lösungsverfahren für Flugbahnmodelle analysiert.
Die ausführlich diskutierte Theorie der Fehler- und Effizenzbestimmung von Lösungsmethoden für Anfangswertprobleme kann aber auf jede beliebige ODE übertragen werden.
Für das Anfangswertproblem des Flugbahnmodelles 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.
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.
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 PunktMasse-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
[...]
- Arbeit zitieren
- Fabian Suhrke (Autor:in), 2009, Fehler und Effizienz von Lösungsmethoden für Anfangs- und Randwertproblemen aus Flugbahnmodellen, München, GRIN Verlag, https://www.grin.com/document/186617
-
-
-
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen. -
-
Laden Sie Ihre eigenen Arbeiten hoch! Geld verdienen und iPhone X gewinnen.