zuletzt geändert: 22.03.98

Die Phygoide: Wellenflug-Simulation und Wirklichkeit.

von Peter Rother, Feb. 1998

Die Simulation

Das Programm
Beispiel
Was lernt man daraus?
Ergebnis der Simulation der Wellenkompensation

Flugparametermessungen an echten Modellflugzeugen

Thermikflug im Winter ?

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.
 
 

Allen Interessierten, Nachbauern und "Weiterbauern" viel Glück!

Peter Rother, am 24.2.1998. Fasching !! Muß raus.