Ein Verfahren zur Schätzung von Formfaktoren, Magnus-Koeffizienten usw. aus Flugbahndaten Prof. Dr. Andreas Rudolph Universität der Bundeswehr München Institut für Mathematik und Informatik FB BW Werner-Heisenberg-Weg 39 85577 Neubiberg Email: [email protected] 8. September 2015 Inhaltsverzeichnis 1 Überblick 1.1 Einführung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Mathematische Grundlagen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Diverse Algorithmen 2.1 Ein einführendes Beispiel - das Kaliber .308 . . . . . . . . . . . . . . . 2.2 Das Kaliber .300 Win Mag . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Das Kaliber .338 Lapua Magnum . . . . . . . . . . . . . . . . . . . . . 2.4 Schätzung des Formfaktors mithilfe der Geschwindigkeit . . . . . . . . 2.5 Schätzung des Formfaktors mithilfe der Energie . . . . . . . . . . . . . 2.6 Schätzung des Formfaktors mithilfe der Zeiten . . . . . . . . . . . . . . 2.7 Verallgemeinerungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7.1 Schätzung der Parameter aus Flughöhe und Seitenabweichung 2.7.2 Schätzung der Parameter aus der Geschwindigkeit . . . . . . . 2.7.3 Schätzung der Parameter aus der Geschossenergie . . . . . . . 2.8 Bereitstellung der Daten . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Auswahl eines adäquaten Luftwiderstandsmodells . . . . . . . . . . . . 2.10 Erweiterungen des Verfahrens . . . . . . . . . . . . . . . . . . . . . . . v . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 2 11 11 15 17 18 21 21 23 26 27 28 29 29 30 vi INHALTSVERZEICHNIS Kapitel 1 Überblick 1.1 Einführung Bei allen Modellen zur Beschreibung der Flugbahn eines Geschosses liegt die Newton’sche Mechanik zugrunde, d. h. die Gesamtheit aller auf das Geschoss einwirkenden Kräfte F ist gleich deren Masse m multipliziert mit seiner Beschleunigung a: F = m·a = m d2 x dt 2 wenn x = x(t) die Ortskoordinate des Geschosses zum Zeitpunkt t darstellt. Hierbei unterstellt man, dass das Projektil als Massenpunkt dargestellt werden kann - man beschreibt strenggenommen die Bahn seines Schwerpunktes1 . In erster Näherung vernachlässigt man zuerst die Ausdehnung des Geschosses - Rotationseffekte wie die Präzession, die Nutation, den Magnus-Effekt und den Spin-Dämpfungseffekt betrachtet man dann in einem zweiten Schritt. Formuliert man nun die Kräfte, die auf das Geschoss einwirken, und löst das obige Differentialgleichungssystem zweiter Ordnung F = m·a = m d2 x dt 2 dann erhält man damit die Flugbahn des Geschosses in sogenannter parametrisierter Form damit kann man in Abhängigkeit von der Zeit alle Größen bestimmen, die für das sportliche, aber auch für das militärische Schießen von Interesse sind. 1 Anscheinend wird der Lagrange-Formalismus und die Hamilton-Jacobi Theorie bisher nicht verwendet - die Möglichkeit allgemeinerer Koordinatensysteme bleibt unbenommen. 1 KAPITEL 1. ÜBERBLICK 2 1.2 Mathematische Grundlagen Betrachtet man das vorliegende Differentialgleichungssystem zweiter Ordnung F = m·a = m d2 x dt 2 dann stellt sich die Frage, ob und unter welchen Bedingungen ein derartiges System lösbar ist, und wenn es lösbar ist, wann die Lösung eindeutig ist. Dazu schreibt man das System um in ein System erster Ordnung - man erhält dadurch prinzipiell ein System der Gestalt y′ = f(t, y) dem noch entsprechende Anfangsbedingungen y(t0 ) = y0 hinzugefügt werden. Hierbei ist nun y1 (t) f1 (t, y1 , . . . , yn ) .. y(t) = ... und f(t, y) = . yn (t) fn (t, y1 , . . . , yn ) Gesucht ist damit ein Vektor y(t) von Funktionen y1 (t), . . . , yn (t), der das obige Differentialgleichungssystem erster Ordnung samt Anfangsbedingungen erfüllt. Beispiel 1 Wir betrachten das sogenannte Vakuum-Modell, d. h. der Luftwiderstand wird vernachlässigt. Damit haben wir als Differentialgleichungssystem mẍ = 0 mÿ = −mg wenn nur die Flugweite x = x(t) und die Flughöhe y = y(t) betrachtet werden sollen. Der Punkt soll dabei die Ableitung nach der Zeit bedeuten. Zusätzlich stellt man Anfangsbedingungen x(0) y(0) vx (0) vy (0) = = = = x0 y0 v0 cos(θ ) v0 sin(θ ) wenn v0 die Anfangsgeschwindigkeit des Geschosses ist und g die Erdbeschleunigung. Dieses Differentialgleichungssystem zweiter Ordnung überführen wir durch Division mit m in 1.2. MATHEMATISCHE GRUNDLAGEN 3 ẍ = 0 ÿ = −g Wir definieren die Vektoren x y y= ẋ , ẏ ẋ ẏ f(t, y) = 0 −g und erhalten damit x ẋ ẏ d y = dt ẋ 0 ẏ −g oder ẏ = f(t, y) ein Differentialgleichungssystem erster Ordnung. Damit stellt sich die Frage nach der Existenz und Eindeutigkeit einer Lösung eines derartigen Differentialgleichungssystems erster Ordnung. Ein Satz, der die Existenz und Eindeutigkeit von Lösungen von Differentialgleichungssystemen erster Ordnung generell sicherstellt, ist der bekannte Satz von Picard und Lindelöf - dieser besagt: Satz 1 Wir betrachten das Differentialgleichungssystem y′ = f(x, y), y(u) = v Die Abbildung f sei in dem (n + 1)-dimensionalen Zylinder | x − u |< a, k y − v k< b stetig und beschränkt, sie soll weiter in diesem eine Lipschitz-Bedingung erfüllen, d. h. es gelte k f(x, y) k≤ A, k f(x, y1 ) − f(x, y2 ) k≤ M k y1 − y2 k Dann gibt es für das Differentialgleichungssystem y′ = f(x, y) KAPITEL 1. ÜBERBLICK 4 genau eine Lösung y = y(x), welche die Anfangsbedingung y(u) = v erfüllt. und für | x − u |< α existiert, wobei b α = min{a, } A ist. B EWEIS Siehe zum Beispiel Kamke (1969), S. 74, Satz 1 und S. 104, Satz 3. Beispiel 2 Betrachten wir nochmals das obige Beispiel des Vakuum-Modells, dann besagt der Satz von Picard-Lindelöf, dass genau eine Lösung existiert. Dies bedeutet, dass ein Vektor y(t) existiert mit x(t) y(t) y(t) = vx (t) vy (t) der die Anfangsbedingungen x(0) = x0 , y(0) = y0 , vx (0) = v0 cos(θ ) und vy (0) = v0 sin(θ ) erfüllt, wenn v0 die Anfangsgeschwindigkeit darstellt. Diese Lösung kann man durch entsprechende Integrationen oder durch einen Potenzreihenansatz erhalten2 - es ergibt sich: x(t) = v0 cos(θ )t y(t) = y0 + v0 sin(θ )t − 21 gt 2 und durch Differentiation vx (t) = v0 cos(θ ) vy (t) = v0 sin(θ ) − gt Damit hängt die Lösung von einem Parameter ab - dem Erhebungswinkel θ . Fordert man nun, dass in einer vorgegebenen Entfernung x(tZiel ) gelten soll y(tZiel ) = 0 dann stellt sich damit die Frage, wie der Winkel θ gewählt werden muss, um dies zu erreichen. Insbesondere stellt sich die Frage, ob in jedem Fall ein entsprechender Winkel existiert. 2 Es gibt auch die Möglichkeit, die Picard-Lindelöf-Iterationen durchzuführen - eine Frage des Geschmacks 1.2. MATHEMATISCHE GRUNDLAGEN 5 Generell gilt hierzu folgender Satz 2 Sei G(x, y) ein Gebiet, c = (c1 , . . . , ck ) ein k-dimensionaler Vektor, C sei eine offene Menge in dem c-Raum und G × C das Kartesische Produkt von G und C - d. h. der Bereich der Punkte x, y, c eines (n + k + 1)-dimensionalen Raumes mit x, y ∈ G und c ∈ C. Die Abbildung f(x, y, c) sei in G × C stetig, so dass für jedes c ∈ C durch jeden Punkt x, y ∈ G mindestens eine Integralkurve der Differentialgleichung y′ = f(x, y, c) geht. Es wird weiter vorausgesetzt, dass es auch nur eine solche von Rand zu Rand laufende Integralkurve gibt, die mit y = Y(x, u, v, c) bezeichnet wird. Dann gilt folgendes: (a) Wenn die Funktion y = Y(x, u, v, c) an einer Stelle x0 , u0 , v0 , c0 existiert und wenn α ≤ x ≤ β ein abgeschlossenes Intervall ist, dass x0 , u0 als innere Punkte enthält und in dem Y(x, u0 , v0 , c0 ) existiert, dann existieren die Funktionen y = Y(x, u, v, c) in diesem Intervall auch für alle hinreichend nahe an u0 , y0 , c0 gelegenen u, y, c. (b) Die Funktionen y = Y(x, u, v, c) sind an jeder Stelle ihres Existenzbereiches stetige Funktionen ihrer Argumente x, u, v, c. B EWEIS Siehe Kamke (1969), S. 116, Satz 2. Beispiel 3 Wir betrachten wiederum das Vakuum-Modell mit der Lösung x(t) = v0 cos(θ )t y(t) = y0 + v0 sin(θ )t − 21 gt 2 Offensichtlich hängen die Lösungen in stetiger Weise von t und θ ab. Gibt man nun ein x(tZiel ) = x0 vor und fordert, dass y(tZiel ) = 0, dann liegt es nahe, das folgende Funktional zu minimieren: 1 2 ))2 h(θ ) = (x(tZiel ) − v0 cos(θ )tZiel )2 + (y(tZiel ) − (y0 + v0 sin(θ )tZiel − gtZiel 2 Da es eine maximale Flugzeit gibt und ebenso einen maximalen Winkel für θ , wird die Funktion h(.) auf einem kompakten Intervall zu betrachten sein. KAPITEL 1. ÜBERBLICK 6 Aus Stetigkeitsgründen wird somit ein entsprechender Winkel existieren müssen. Um den notwendigen Erhebungswinkel θ zu berechnen - will man in einer vorgegebenen Entfernung x0 treffen - eliminieren wir zuerst den Parameter t, indem wir die erste Gleichung nach t auf lösen t= x(t) v0 cos(θ ) und dann t in der zweiten Gleichung eliminieren: 1 x2 y(x) = y0 + tan(θ )x − g 2 2 2 v0 cos (θ ) Wir beachten nun die Beziehung 1 cos2 (θ ) = 1 + tan2 (θ ) und erhalten 0 = y(x0 ) = y0 + tan(θ )x0 − gx20 gx20 2 − tan (θ ) 2v20 2v20 Durch Umstellen ergibt sich hieraus eine quadratische Gleichung in tan(θ ) gx20 gx2 tan2 (θ ) − tan(θ )x0 = y0 − 20 2 2v0 2v0 oder tan2 (θ ) − 2 v20 x0 2v20 y0 θ ) = tan( −1 gx20 gx20 Mittels einer quadratischen Ergänzung ergibt sich hieraus 2 2 2 2 v0 v20 v0 2v2 y0 tan (θ ) − 2 = + 02 − 1 tan(θ ) + gx0 gx0 gx0 gx0 2 und damit tan(θ ) = v20 gx0 ± s v20 gx0 2 + 2v20 y0 −1 gx20 1.2. MATHEMATISCHE GRUNDLAGEN 7 woraus sich mithilfe der Umkehrfunktion des Tangens der gesuchte Winkel θ ergibt - sofern man fordert, dass der Winkel zwischen 0 und 45◦ liegen soll. Man nimmt hierzu die negative Wurzel. In dem obigen Beispiel konnten wir den notwendigen Winkel θ analytisch bestimmen - dies ist sicherlich die Ausnahme. Man kann aber das obige Beispiel verallgemeinern. Sei hierzu ein Vektor yZiel gegeben, dann ist ein xziel und cZiel gesucht mit h(xZiel , cZiel ) =k yZiel − Y(xZiel , u, v, cZiel ) k2 = min Das Funktional h(x, c) =k yZiel − Y(x, u, v, c) k2 soll somit minimiert werden. Unterstellen wir nun, dass - das ballistische Problem auf einer kompakten Teilmenge D von G × C betrachtet werden kann - yZiel aus dem Bildraum Y(x, u, v, c)(D) entstammt dann muss es nach dem Satz von Weierstrass vom Minimum und Maximum ein xZiel und ein cZiel geben, welches das obige Funktional minimiert - das obige Funktional h(c) ist ja nach dem vorhergehenden Satz stetig. Das wesentliche Problem besteht damit darin, wie die Größe cZiel berechnet werden kann. Dieses Optimierungsproblem kann aber mit Hilfe eines Verfahrens wie dem BFGS-Verfahren oder dem DFP-Verfahren gelöst werden - siehe hierzu die Funktion fminunc aus der Programmersprache Octave. Für die weiteren Betrachtungen wollen wir das vorliegende Problem etwas verallgemeinern hierbei übernimmt der Rohrerhebungswinkel (und gegebenenfalls der Seitenwinkel) die Rolle eines Hilfsparameter (im folgenden Satz der Parametervektor b) und der Formfaktor (und gegebenenfalls weitere Parameter wie der Magnus-Koeffizient, Lift-Koeffizient etc.) die Rolle des eigentlichen Parameters (im folgenden Satz der Parametervektor c). Der folgende Satz stellt dabei auch noch die Differenzierbarkeit der Integralkurven y = Y(x, u, v, c) sicher. Dabei sei c = (c1 , . . . , ck ) ein k-dimensionaler Vektor, b = (b1 , . . . , bl ) ein l-dimensionaler Vektor, C eine offene Menge im c-Raum und B eine offene Menge im b-Raum. KAPITEL 1. ÜBERBLICK 8 Ist G(x, y) ein Gebiet, dann sei G × C × B die Menge aller Punkte x, y, c, b mit x, y ∈ G, c ∈ C, b ∈ B. Es gilt dann Satz 3 Die Abbildung f(x, y; c, b) sei in G × C × B stetig und besitze stetige partielle Ableitungen bzgl y1 , . . . , yn , c1 , . . . , ck . Für jeden Punkt c, b von C × B gibt es dann durch den Punkt u, v ∈ G eine in G von Rand zu Rand verlaufende Integralkurve y = Y(x, u, v, c, b) der Differentialgleichung y′ = f(x, y; c, b) Die Integralkurve y = Y(x, u, v, c, b) ist an jeder Stelle x, u, v, c, b ihres Existenzbereiches nach x, u, v1 , . . . , vn , c1 , . . . , ck partiell differenzierbar, und diese partiellen Ableitungen sind stetige Funktionen von x, u, v, c, b. Ferner ist die partielle Ableitung wp = ∂Y , ∂ cp 1≤ p≤k die Lösung der linearen Differentialgleichung w′ = fc p (x, Y(x, u, v, c, b)) + F(x, Y(x, u, v, c, b), c, b)w mit w p (u) = 0 Hierbei ist F die Matrix F(x, y, c, b) = ∂ fp , ∂ yq p, q = 1, . . . , n B EWEIS Siehe Kamke (1969), S. 122, Satz 2. Anmerkung Ist die Abbildung f(x, y, c, b) entsprechend oft partiell differenzierbar, dann überträgt sich dies auf die Integralkurven y = Y(x, u, v, c, b) siehe hierzu Kamke (1969), S. 124, Satz 3, und S. 125, Satz 4. 1.2. MATHEMATISCHE GRUNDLAGEN Die angeführten Sätze bilden nun die Arbeitspferde für alle weiteren Betrachtungen. 9 10 KAPITEL 1. ÜBERBLICK Kapitel 2 Diverse Algorithmen 2.1 Ein einführendes Beispiel - das Kaliber .308 Wir betrachten das Kaliber .308 und hier das Produkt D46 vom Hersteller Lapua mit den folgenden ballistischen Daten (zuerst die Geschwindigkeit, dann die kinetische Energie): 0m 100m 200m 300m 600m 800 760 704 650 597 452 369 3466 2974 2533 2141 1228 815 Das Geschossgewicht beträgt hier 185 Grains oder 12 Gramm. Für dieses Produkt ergeben sich folgende Flughöhen (nach Herstellerangaben): 100m 200m 300m 600m 800m 100m 0 −157 −547 −3588 −7908 300m 182 207 0 −2495 −6450 600m 598 1039 1247 0 −3123 Als Luftwiderstandsgesetz unterstellen wir G7 - siehe zum Beispiel Peelen (2014) - und gesucht sei der zugehörige Formfaktor. Das zugehörige Differentialgleichungssysstem lautet wie folgt mit der Abkürzung Ĉd∗ = ρπ Cd 1 · ρ · S ·Cd = , 2m 8 C in Komponentenschreibweise v˙x = −Ĉd∗ vvx v˙y = −Ĉd∗ vvy − g v˙z = −Ĉd∗ vvz 11 C= m d2 KAPITEL 2. DIVERSE ALGORITHMEN 12 Abstrahieren wir von der z-Komponente und separieren wir Ĉd∗ in ein Produkt aus Formfaktor und Luftwiderstandsbeiwert Cd , dann ergeben sich mindestens folgende Möglichkeiten: - es werden Zeitpunkte t1 , . . . , tn vorgegeben (diese Zeitpunkte müssen dann sinnvoll gewählt werden) und zu diesen Zeitpunkte ti werden die zugehörigen Flugweiten x(ti ) und die zugehörigen Flughöhen y(ti ) gemessen. Der Formfaktor wird dann so bestimmt, dass an den vorgegebenen Zeitpunkten ti die Flugweiten x(ti ) und die Flughöhen y(ti ) möglichst gut übereinstimmen. - es werden Flugweiten x1 , . . . , xn vorgegeben und zu diesen Flugweiten werden Flughöhen y1 , . . . , yn gemessen. Der Formfaktor wird dann so bestimmt, dass die Lösungskomponente für die Flughöhe möglichst gut an die gemessenen Flughöhen angepasst wird. Da die Komponente für die Flugweite eine streng monotone und differenzierbare Funktion ist, muss diese Komponente umkehrbar sein, d. h. zu vorgegebenen Flugweiten xi müssen sich die zugehörigen Zeitpunkte ti berechnen lassen. Demzufolge können die Flughöhen als Funktion der Flugweite prinzipiell dargestellt werden. - wieder werden sinnvoll gewählte Zeitpunkte t1 , . . . , tn vorgegeben und die zugehörigen Geschwindigkeiten vi gemessen. Der Formfaktor wird dann so bestimmt, dass die gemessenen Geschwindigkeiten möglichst gut an die durch die Lösung des zugehörigen Differentialgleichungssystems bestimmten Beträge des Geschwindigkeitsvektors angepasst sind. - es werden Flugweiten x1 , . . . , xn vorgegeben und die zugehörigen Geschwindigkeiten v1 , . . . , vn gemessen. Der Formfaktor wird dann so bestimmt, dass die gemessenen Geschwindigkeiten möglichst gut an die durch die Lösung des zugehörigen Differentialgleichungssystems bestimmten Beträge des Geschwindigkeitsvektors angepasst sind. Auch hier gilt: Da die Komponente für die Flugweite eine streng monotone und differenzierbare Funktion ist, muss diese Komponente umkehrbar sein, d. h. zu vorgegebenen Flugweiten xi müssen sich die zugehörigen Zeitpunkte ti berechnen lassen. Demzufolge können die Geschwindigkeiten als Funktion der Flugweite prinzipiell dargestellt werden. - Da die kinetischen Energien Funktionen der Geschwindigkeit sind, ist auch der Weg über die kinetischen Energien denkbar. Wir erhalten damit folgendes Differentialgleichungssystem 2.1. EIN EINFÜHRENDES BEISPIEL - DAS KALIBER .308 13 1 · ρ · S · i ·Cd vvx 2m 1 · ρ · S · i ·Cd vvy − g = − 2m v˙x = − v˙y wobei i der zu bestimmende Formfaktor ist (siehe auch das Buch von Peelen (2014)), ρ die Luftdichte und S die Referenzfläche. v Cd ist der von der Mach-Zahl abhängige Luftwiderstandsbeiwert, demzufolge Cd = Cd ( ), wenn a a die Schallgeschwindigkeit ist. Wir überführen nun das Differentialgleichungssystem zweiter Ordnung in ein System erster Ordnung, indem wir definieren: x y ′ y= vx ⇒ y = vy vx vy v′x v′y und damit ′ vx vy x y = 1 y′ = vx − · ρ · S · i ·Cd vvx 2m 1 vy − · ρ · S · i ·Cd vvy − g 2m wobei v = q v2x + v2y , m die Geschossmasse und ρ die Luftdichte ist. Nehmen wir zuerst an, wir geben eine Fleckschussentfernung von zum Beispiel x(t100 ) = 100 Metern vor, und die Flughöhen sind gegeben. Ebenso haben wir Anfangsbedingungen x(0) = 0, y(0) = y0 , vx (0) = v0 cos(θ ) und vy (0) = v0 sin(θ ), wobei der Rohrerhebungswinkel θ als ein Parameter fungiert. Die Lösung des Differentialgleichungssystems besteht dann aus zwei Funktionen x(t) und y(t) wie auch Geschwindigkeiten vx (t) und vy (t). Diese Funktionen hängen damit von dem Formfaktor i ab, d. h. x(t) y(t) vx (t) vy (t) = = = = Gegeben seien nun folgende Daten vom Anfang: x(t; i) y(t; i) vx (t; i) vy (t; i) KAPITEL 2. DIVERSE ALGORITHMEN 14 100m 200m 300m 600m 800m 100m 0 −157 −547 −3588 −7908 300m 182 207 0 −2495 −6450 600m 598 1039 1247 0 −3123 Wollen wir somit die Flughöhen möglichst genau vorhersagen, müssen wir den Formfaktor so bestimmen, dass das folgende Funktional minimal wird: h(i) = (0 − y(t100 ; i))2 + (−0, 157 − y(t200 ; i))2 + (0, 547 − y(t300 ; i))2 +(−3, 588 − y(t600 ); i))2 + (−7, 908 − y(t800 ); i))2 Dabei muss zuerst der Rohrerhebungswinkel für die Fleckschussweite bestimmt und danach das angegebene Funktional minimiert werden. Eine geeignete Funktion in Octave ist die Funktion fminbnd - hier liegt ja eine eindimensionale Funktion vor, sie basiert auf dem Verfahren des goldenen Schnitts. Aufgrund der vorhergehenden Sätze muss eine optimale Lösung für den Formfaktor existieren es ist dann sicher interessant, einen Vergleich anzustellen, welche Ergebnisse mit der Angabe für den Formfaktor aus der Literatur erhalten werden und welche Ergebnisse mit dem optimierten Formfaktor erzielt werden. Zu den oben angegebenen Daten ergibt sich durch Rechnung mit dem Formfaktor 1, 081 aus dem Buch von Litz (2011), S. 530, und dem Luftwiderstandsgesetz G7: 100m 200m 300m 600m 800m 100m 0 −158 −549 −3584 −7863 300m 183 208 0 −2486 −6400 600m 597 1068 1243 0 −3085 Die erste Spalte beinhaltet die Fleckschussweiten. Berechnen wir gemäß der vorhergehenden mathematischen Sätze den optimalen Formfaktor zu 1, 08944, so erhalten wir die folgende Tabelle: 100m 200m 300m 600m 800m 100m 0 −158 −550 −3597 −7905 300m 183 208 0 −2498 −6439 600m 600 1041 1249 0 −3109 Hier nochmals zum Vergleich die Herstellerdaten: 2.2. DAS KALIBER .300 WIN MAG 15 100m 200m 300m 600m 800m 100m 0 −157 −547 −3588 −7908 300m 182 207 0 −2495 −6450 600m 598 1039 1247 0 −3123 Wir betrachten nun die Daten mit dem Produkt-Code GB491 von Lapua (Kaliber .308). Hier haben wir folgende ballistische Daten (zuerst die Geschwindigkeit, dann die kinetische Energie): 0m 100m 200m 300m 600m 800 1000m 860 794 730 669 503 408 331 3698 3151 2663 2235 1264 833 548 Das Geschossgewicht beträgt hier 155 Grains oder 10 Gramm. Für dieses Produkt ergeben sich folgende Flughöhen (nach Herstellerangaben, die erste Spalte beinhaltet wieder die Fleckschussweiten): 100m 200m 300m 600m 800m 1000m 100m 0 −115 −414 −2813 −6265 −12106 300m 138 161 0 −1984 −5160 −10726 600m 469 822 992 0 −2515 −7418 und berechnet mit dem Formfaktor 0, 988 aus dem Buch von Litz (2011), S. 527: 100m 200m 300m 600m 800m 1000m 100m 0 −116 −414 −2795 −6206 −11984 300m 138 161 0 −1967 −5101 −10603 600m 466 816 983 0 −2479 −7325 Optimieren wir jetzt gemäß unseren mathematischen Sätzen den Formfaktor und rechnen damit das Differentialgleichungssystem durch, dann erhalten wir folgende Tabelle: 100m 300m 600m 100m 200m 300m 600m 800m 1000m 0 −116 −416 −2811 −6253 −12111 139 161 0 −1979 −5145 −10725 468 821 990 0 −2506 −7427 Der aus den Daten geschätzte Formfaktor ist hierbei 0, 99852. Die Berechnungen wurden mit dem Laufdurchmesser 7, 83 mm durchgeführt (siehe hierzu das Buch von Peelen (2014) auf S. 19, wo eine Tabelle mit den Laufdurchmessern angegeben ist). 2.2 Das Kaliber .300 Win Mag Als weiteres Beispiel betrachten wir das Kaliber .300 Win Mag. KAPITEL 2. DIVERSE ALGORITHMEN 16 RWS gibt in der Serie Target Elite Plus für dieses Kaliber folgende Daten an (Geschossgewicht 13 Gramm, ballistischer Koeffizient 0.555): 0m 200m 400m 600m 800 1000m 870 762 661 568 484 407 4904 3758 2833 2091 1521 1073 Die zweite Zeile umfasst die Geschwindigkeiten, die dritte Zeile die kinetische Energie. Für die Trajektorien gibt RWS folgende Daten an (in Millimetern): 100m 300m 500m 200m 400m 600m 800m 1000m −98 −827 −2446 −5345 −9924 144 −343 −1721 −4377 −8715 507 382 −633 −2927 −6903 Leider macht RWS keine Angaben (?) über die Höhe des Zielfernrohrs über der Laufachse, es steht aber zu vermuten, dass es 4 oder 5 Zentimeter sein könnten (dies ist die übliche Höhe). Bei 5 Zentimetern ergibt sich mit dem vorgestellten Verfahren und Luftwiderstandsgesetz G7 ein Formfaktor von 1, 07755 - rechnet man damit das Differentialgleichungssystem durch, so erhält man folgende Tabelle: 100m 300m 500m 200m 400m 600m 800m 1000m −98 −827 −2463 −5332 −9927 146 −340 −1732 −4358 −8710 508 384 −646 −2909 −6899 Bei einer Höhe des Zielfernrohrs von 4 Zentimetern werden die Ergebnisse schlechter - von daher wurden die Ergebnisse von RWS wohl mit einer Höhe von 5 Zentimetern gewonnen. Der angegebenen ballistische Koeffizient scheint für das Luftwiderstandsgesetz G1 zu gelten - rechnet man das Differentialgleichungssystem mit G1 durch, so erhält man Ergebnisse, die leidlich zu den angegebenen Daten passen. Allerdings - so Litz (2011), S. 13 ff. - ist das Luftwiderstandsgesetz G7 geeigneter für Büchsenpatronen mit entsprechenden Geschossformen. Das Luftwiderstandsgesetz G1 ist eher geeignet für Revolver- und Pistolengeschosse aufgrund ihrer eher zylindrischen Form - dafür sind die ballistischen Koeffizienten größer (?) . Anmerkung Die analoge Rechnung wurde bei gleicher Zielfernrohrhöhe auch für das Kaliber .308 durchgeführt (Geschossgewicht 10,69 Gramm, Mündungsgeschwindigkeit 805 Meter pro Sekunde, ebenfalls Target Elite plus von RWS). Hier der Vergleich - zuerst die Herstellerangaben: 2.3. DAS KALIBER .338 LAPUA MAGNUM 17 0m 200m 400m 600m 800 1000m 805 676 555 447 362 311 3529 2487 1675 1088 715 528 für Geschwindigkeit und Energie, nun die Flughöhen 200m 400m 600m 800m 1000m 100m −131 −1078 −3336 −7644 −15035 300m 186 −443 −2385 −6375 −13449 500m 677 538 −913 −4413 −10996 Der optimale Formfaktor ergibt sich zu 1, 12667 - mit diesem Formfaktor das Differentialgleichungssystem mit dem Luftwiderstandsgesetz G7 durchgerechnet ergibt die Trajektorien: 200m 400m 600m 800m 1000m 100m −131 −1087 −3342 −7608 −15048 300m 188 −450 −2386 −6334 −13456 500m 676 526 −922 −4382 −11015 2.3 Das Kaliber .338 Lapua Magnum Die Munition mit dem Produkt-Code GB528 von Lapua hat folgende ballistische Daten (zuerst die Geschwindigkeit, dann die kinetische Energie): 0m 100m 200m 300m 600m 800 1000m 830 789 750 712 605 539 478 6696 6054 5465 4925 3558 2827 2223 Das Geschossgewicht beträgt hier 300 Grains oder 19,4 Gramm. Für dieses Produkt ergeben sich folgende Flughöhen (nach Herstellerangaben): 100m 300m 600m 100m 200m 300m 600m 800m 1000m 0 −117 −413 −2583 −5231 −9289 136 155 0 −1710 −4146 −7932 421 725 855 0 −1865 −5082 100m 300m 600m 100m 200m 300m 600m 800m 1000m 0 −117 −406 −2506 −5178 −9169 135 154 0 −1694 −4094 −7814 418 719 847 0 −1836 −4992 und berechnet: KAPITEL 2. DIVERSE ALGORITHMEN 18 Auch hier beinhaltet die erste Spalte die Fleckschussweiten. Der Formfaktor wurde dem Buch von Litz (2011), Seite 557, entnommen zu 0, 956. Optimieren wir nun den Formfaktor, so erhalten wir hierfür den Wert 0, 98351, und es ergibt sich folgende Tabelle nach dem Lösen des Differentialgleichungssystems: 100m 200m 300m 600m 800m 1000m 100m 0 −117 −408 −2526 −5233 −9297 300m 136 155 0 −1710 −4145 −7936 600m 421 725 855 0 −1865 −5086 Hier wurde mit dem Laufdurchmesser 8, 59 mm gearbeitet. 2.4 Schätzung des Formfaktors mithilfe der Geschwindigkeit Wir betrachten nochmals das Kaliber .308 und hier das Produkt D46 vom Hersteller Lapua mit den folgenden ballistischen Daten (zuerst die Geschwindigkeit, dann die kinetische Energie): 0m 100m 200m 300m 600m 800 760 704 650 597 452 369 3466 2974 2533 2141 1228 815 Das Geschossgewicht beträgt hier 185 Grains oder 12 Gramm. Falls es schwierig, wenn nicht unmöglich ist, die Flughöhen zu messen, liegt es nahe, stattdessen die Geschwindigkeiten für die Schätzung des gesuchten Formfaktors zu verwenden. Die Waffe sei hierzu auf 100 Meter eingeschossen. Wieder betrachten wir das Differentialgleichungssystem ′ vx x vy y 1 y′ = vx = − · ρ · S · i ·Cd vvx 2m 1 vy − · ρ · S · i ·Cd vvy − g 2m wobei v = q v2x + v2y , m die Geschossmasse und ρ die Luftdichte ist. Die Lösungen seien wie vorhin x(t) y(t) vx (t) vy (t) = = = = x(t; i) y(t; i) vx (t; i) vy (t; i) 2.4. SCHÄTZUNG DES FORMFAKTORS MITHILFE DER GESCHWINDIGKEIT 19 sie hängen somit von dem Formfaktor i ab. Wir bilden nun die Geschwindigkeit v = v(t) = tional p vx (t; i)2 + vy (t; i)2 und damit das folgende Funk- h(i) = (704 − v(t100 ; i))2 + (650 − v(t200 ; i))2 + (597 − v(t300 ; i))2 +(452 − v(t600 ); i))2 + (369 − v(t800 ); i))2 Dabei werden die Zeitpunkte t j , j = 100, j = 200, j = 300, j = 600 und j = 800, so bestimmt, dass x(t j ; i) = j. Dieses Funktional muss nun bzgl. des Formfaktors i minimiert werden - es liegt nahe, wie in den vorhergehenden Beispielen vorzugehen und die Funktion fminbnd aus dem Programmpaket Octave zu verwenden. Wir erhalten dann wiederum den optimalen Formfaktor, mit dem dann die Flugbahn berechnet werden kann - hier ergibt sich für den Formfaktor der Wert 1, 09748. Rechnet man hierfür das Differentialgleichungssystem durch, dann erhält man folgende Tabelle: 100m 200m 300m 600m 800m 100m 0 −159 −551 −3610 −7945 300m 184 209 0 −2508 −6476 600m 602 1045 1254 0 −3132 Hier nochmals zum Vergleich die Herstellerdaten: 100m 200m 300m 600m 800m 100m 0 −157 −547 −3588 −7908 300m 182 207 0 −2495 −6450 600m 598 1039 1247 0 −3123 Betrachten wir nun noch die zweite Laborierung des Kalibers .308 - hier hatten wir folgende Angaben für die Geschwindigkeit und die kinetische Energie: 0m 100m 200m 300m 600m 800 1000m 860 794 730 669 503 408 331 3698 3151 2663 2235 1264 833 548 Für dieses Produkt ergeben sich folgende Flughöhen (nach Herstellerangaben, die erste Spalte beinhaltet wieder die Fleckschussweiten): KAPITEL 2. DIVERSE ALGORITHMEN 20 100m 200m 300m 600m 800m 1000m 100m 0 −115 −414 −2813 −6265 −12106 300m 138 161 0 −1984 −5160 −10726 600m 469 822 992 0 −2515 −7418 und berechnet mit dem Formfaktor 0, 988 aus dem Buch von Litz (2011), S. 527: 100m 200m 300m 600m 800m 1000m 100m 0 −116 −414 −2795 −6206 −11984 300m 138 161 0 −1967 −5101 −10603 600m 466 816 983 0 −2479 −7325 Optimieren wir jetzt gemäß unseren mathematischen Sätzen den Formfaktor und rechnen damit das Differentialgleichungssystem durch, dann erhalten wir folgende Tabelle: 100m 300m 600m 100m 200m 300m 600m 800m 1000m 0 −116 −416 −2811 −6253 −12111 139 161 0 −1979 −5145 −10725 468 821 990 0 −2506 −7427 Optimieren wir nun den Formfaktor mittels der Geschwindigkeit, so erhalten wir für den Formfaktor den Wert 0, 99847 und damit folgende Tabelle 100m 300m 600m 100m 200m 300m 600m 800m 1000m 0 −116 −416 −2810 −6253 −12110 139 161 0 −1979 −5145 −10725 468 821 990 0 −2506 −7426 Es verbleibt jetzt noch, die Daten für das Kaliber .338 Lapua Magnum durchzurechnen. Hier hatten wir folgende Herstellerdaten (zuerst die Geschwindigkeit, dann die kinetische Energie): 0m 100m 200m 300m 600m 800 1000m 830 789 750 712 605 539 478 6696 6054 5465 4925 3558 2827 2223 Das Geschossgewicht beträgt hier 300 Grains oder 19,4 Gramm. Für dieses Produkt ergeben sich folgende Flughöhen (nach Herstellerangaben): 100m 300m 600m 100m 200m 300m 600m 800m 1000m 0 −117 −413 −2583 −5231 −9289 136 155 0 −1710 −4146 −7932 421 725 855 0 −1865 −5082 2.5. SCHÄTZUNG DES FORMFAKTORS MITHILFE DER ENERGIE 21 Nach dem Optimieren erhielten wir seinerzeit einen Formfaktor mit dem Wert 0, 98351, und folgende Tabelle nach dem Lösen des Differentialgleichungssystems: 100m 200m 300m 600m 800m 1000m 100m 0 −117 −408 −2526 −5233 −9297 300m 136 155 0 −1710 −4145 −7936 600m 421 725 855 0 −1865 −5086 Optimieren wir den Formfaktor über die Geschwindigkeit, so erhalten wir für den Formfaktor den Wert 0, 98521, und folgende Tabelle nach dem Lösen des Differentialgleichungssystems: 100m 300m 600m 100m 200m 300m 600m 800m 1000m 0 −117 −408 −2528 −5237 −9305 136 155 0 −1711 −4148 −7944 421 725 856 0 −1867 −5092 2.5 Schätzung des Formfaktors mithilfe der Energie Die kinetische Energie ist gegeben durch 1 · m · v2 2 von daher besteht die Möglichkeit, die kinetische Energie - zu vorgegebenen Zeitpunkten t1 , . . . , tn zu messen und ein entsprechendes Funktional bzgl. des Formfaktors zu minimieren - zu vorgegebenen Flugweiten x1 , . . . , xn zu messen und ebenfalls ein entsprechendes Funktional bzgl. des Formfaktors zu minimieren Da die kinetische Energie eine Funktion der Geschwindigkeit ist, steht allerdings zu erwarten, dass sich bei dieser Vorgehensweise keine neuen Erkenntnisse gewinnen lassen. Aus praktischer Sicht ist es eine Frage der Bequemlichkeit, ob man lieber mit der Geschwindigkeit oder der kinetischen Energie arbeiten möchte. 2.6 Schätzung des Formfaktors mithilfe der Zeiten Wir betrachten wieder das Kaliber .308 und hier das Produkt D46 vom Hersteller Lapua mit den folgenden ballistischen Daten (zuerst die Geschwindigkeit, dann die kinetische Energie): 0m 100m 200m 300m 600m 800 760 704 650 597 452 369 3466 2974 2533 2141 1228 815 KAPITEL 2. DIVERSE ALGORITHMEN 22 Das Geschossgewicht beträgt hier 185 Grains oder 12 Gramm. Nun gehen wir davon aus, dass die Zeiten gemessen werden, die zum Erreichen der Entfernungen 100, 200, 300, 600 und 800 Meter benötigt werden - diese seien mit t˜100 , t˜200 , t˜300 , t˜600 und t˜800 bezeichnet. Wiederum soll der Formfaktor i optimal geschätzt werden. Die Waffe sei hierzu auf 100 Meter eingeschossen. Wieder betrachten wir das Differentialgleichungssystem ′ vx vy x y = 1 y′ = vx − · ρ · S · i ·Cd vvx 2m 1 vy − · ρ · S · i ·Cd vvy − g 2m wobei v = q v2x + v2y , m die Geschossmasse und ρ die Luftdichte ist. Die Lösungen seien wie vorhin x(t) y(t) vx (t) vy (t) = = = = x(t; i) y(t; i) vx (t; i) vy (t; i) sie hängen somit von dem Formfaktor i ab. Wir berechnen nun zu den Entfernungen 100, 200, 300, 600 und 800 die Zeiten aus den Gleichungen 100 = x(t100 ; i), 200 = x(t200 ; i), 300 = x(t300 ; i), 600 = x(t600 ; i) und 800 = x(t800 ; i) - die so gewonnenen Zeitpunkte t j hängen dann natürlich noch von dem Formfaktor i ab. Wir betrachten dann folgendes Funktional h(i) = (t˜100 − t100 )2 + (t˜200 − t200 )2 + (t˜300 − t300 )2 +(t˜600 − t600 )2 + (t˜800 − t800 )2 Dieses Funktional muss dann bzgl. i minimiert werden - dies ergibt dann die optimale Schätzung für den gesuchten Formfaktor. Allerdings sind in dieser Situation keine Daten verfügbar - Auswertungen sind somit nicht möglich. 2.7. VERALLGEMEINERUNGEN 23 2.7 Verallgemeinerungen Der Ausgangspunkt des einführenden Beispiels war das Differentialgleichungssystem 1 · ρ · S · i ·Cd vvx 2m 1 = − · ρ · S · i ·Cd vvy − g 2m v˙x = − v˙y wobei i der zu bestimmende Formfaktor ist (siehe auch das Buch von Peelen (2014)). Die erste Erweiterung besteht darin, anstatt des zweidimensionalen Systems ein dreidimensionales System zu betrachten. 1 · ρ · S · i ·Cd vvx 2m 1 · ρ · S · i ·Cd vvy − g = − 2m 1 = − · ρ · S · i ·Cd vvz 2m v˙x = − v˙y v˙z Dieses System zweiter Ordnung kann man natürlich analog in ein System erster Ordnung überführen, indem wir definieren: y= x y z vx vy vz ⇒ y′ = vx vy vz v′x v′y v′z und damit ′ y = wobei v = q x y z vx vy vz ′ vx vy vz 1 = · ρ · S · i ·Cd vvx − 2m − 1 · ρ · S · i ·Cd vvy − g 2m 1 − · ρ · S · i ·Cd vvz 2m v2x + v2y , m die Geschossmasse und ρ die Luftdichte ist. In dieser Formulierung ist es jetzt möglich, auch den Einfluss des Seitenwindes zu modellieren. KAPITEL 2. DIVERSE ALGORITHMEN 24 Dieses Differentialgleichungssystem kann man nun um Rotationseffekte erweitern - siehe hierzu McCoy (1999) oder Carlucci und Jacobson (2008). Diese Erweiterung beruht auf dem Drehimpulssatz - hierzu müssen wir folgende Größen definieren - siehe McCoy (1999), S. 189: Cd CL CNp CMp CMq Cl p CMα = = = = = = = Luftwiderstandskoeffizient Auftriebskoeffizient Koeffizient für die Magnus-Kraft Koeffizient für das Magnus-Moment Koeffizient für die Kraft der Pitch-Dämpfung Koeffizient der Spin-Dämpfung Koeffizient des Überschlagsmoments Um in den Differentialgleichungen nicht ständig die Masse wie auch die Spin-Rate mitschleppen zu müssen, bildet man folgende abgekürzte Größen: C̃d = C̃L = C̃Np = C̃Mp = C̃Mq = C̃l p = C̃Mα = ρ vSCd 2m ρ vSCL 2m ρ vSCNp p 2m ρ Sd 2CMp p 2Iy ρ vSdCMq 2m ρ vSd 2Cl p p 2Iy ρ vSdCMα 2Iy wobei p= Iy (h · x) Ix die axiale Spin-Rate des Projektils v1 = V1 −W1 , v2 = V2 −W2 , v3 = V3 −W3 , v = cos(α ) = q v21 + v22 + v23 v1 x1 + v2 x2 + v3 x3 v der Kosinus des Anstellwinkels des Projektils ist, Wi , i = 1, 2, 3 die Komponenten der Windgeschwindigkeit und ẋi = vi . 2.7. VERALLGEMEINERUNGEN 25 Vi sind die Komponenten des Geschwindigkeitsvektors des Projektils bzgl. eines erdfesten Koordinatensystem, vi die Komponenten des Geschwindigkeitsvektors des Projektils bzgl der Luft. Ix und Iy sind die Trägheitsmomente des Projektils in x- bzw. y-Richtung des Projektils. Hiermit können wir dann folgende sechs gewöhnliche Differentialgleichungen aufstellen: v̇1 = −C̃d v1 + C̃L v2 x1 − vv1 cos(α ) + C̃Np (x2 v3 − x3 v2 ) + C̃Nq (h2 x3 − h3 x2 ) + g1 v̇2 = −C̃d v2 + C̃L v2 x2 − vv2 cos(α ) + C̃Np (x3 v1 − x1 v3 ) + C̃Nq (h3 x1 − h1 x3 ) + g2 v̇3 = −C̃d v3 + C̃L v2 x3 − vv3 cos(α ) + C̃Np (x1 v2 − x2 v1 ) + C̃Nq (h1 x2 − h2 x1 ) + g3 ḣ1 = C̃l p x1 + C̃Mα (v2 x3 − v3 x2 ) + C̃Np (v1 − vx1 cos(α )) + C̃Mq ḣ2 = C̃l p x2 + C̃Mα (v3 x1 − v1 x3 ) + C̃Mp (v2 − vx2 cos(α )) + C̃Mq Ix p x2 h2 − Iy ḣ3 = C̃l p x3 + C̃Mα (v1 x2 − v2 x1 ) + C̃Mp (v3 − vx3 cos(α )) + C̃Mq Ix p x3 h3 − Iy Ix p h1 − x1 Iy wobei g1 , g2 und g3 die Komponenten des Gravitationsvektors sind. Dabei ist h = (h1 , h2 , h3 )T und H der Drehimpulsvektor - h ist der Drehimpulsvektor H durch Iy devidiert. Beachtet man noch, dass v̇i = V̇i , dann nehmen die sechs Differentialgleichungen die folgende Form an: V̇1 = −C̃d v1 + C̃L v2 x1 − vv1 cos(α ) + C̃Np (x2 v3 − x3 v2 ) + C̃Nq (h2 x3 − h3 x2 ) + g1 V̇2 = −C̃d v2 + C̃L v2 x2 − vv2 cos(α ) + C̃Np (x3 v1 − x1 v3 ) + C̃Nq (h3 x1 − h1 x3 ) + g2 V̇3 = −C̃d v3 + C̃L v2 x3 − vv3 cos(α ) + C̃Np (x1 v2 − x2 v1 ) + C̃Nq (h1 x2 − h2 x1 ) + g3 ḣ1 = C̃l p x1 + C̃Mα (v2 x3 − v3 x2 ) + C̃Np (v1 − vx1 cos(α )) + C̃Mq ḣ2 = C̃l p x2 + C̃Mα (v3 x1 − v1 x3 ) + C̃Mp (v2 − vx2 cos(α )) + C̃Mq Ix p x1 h1 − Iy Ix p x2 h2 − Iy KAPITEL 2. DIVERSE ALGORITHMEN 26 ḣ3 = C̃l p x3 + C̃Mα (v1 x2 − v2 x1 ) + C̃Mp (v3 − vx3 cos(α )) + C̃Mq Ix p x3 h3 − Iy Dieses System zweiter Ordnung kann man aber leicht umschreiben in ein System erster Ordnung. Dazu definieren wir (plus geeignete Anfangsbedingungen) y= X1 X2 X3 h1 h2 h3 ⇒ ẏ = Ẋ1 Ẋ2 Ẋ3 ḣ1 ḣ2 ḣ3 = ẏ = V1 V2 V3 ḣ1 ḣ2 ḣ3 Fasst man nun die rechten Seiten der sechs gewöhnlichen Differentialgleichungen zu einer Funktion f(y) zusammen, so erhält man offensichtlich ein Differentialgleichungssystem ẏ = f(y) dem noch geeignete Anfangsbedingungen hinzuzufügen sind. Dies bedeutet, dass eine Fleckschussentfernung gewählt worden sein muss, auf die sich die Anfangsbedingungen beziehen. Bei dem Differentialgleichungssystem ist nun f(y) = f(y; C̃d , C̃L , C̃Np , C̃Mp , C̃Mq , C̃l p , C̃Mα ), wobei C̃d eine Funktion des Formfaktors i ist. Diese sieben Parameter C̃d , C̃L , C̃Np , C̃Mp , C̃Mq , C̃l p , C̃Mα müssen somit aus geeigneten Daten geschätzt werden. Setzen wir nun voraus, dass geeignete Daten vorliegen und dass y(t; C̃d , C̃L , C̃Np , C̃Mp , C̃Mq , C̃l p , C̃Mα ) die Lösung des Differentialgleichungssystems darstellt. Um die Parameter zu schätzen, muss ein geeignetes Funktional definiert werden - analog zu den vorhergehenden Beispielen. 2.7.1 Schätzung der Parameter aus Flughöhe und Seitenabweichung Nehmen wir an, dass die zweite Komponente im Vektor y die Flughöhe und die dritte Komponente im Vektor y die seitliche Abweichung darstellen. 2.7. VERALLGEMEINERUNGEN 27 Sind nun y(ti ) für i = 1, . . . , n die Daten für die Flughöhen und z(ti ) für i = 1, . . . , n die Daten für die seitlichen Abweichungen bei bestimmten Flugweiten x(ti ), dann liegt es nahe, folgendes Funktional zu minimieren: n n i=1 i=1 h(C̃d , C̃L , C̃Np , C̃Mp , C̃Mq , C̃l p , C̃Mα )) = ∑ (y(ti ) − X2 (ti ))2 + ∑ (z(ti ) − X3(ti ))2 wobei X2 (t) und X3 (t) jeweils von C̃d , C̃L , C̃Np , C̃Mp , C̃Mq , C̃l p , C̃Mα abhängen. Dabei sind die ti aus den zugehörigen Flugweiten zu berechnen. Die Minimierung dieses Funktionals entspricht prinzipiell der Lösung eines nichtlinearen Ausgleichsproblems - die Funktion fminunc aus dem Programm Octave oder Matlab löst diese Aufgabe. 2.7.2 Schätzung der Parameter aus der Geschwindigkeit Wir kürzen den Parametervektor durch C̃ = (C̃d , C̃L , C̃Np , C̃Mp , C̃Mq , C̃l p , C̃Mα )T ab. Analog zu der Vorgehensweise in den vorausgegangenen Beispielen sei v(t j ; C̃) = q vx (t j ; C̃)2 + vy (t j ; C̃)2 + vz (t j ; C̃)2 der Betrag des Geschwindigkeitsvektors zur Zeit t j , j = 1, . . . , n. Da die Geschwindigkeit in gewissen Entfernungen j = x(t j ; i) gemessen wird, ergeben sich hieraus die zugehörigen Beträge der Geschwindigkeitsvektoren. Seien nun in bestimmten Entfernungen o. B. d. A. z. B. j = 100, j = 200, j = 300, j = 400, j = 500, j = 600, j = 800 und j = 1000 Meter die zugehörigen Geschwindigkeiten v100 , v200 , v300 , v400 , v500 , v600 , v800 und v1000 gemessen. Wir betrachten dann das Funktional h(C̃) = (v100 − v(t100 ; C̃))2 + (v200 − v(t200 ; C̃))2 + (v300 − v(t300 ; C̃))2 + (v400 − v(t400 ; C̃))2 +(v500 − v(t500 ; C̃))2 + (v600 − v(t600 ; C̃))2 + (v800 − v(t800 ; C̃))2 + (v1000 − v(t1000 ; C̃))2 KAPITEL 2. DIVERSE ALGORITHMEN 28 oder ganz allgemein n h(C̃) = ∑ (v j − v(t j ; C̃))2 i=1 Dieses Funktional muss nun in Abhängigkeit von C̃ minimiert werden - eine derartige Aufgabe löst die Funktion fminunc aus dem Programm Octave. 2.7.3 Schätzung der Parameter aus der Geschossenergie Analog zu der Vorgehensweise in den vorausgegangenen Beispielen sei v(t j ; C̃) = q vx (t j ; C̃)2 + vy (t j ; C̃)2 + vz (t j ; C̃)2 der Betrag des Geschwindigkeitsvektors zur Zeit t j , j = 1, . . . , n. Hieraus kann die kinetische Energie 1 e(t j ; C̃) = mv(t j ; C̃)2 2 berechnet werden. Seien nun in bestimmten Entfernungen o. B. d. A. z. B. j = 100, j = 200, j = 300, j = 400, j = 500, j = 600, j = 800 und j = 1000 Meter die zugehörigen kinetischen Energien e100 , e200 , e300 , e400 , e500 , e600 , e800 und e1000 gemessen. Wir betrachten dann das Funktional h(C̃) = (e100 − e(t100 ; C̃))2 + (e200 − e(t200 ; C̃))2 + (e300 − e(t300 ; C̃))2 + (e400 − e(t400 ; C̃))2 +(e500 − e(t500 ; C̃))2 + (e600 − e(t600 ; C̃))2 + (e800 − e(t800 ; C̃))2 + (e1000 − e(t1000 ; C̃))2 oder ganz allgemein n h(C̃) = ∑ (e j − e(t j ; C̃))2 i=1 Dieses Funktional muss nun in Abhängigkeit von C̃ minimiert werden - eine derartige Aufgabe löst die Funktion fminunc aus dem Programm Octave. 2.8. BEREITSTELLUNG DER DATEN 29 2.8 Bereitstellung der Daten Bisher wurde vorausgesetzt, dass die Daten zur Schätzung des Formfaktors etc. bereits vorliegen. Normalerweise müssen diese aber erst erzeugt werden - es muss zum Beispiel auf bestimmten Entfernungen erst eine bestimmte Anzahl von Schüssen abgegeben werden. Damit läuft man aber sofort in eine Stichprobenproblematik hinein, d. h. es stellt sich die Frage, wieviele Treffer muss man aufnehmen in welcher Entfernung. Eine Möglichkeit besteht darin, ein Konfidenzintervall vorgegebener Länge bzw. ein Konfidenzellipsoid vorgegebener Länge zu betrachten und hierüber den notwendigen Stichprobenumfang zu berechnen. Diese Problematik ist aber in der Literatur bekannt, siehe zum Beispiel das Buch von Morrison (1976) oder das Buch von Seber (1984). 2.9 Auswahl eines adäquaten Luftwiderstandsmodells Die vorgestellten Verfahren ermöglichen offensichtlich, unter gegebenen Modellen ein optimales Luftwiderstandsmodell auszuwählen. Berechnet man hierzu für optimalen Formfaktor die Standardabweichungen des jeweiligen Modells, so liegt es nahe, dasjenige Modell zu wählen, bei dem die Standardabweichung minimal ist. Sind beispielsweise G1 bis G8 die zur Auswahl stehenden Luftwiderstandsmodelle mit zugehörigen Standardabweichungen s1 bis s8 , dann würde das Modell gewählt werden, deren zugehörige Standardabweichung minimal ist. Eine weitere Möglichkeit besteht darin, das Funktional selbst zur Entscheidung heranzuziehen - je besser die Anpassung durch das Luftwiderstandsmodell, umso kleiner muss das Funktional zum Beispiel h(i) = (0 − y(t100 ; i))2 + (−0, 157 − y(t200 ; i))2 + (0, 547 − y(t300 ; i))2 +(−3, 588 − y(t600 ); i))2 + (−7, 908 − y(t800 ); i))2 (siehe das Beispiel vom Anfang) in Abhängigkeit von dem Luftwiderstandsgesetz sein. Ein derartiges Funktional beschreibt ja gerade eines Teils der Lösung des Differentialgleichungssystems an die Daten - im vorliegenden Fall an die Flughöhen. Betrachten wir hierzu nochmals die Daten des Kalibers .308 der Target Elite Plus von RWS: KAPITEL 2. DIVERSE ALGORITHMEN 30 200m 400m 600m 800m 1000m 100m −131 −1078 −3336 −7644 −15035 300m 186 −443 −2385 −6375 −13449 500m 677 538 −913 −4413 −10996 Mit der Fleckschussweite von 100 Metern erhalten wir für das Luftwiderstandsgesetz G7 einen Formfaktor von 0, 12667 und für das Funktional einen Wert von 0, 0015967 - bei einer Zielfernrohrhöhe von 5 Zentimetern. Bei gleicher Fleckschussweite und gleicher Ziefernrohrhöhe ergibt das Funktional für das Luftwiderstandsgesetz G1 einen Wert von 0, 0038125 und für den Formfaktor einen Wert von 1, 13100. Der Wert des Funktionals ist für das Luftwiderstandsgesetz G1 mehr als doppelt so groß, damit passt das Luftwiderstandsgesetz G7 besser zu den Daten als das Luftwiderstandsgesetz G1, was die Behauptung von Bryan Litz stützt. 2.10 Erweiterungen des Verfahrens Das vorgestellte Verfahren wurde bisher nur für Munition für Handfeuerwaffen konzipiert. Erweiterungen sind vorstellbar: - Wählt man ein geeignetes Luftwiderstandsgesetz (ein Luftwiderstandsgesetz, welches Artilleriegeschosse adäquat beschreibt), dann eröffnet sich die Möglichkeit, die entsprechenden Parameter C̃ = (C̃d , C̃L , C̃Np , C̃Mp , C̃Mq , C̃l p , C̃Mα )T optimal zu schätzen. Damit wäre auch eine potentielle Anwendung für die Artillerie aufgezeigt. - Ebenso liegt es nahe, mit weiteren Luftwiderstandsgesetzen zu arbeiten - das vorgeschlagene Verfahren ist ja unabhängig von den Gesetzen G1 bis G7. Literaturverzeichnis [Carlucci und Jacobson 2008] CARLUCCI, D. E. ; JACOBSON, S. S.: Ballistics Theory and Design of Guns and Ammunition. CRC Press, 2008 [Kamke 1969] [Litz 2011] KAMKE, E.: Differentialgleichungen. Akademische Verlagsgesellschaft, 1969 L ITZ, B.: Applied Ballistics for Long-Range Shooting. Applied Ballistics, LLC, 2011 [McCoy 1999] M C C OY, R. L.: Modern Exterior Ballistics. Schiffer Military History, 1999 [Morrison 1976] pany, 1976 M ORRISON, D. F.: Multivariate Statistical Methods. McGraw-Hill Book Com- [Peelen 2014] P EELEN, J.: Ballistik kleiner Kaliber: der Luftwiderstand. Verlagshaus Monsenstein und Vannerdat OHG Münster, 2014 [Seber 1984] S EBER, G. A. F.: Multivariate Observations. Wiley, 1984 31 32 LITERATURVERZEICHNIS
© Copyright 2025 ExpyDoc