zuletzt geändert: 22.03.98
Die Phygoide: Wellenflug-Simulation
und Wirklichkeit.
Das
Programm
Beispiel
Was lernt man daraus?
Ergebnis der Simulation der Wellenkompensation
Flugparametermessungen an echten Modellflugzeugen
Die Phygoide: ist sie meßbar an unseren Modellen ?
What next ? ..... Was kommt demnächst
?
Die Simulation
Um meinem Traum eines Thermik-Autopiloten ein bißchen näher zu kommen, muß man zuerst die Fluglage eines Flugzeugs gut unter Kontrolle halten. Der allererste Schritt in dieser Richtung ist die Beseitigung des Schaukelns um die Querachse, die Beseitigung oder Verkleinerung der Phygoide.
Helmut Schenk hat sehr schön die Phygoide beschrieben.
Hier möchte ich nicht nur die graphischen Plots der Lösungen
der Gleichungen zeigen, sondern auch die Vorgehensweise darstellen,
die ähnlich in anderen modernen Mathematikprogrammen verwirklicht
werden kann.
Das gesamte Programm ist in MATLAB Ver.4.2c (1994) geschrieben
und somit lauffähig auf fast allen MATLAB Versionen. Eine
Umschreibung des Programms auf MATHCAD ist schnell möglich.
Bestechend ist die Einfachheit so einer Vorgehensweise. Im Programm
wird die Funktion ode45 benutzt, deren engl. Hilfe folgt aussieht:
ODE45 Solve differential equations, higher order method.
ODE45 integrates
a system of ordinary differential equations using 4th and 5th
order Runge-Kutta formulas.
[T,Y] = ODE45('yprime', T0, Tfinal, Y0) integrates the
system of ordinary differential equations described by the M-file
YPRIME.M, over the interval T0 to Tfinal, with initial conditions
Y0.
Um unser Phygoidenproblem zu beschreiben,
brauchen wir vier Zustandsvariablen: x, vx, y, vy; was in dem
Programm als z[1:4] Vektor wiederzufinden ist. Als Anfangswerte
z0 setze ich folgendes: x=y=0 (Anfang des Koordinatensystems).
Vx und vy sind die Komponenten der Startgeschwindigkeit. Die
ganze Lösung wird dann in einer Zeile berechnet:
[t,z]=ode45('phypri',0,50,z0)
wobei die 'phypri' die Ableitungen der 4 ordinaren Gleichungen
bereitstellt. Die Zeilen mit figure(); plot() generieren dann
die 4 Bilder.
Das Programm:
================= snipp =======================
% phygo
% Dieses Programm simuliert die Phygoide nach den Gleichungen
% die aus der Gleitzahl E=1/e und der mittleren Geschwindigkeit
w
% das Verhalten nach einer Störung beschreiben.
% Parameter:
% e= 1/E 1/Gleitzahl
% w= Mittlere Geschwindigkeit
% g= Erdbeschleunigung
global w e g
% folgend Beispielwerte
e=0.063; % 1/Gleitzahl
w=11.7; % Grundgeschwindigkeit
m/s
g=9.81;
% Erdbeschleunigung g=9.81m/s^2
z0=[0 0 4.91 0.95]; % Startvektor V=5m/s viel zu
langsam, Vx=4.91m/s Vy=0.95
[t,z]=ode45('phypri',0,50,z0); % die Lösung
figure(1);clf;plot(t,sqrt(z(:,3).^2+z(:,4).^2));grid;xlabel('Sek');ylabel('m/s');
title('Fluggeschwindigkeit');set(gca,'xtick',[0:5:100]);
figure(2);clf;plot(t,atan2(z(:,4),z(:,3))*180/pi);grid;xlabel('Sek');ylabel('Grad');
title('Bahnwinkel');set(gca,'xtick',[0:5:100]);
figure(3);clf;plot(t,-z(:,2));grid;xlabel('Sek');ylabel('Höhe
[m]');
title('Höhe');set(gca,'xtick',[0:5:100]);
figure(4);clf;plot(z(:,1),-z(:,2));grid;xlabel('x [m]');ylabel('y
[m]');
title('Flugort');
================= snipp =======================
Die Funktion 'phypri'
================= snipp =======================
function zp=phypri(t,z)
% phygoide prim aus Phygoide/Helmut Schenk Seite 12
% z(1)==x
% z(2)==x'
% z(3)==y
% z(4)==y'
global w e g % grundgeschw. w,
e=1/E Gleitzahl, g=9.81m/s^2
zp(1)=z(3);
zp(2)=z(4);
zp(3)=g*sqrt(z(3)^2+z(4)^2)/w^2*(z(4)-e*z(3));
% Gleichung (25)
zp(4)=g*(1-sqrt(z(3)^2+z(4)^2)/w^2*(z(3)+e*z(4))); % Gleichung
(26)
================= snipp =======================
Beispiel:
Ein Modellflugzeug mit einer Grundgeschwindigkeit von 11.7m/s und einer Gleitzahl von 15.87 (e=1/15.87=0.063) wird viel zu langsam mit v=5m/s (vx=4.91 vy=0.95) leicht nach unten (alpha=11 Grad) freigegeben. Seine Fluggeschwindigkeit
wächst schnell von 5m/s auf über 16m/s, indem es den Bahnwinkel fast bis 45 Grad vergrößert.
Es verliert rasant über 15m Höhe.
Nach ca. 2,7 Sekunden hat es aber so viel kinetische Energie aus der potentiellen Energie der Höhe aufgenommen, daß es die Überfahrt (16m/s anstatt 11.7m/s) in Höhe zurück umsetzen muß. Der Bahnwinkel wird negativ, das Flugzeug steigt bis zum wiederholten lokalen Maximum an Höhe und Minimum an Geschwindigkeit (unter 8m/s bei ca. 5.4 Sekunden). Dieses Spiel wiederholt sich mit einem Dämpfungsfaktor, der von der Gleitzahl (cw und ca. eingeschlossen) abhängt. Der Flugort ist im folgenden sichtbar:
Was lernt man daraus:
Diese Simulation verrät, daß man entweder die Fluggeschwindigkeit und/oder den Bahnwinkel zu Steuerung des Höhenruders heranziehen könnte, um den Wellenflug zu "glätten". Diese Vermutung kann man mit dem vorgestellten Programm sofort simulieren, in dem man die Funktion 'phypri' gegen 'phypri2' austauscht.
================= snipp =======================
function zp=phypri(t,z)
% phygoide prim aus Phygoide/Helmut Schenk Seite 12
% z(1)==x
% z(2)==y
% z(3)==x'
% z(4)==y'
global w e g k2% grundgeschw. w, e=1/E Gleitzahl, g=9.81m/s^2
zp(1)=z(3);
zp(2)=z(4);
zp(3)=g*sqrt(z(3)^2+z(4)^2)/w^2*(z(4)-e*z(3));
% Gleichung (25)
for i=1:3;
w=11.7-k2*zp(3);
zp(3)=g*sqrt(z(3)^2+z(4)^2)/w^2*(z(4)-e*z(3));
% Gleichung (25)
end;
zp(4)=g*(1-sqrt(z(3)^2+z(4)^2)/w^2*(z(3)+e*z(4))); % Gleichung
(26)
================= snipp =======================
In dieser Funktion kommt die Aufschaltung der Fluggeschwindigkeit
(vx) mit einem Faktor k2=0.1, 0.2, 0.4 auf die Grundgeschwindigkeit
des Flugzeug in folgender Form:
w=11.7-k2*vx
Leider ist der Zustand der dadurch perturbierenden Variablen
nicht gleich ausgeglichen, deshalb mache ich eine zusätzliche
Iteration von 3 für jeden Zeitpunkt.
Ergebnis der Simulation der Wellenkompensation.
Das Ergebnis ist positiv. Durch die Variation der Grundgeschwindigkeit w durch die Fluggeschwindigkeit kann man die Abklingzeit der Phygoide nach einer Bahnstörung wesentlich verkürzen. Die Grundgeschwindigkeit eines Modellflugzeugs kann durch das Höhenruder leicht beeinflußt werden. Wir tun es durch die Trimmung des Ruders. Für die Thermik wählen wir niedrige Fluggeschwindigkeit mit einem hohen Ca-Wert, für den Hangflug oder Streckenflug trimmen wir "tiefer", um eine höhere Fluggeschwindigkeit mit niedrigeren Ca-Wert zu erreichen. Das gleiche kann der Flugregler tun, um den Bahnwinkel und die Geschwindigkeit zu stabilisieren.
Durch diese Aufschaltung stabilisiert sich die Fluggeschwindigkeit schon nach der kurzen Zeit von 10 Sek., wo ohne Aufschaltung sogar nach 50 Sek. noch keine konstante Fluggeschwindigkeit vorhanden war.
Auch der Bahnwinkel wird schnell konstant
und die Höhenkurve gleicht jetzt einer Geraden:
Das Flugzeug fliegt mit beruhigter Phygoide.
Eine ähnliche proportionale Aufschaltung des Bahnwinkels
führt zu fast identischen Ergebnissen, deshalb verzichte
ich hier auf die Präsentation der Bilder. Jeder kann also
entscheiden, welchen Sensor er lieber nimmt, um den Wellenflug
unserer Modelle zu verkleinern.
Flugparametermessungen an echten Modellflugzeugen.
Das Ergebnis der Simulation veranlaßte mich, auch in meinem 2.5m Flugzeug beide Sensoren einzubauen:
- für die Fluggeschwindigkeit ein selbstgebasteltes Flügelrad aus GfK, montiert in einem 12 cm hohen CfK-U-Rahmen (Masse zusammen 1.4g), mit einem Hallsensor (von Conrad-Elektronik) auf der Achse;
- für den Bahnwinkel einen schnellen (bis 4000Meß/Sek) und genauen (0.1 Grad/100ms) Neigungssensor (Nicksensor) basierend auf dem ADXL05 von Analog Devices.
Um die Höhe zu beurteilen, kam noch
ein präziser, elektronischer Höhenmesser eigener Konstruktion
hinzu. Er basiert nicht auf einem günstigeren, monolithischen
Drucktransducer (z.B. Motorola, Sensym, etc.), sondern auf einem
genaueren Hybrid-Dickschicht-Transducer, komplett Laser-abgeglichen,
temperaturstabilisiert und mit strombestimmenden Elementen auf
dem Hybrid versehen. Das Signal aus dem Sensor wird mit einem
ausgezeichneten, extrem rauscharmen Instrumenten-Verstärker
AD620 von Analog Devices verarbeitet. Die ganze Höhenmeter-Einheit
ist ca. 14x14x7mm groß. Das rms-Rauschen auf meinem Labortisch
betragt ca. 2cm Höhe in der Sekunde. Im Flugzeug ist der
Einbauort schon wichtig, will man genaue Messungen durchziehen.
Die Flugaufzeichnungen wurden mittels 14bit 8-kanaligen Telemetrie
gemacht. Um später zu wissen, wie das Flugzeug gegen den
Wind stand, benutzte ich einen elektronischen Kompaß (Vector2x
von Precision Navigation/Kalifornien). Die Anzahl der Datenpunkte,
obwohl 32/Sek gemessen, müßte auch reduziert werden,
um mit handlichen Datensätzen von 3000-6000 Datenpunkten
operieren zu können. 16 Meßpunkte werden zu einem
0.5-sekundigen Wert zusammengefaßt und per Funk zu "Erdstation"
übermittelt. Ich kann auch 64/Sek. Messungen einordnen,
aber die 0.25-sekundigen Werte bringen wegen der 5- bis 6-sekundigen
Phygoide keine Qualitätsverbesserung. Der enorme Vorteil
der riesigen Datenmengen liegt darin, daß ich in beliebige
Abschnitte eines 30-Minutigen Fluges mehrfach hineinzoomen kann,
ohne auf die Gesamtauflösung verzichten zu müssen.
Thermikflug im Winter?
Wie so eine Höhenaufnahme eines einstündigen Thermikflugs, gestartet um 16:10h (beendet fast in Dämmerung um 17:10h) im Winter auf einer Wiese aussieht, zeigt das folgende Bild:
(Nellinger Siedlung sollte Nelliger Siedlung Straße heißen, sorry). Unten sieht man, wann der Motor eingeschaltet war (grüne Linie). Die ersten 1100 Sekunden hatte ich in 6 Aufstiege, mit ca. 2.5m/s Steigen, kein Glück. Erst der siebte Aufstieg mit fast leeren Batterien (8 Zellen) brachte mich in die mühsame Thermik mit Steigen von 0.16m/s über 400 Sekunden lang, aber dann mit 0.58m/s bis auf die Höhe von 300m, wo ich den Bart verlassen habe und im Gleitflug (ca. 1600 Sek) zurück zu meinem Standort geflogen bin. Dort konnte ich mich noch ca. 500 Sek. im Nullschieber halten, bevor es mit 0.56m/s abwärts ging. Die Aufstiege 8 und 9 brachten kein Glück mehr. Hier sehen Sie also 6000 1/2-sekündliche Meßpunkte nur der Höhe. Dazu gehören natürlich 6000 Nasenneigungs-, Querneigungs-, Geschwindigkeits- und Kompaßkurswerte, die dann zusammen in kürzeren Abschnitten analysiert werden können. Das nächste Bild
zeigt nur die Thermikflugphase eines Mittagsfluges um 13:25h.
Man ahnt schon, daß dort ein starker Wellenflug mit einer
Periode von 5 Sek. vorhanden ist. Die Steigraten wieder kaum
sichtbar: 0.3m/s oder 0.11m/s.
Die Phygoide: ist sie meßbar an unseren Modellen ?
Um eine einigermaßen ungestörte Phygoide zu sehen, muß man einen ruhigen Gleitflug wählen. Natürlich ist die Phygoide auch bei allen Thermikflügen vorhanden, aber dann überlagert mit den kämpferischen Knüppelbewegungen des Piloten (Knüppelthermik).
Hier ein Sonnenuntergangsflug im Winter mit vier Steigphasen (Steigen 1.87m/s, 1.68m/s, 1.4m/s, 1.28m/s), wo es nur "rauf und runter geht", wird der erste Sinkflug im folgenden Bild analysiert:
Es ist also ein ruhiger Flug. Man erkennt nach der Beendigung der Steigphase, die eine Bahnstörung hervorruft, den Wellenflug, die Phygoide. Das nächste Bild zeigt die zugehörigen Neigungsoszillationen:
Die ersten 7 Oszillationen haben eine mittlere Periode von 4.84 Sek. mit einer Standardabweichung von 0.34 Sek., also 7% des Wertes.
Um die beiden Kurven für die Eigenschaft einer Phygoide zu vergleichen, will ich alle Signale entfernen, die wesentlich langsamer als die Phygoide selbst sind. Ich entferne alle Signale mit einer Periode größer als 10 Sek. Für diesen Zweck entwerfe ich ein digitales Butterworth-Hochpaßfilter 4. Ordnung mit einer Grenzfrequenz von 0.1Hz (T=10Sek) und der Form:
y(n) = b(1)*x(n)
+ b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
- a(2)*y(n-1) - ... - a(na+1)*y(n-na)
Die Koefizienten sind:
b = 0.6620 -2.6481 3.9721
-2.6481 0.6620
a = 1.0000 -3.1806 3.8612
-2.1122 0.4383
Um keine Phasenverschiebung auf die Kurven zu bringen, wende ich diesen Filter vorwärts und rückwärts in der Zeit. Damit ist die Operation gleich einem rekursiven Filter 8. Ordnung (IIR) ohne Zeit- und Phasenverschiebung. Das Ergebnis zeigt für den stabilen Sinkflug ab der 440. Sek. sehr gute Übereinstimmung zwischen Höhenkurve (gelb) und Neigungskurve (rot):
Beim Übergang von der Steigphase in den Sinkflug treten
starke Winkeländerungen (plötzlich von +18 Grad auf
0) auf, so daß die Neigungsdifferenz bei 435 Sek. natürlich
zusätzlich einen Hochpaß-gefilterten Sprung von 18
Grad auf Null zeigt. Die Amplituden der Schwingungen korrelieren
nicht sehr gut, das muß noch geklärt werden. Grundsätzlich
gilt jedoch, daß die Neigung als Aufschaltgröße
für die Beseitigung des Wellenfluges tatsächlich geeignet
wäre.
What next?..... Was kommt demnächst ?
Um die Möglichkeit der Beseitigung des Wellenflugs zu verifizieren, braucht man im Flugzeug eine aktive Steuerung durch einen PD-Regler des Höhenruders und telemetrische Aufnahmen einer Bahnstörung ohne und mit dem Regler in der Steuerschleife. Ich habe zwar die Simulation für einen Proportional-Regler in der Form w=A-k2*vx gezeigt, aber ein um den differentialen Teil erweiterter Regler eignet sich besser (Simulation und Ergebnisse vorhanden).
Dafür baue ich eine aktive Steuerung meines 2.5m-Modells. In diesem Falle greift man in die Sicherheit des Fluges ein. Deshalb muß die Steuerung durch mehrere Maßnahmen gesichert sein. Ich wende softwaremäßige Maßnahmen (Safety Guards) an, die alle beteiligten Programmsegmente laufend überprüfen und die Variablen auf deren Bereich (AD-Umsetzer der Sensoren liefert falsche Werte) und deren Änderungen testen. Erst nach Freigabe durch unabhängige Programmwächter werden die Meßgrößen in die Rechnung aufgenommen und die Steuersignale für die Servos generiert. Jede kleinste Fehlermeldung des internen Überwachungssystems schaltet die Steuerung auf den Piloten um. Das sind die Softwaremaßnahmen.
Um den hardwaremäßigen Ausfall (Oszillator steht,
oder der Prozessor spinnt) auch zu verhindern, wacht ein Watchdog
(Wachhund) über den Schalter (Multiplexer), der die Servokabel
von dem Empfänger (Pico2000) auf das Mikroprozessorsystem
umschaltet. Erst der Wachhund kann, nach entsprechender "Bitte"
des Mikroprozessorsystems, die Steuerung des Flugzeugs
dem System übergeben. Natürlich kann der Pilot per
RC-Fernbedienung jederzeit und ohne jegliche Schalterbetätigung
die Steuerung übernehmen. Diese Maßnahmen verlangen
einige Arbeit, wovon jedoch ein großer Teil schon fertig
vorliegt. Deshalb ist es nicht sehr weit, bis ich die Ergebnisse
einer aktiven und autonomen Beseitigung des Wellenflugs präsentiere.
Zuerst jedoch steht ein zweiwöchiger Urlaub "in warmen
Gefilden" bevor.